cannin/biostars_qa huggingface dataset:
Question: Hi @ all,
I have cel files from 2 groups, 3 replicates each group. I want to use limma for the differential gene expression analysis.
My script looks like this

prefix = "101053-0"
cels <- ReadAffy ( filenames = paste ( paste (prefix, c(4,5,6,1,2,3), sep=""), ".CEL", sep="" ),
sampleNames = c("N402_rep_1", "N402_rep_2", "N402_rep_3", "B36_rep_1", "B36_rep_2", "B36_rep_3") )

results <- expresso (cels, bgcorrect.method = "rma", normalize.method = "quantiles", pmcorrect.method = "pmonly", summary.method="medianpolish")
design <- model.matrix(~factor(c(1,1,1,2,2,2)))
fit <- lmFit (results, design=design)
ebayes <- eBayes (fit)
topTable (ebayes)

The first gene of the list is this

logFC AveExpr t P.Value adj.P.Val B
An00g11929_at -7.347305 8.465857 -110.96206 6.249263e-18 9.095178e-14 27.44659

So logFC is reported as -7.34
But when I look at the absolute expression values

> exprs(results)["An00g11929_at", ]
N402_rep_1 N402_rep_2 N402_rep_3 B36_rep_1 B36_rep_2 B36_rep_3
12.127017 12.148645 12.142866 4.683787 4.842077 4.850748
The logFC is more about -1.5 than -7.34
How can this be?

I read that limma expect log2-transformed values. I also read that using median polish to summarize already produces log2 transformed expression values.
Anyway, I transformed the value by my self

exprs(results) <- log2(exprs(results))

And now the logFC is what I would expect from the un-transformed values

logFC AveExpr t P.Value adj.P.Val B
An00g11929_at -1.3411349 2.931070 -80.01571 2.751338e-12 4.004298e-08 16.27317

But now the log2-transformed values do not fit to the logFC

N402_rep_1 N402_rep_2 N402_rep_3 B36_rep_1 B36_rep_2 B36_rep_3
3.600153 3.602724 3.602037 2.227675 2.275626 2.278207

So I'm very confused. Why do the absolute expression values do not match the logFC? Do I have to apply the log2-transformation on my own?

Would be very nice if someone could help :)

Answer: The data stored in `exprs` is already log-transformed. The intensities on the array will be in the thousands.

A _back-of-the envelope_ calculation is as follows:

estimate of fold-change from the de-transformed array data:

geometric_mean(array_data(group_2)) / geometric_mean(array_data(group_1))
~ typical_value(group_2) / typical_value(group_1)
~ (2 ^ log_transformed_typical_value (group_2)) / (2 ^ log_transformed_typical_value(group_1))
~ (2^5) / (2^12)
~ 0.008
~ 1 / 125

So you get (an inaccurately calculated) ~ 125-fold change

estimate of log-fold-change from the transformed array data (as stored in exprs)

arithmetic_mean(exprs(group_2)) - arithmetic_mean(exprs(group_1))
~ typical_value(group_2) - typcial_value(group_1)
~ 5 - 12
~ -7

So your estimate of fold-change is `2^-7 ~ 0.008`; again 1/125
Question: (1) I have uploaded the genome fasta and GFF annotation in the mysql data using following command of `Bio::DB::SeqFeature::Store` database

./bp_seqfeature_load.pl -c -d WS240 -u xxx -p xxx c_elegans.PRJNA13758.WS240.genomic.fa c_elegans.PRJNA13758.WS240.annotations.gff3

(2) Following tables are inserted into my database WS240

```
+-----------------+
| Tables_in_WS240 |
+-----------------+
| attribute |
| attributelist |
| feature |
| locationlist |
| meta |
| name |
| parent2child |
| sequence |
| typelist |
+-----------------+
```

(3) Now I want to make query in the mysql to extract the sequence and its feature using following commands:

```
#!/usr/bin/perl
use strict;

use Bio::DB::GFF;
my $db = Bio::DB::GFF->new(-dsn => 'dbi:mysql:database=WS240',
-user => 'XXX',
-pass => 'XXX',
-aggregators => 'gene_model{coding_exon,5_UTR,3_UTR/CDS}');

my $gene_stream = $db->get_seq_stream('gene_model:curated');

while (my $gene = $gene_stream->next_seq) {
print $gene->name,"\n";
for my $part ($gene->get_SeqFeatures) {
print "\t",join("\t",$part->method,$part->start,$part->end),"\n";
}
print "\n";
}
```

(4) I am getting following error:

```
------------- EXCEPTION: Bio::Root::Exception -------------
MSG: Couldn't execute query SELECT fref,fstart,fstop,fsource,fmethod,fscore,fstrand,fphase,gclass,gname,ftarget_start,ftarget_stop,fdata.fid,fdata.gid
FROM fdata,ftype,fgroup
WHERE fgroup.gid = fdata.gid
AND ftype.ftypeid = fdata.ftypeid
AND ((fmethod = ? AND fsource = ?) OR (fmethod = ? AND fsource = ?) OR (fmethod = ? AND fsource = ?) OR (fmethod = ? AND fsource = ?))
ORDER BY fgroup.gname:
Table 'WS240.fdata' doesn't exist

STACK: Error::throw
STACK: Bio::Root::Root::throw /usr/local/share/perl/5.10.0/Bio/Root/Root.pm:368
STACK: Bio::DB::GFF::Adaptor::dbi::caching_handle::do_query /usr/local/share/perl/5.10.0/Bio/DB/GFF/Adaptor/dbi/caching_handle.pm:123
STACK: Bio::DB::GFF::Adaptor::dbi::range_query /usr/local/share/perl/5.10.0/Bio/DB/GFF/Adaptor/dbi.pm:627
STACK: Bio::DB::GFF::Adaptor::dbi::get_features_iterator /usr/local/share/perl/5.10.0/Bio/DB/GFF/Adaptor/dbi.pm:1008
STACK: Bio::DB::GFF::_features /usr/local/share/perl/5.10.0/Bio/DB/GFF.pm:3476
STACK: Bio::DB::GFF::features /usr/local/share/perl/5.10.0/Bio/DB/GFF.pm:1091
STACK: Bio::DB::GFF::get_seq_stream /usr/local/share/perl/5.10.0/Bio/DB/GFF.pm:1132
STACK: ./aa.pl:12
```

(5) Could you kindly suggest me what is wrong with the script (3), what is the best approach to make query with mysql. I also used the `bp_bulk_load_gff.pl` ([Bio::DB::GFF][1])) to upload the GFF file into mysql but not working well.

Thanks
Firoz

[1]: http://wiki.wormbase.org/index.php/Mining_WormBase_with_Bio::DB::GFF

Answer: See this thread and the answer by Timothy Parnell: http://generic-model-organism-system-database.450254.n5.nabble.com/Getting-gbrowse-syn-to-work-with-MySQL-td5519255.html

You should use `bp_bulk_load_gff.pl` with `Bio::DB::GFF` and `bp_seqfeature_load.pl` with `Bio::DB::SeqFeature::Store`. However your script uses `Bio::DB::GFF` with `bp_seqfeature_load.pl`. You need to change one of them to conform.
Question: Hi,

Making some deepTools heat maps and wondering how I can move the colorbar to the bottom of the image or is this only possible for the multiple color heatmaps shown on the deepTools site? Thanks.

Rob.

Answer: This is only possible if you specify multiple colormaps or scales (e.g., `--colorMap Spectral summer` or `--yMin 0 1 2`). Truth be told, you should be able to specify the same thing multiple times, such as `--yMin 0 0 0` and get the same effect. Note however that I really wouldn't recommend that, as it both busies up the resulting plot and makes it less apparent that the scales are identical everywhere.
Question: In RPPA & MASS-SPEC, which procedure gives greater discovery of proteins in a given sample or which one has higher sensitivity. What is the difference between CPTAC Mass spec data and TCGA protein expression data.?

Answer: According to these authors, RPPA is better than MS:

https://www.google.ru/url?sa=t&rct=j&q=&esrc=s&source=web&cd=1&ved=0ahUKEwiLtcCduJHMAhViJJoKHcQqBe4QFggdMAA&url=http%3A%2F%2Fwww.mdpi.com%2F2076-3905%2F4%2F2%2F98%2Fpdf&usg=AFQjCNEGF0eX_L6LNJ67mkrQPlzoPuHD9w&bvm=bv.119745492,d.bGs&cad=rjt

A few citations from the article:

“Reverse Phase Protein Arrays (RPPA) represent a very promising sensitive and precise high-throughput technology for the quantitative measurement of hundreds of signaling proteins in biological and clinical samples. This array format allows quantification of one protein or phosphoprotein in multiple samples under the same experimental conditions at the same time”.

“Mass spectrometry (MS)-based technologies have rapidly advanced in recent years. Beside other advantages mentioned in Table 1, MS-based methods enable the identification of new biomarker candidates by comparing protein signals obtained from cancer and healthy tissues (de novo discovery platform) [7–10].
However, in author opinion, MS-based technologies are not suitable for the use in clinical routine at the moment. First, due to the complex sample preparation and secondly due to the insufficient profiling of low-abundance signaling proteins.”

Below people discuss TCGA and CPTAC:

https://followthedata.wordpress.com/tag/tcga/

"The Cancer Genome Atlas (TCGA) has become a focal point for a lot of genomics and bioinformatics research. DNA and RNA level data on different tumor types are now used in countless papers to test computational methods and to learn more about hallmarks of different types of cancer.

Perhaps, though, there aren’t as many people who are using the quantitative proteomics data hosted by Clinical Proteomic Tumor Analysis Consortium (CPTAC). There are mass spectrometry based expression measurements for many different types of tumor available at their Data Portal."

The CPTAC Data Portal from above:

https://cptac-data-portal.georgetown.edu/cptacPublic/

Below there is more general discussion and there are more disease types:

http://www.ncbi.nlm.nih.gov/projects/gap/cgi-bin/study.cgi?study_id=phs000892.v2.p1

It's difficult to say what is more important: one approach works with data submitted by another one.

"Renamed the Clinical Proteomic Tumor Analysis Consortium, CPTAC is beginning to leverage its analytical outputs from Phase I to define cancer proteomes on genomically-characterized biospecimens.

The purpose of this integrative approach is to provide the broad scientific community with knowledge that links genotype to proteotype and ultimately phenotype. The data contained in this dataset are derived from samples designed to confirm CPTAC findings from the TCGA samples."

So the TCGA provide samples. CPTAC study them and insert the information into the database.
Question: Hi all

I am new to Biostars and blast.
I am trying to convert a (test)-database of 25 sequences into a blast-database using standard tutorial commands. However, I get the following error:

*Error: mdb_env_open: There is not enough space on the disk.*

In my cd I get a file of >200gb which seems ridiculously high for only 25 sequences. Does anyone have an idea what I am doing wrong?

Kind regards

Michaël

Answer: Response from NLM support solved the issues. Here it is.

Thank you for the report. This is a known issue with the Windows release. The program makeblastdb attempts to allocate a very large amount of virtual memory. You can solve the problem by setting the a new BLAST environment variable
BLASTDB_LMDB_MAP_SIZE=1000000
See the BLAST setup documentation for details on how to set Windows environment variables (https://www.ncbi.nlm.nih.gov/books/NBK52637/). Once you change the variable, you'll need to close and reopen the command window where you were running BLAST for the new setting to take effect.
Question: Hi,

Does anyone know if you should remove the genes that fail independent filtering in deseq2 (I.e have NA padj values). I've read discussion around its purpose for deseq2 but I'm not sure how these genes should be treated when performing GSEA.

Removal of these genes before running fGSEA tends to result in finding more pathways (and higher padj values) which aligns with the idea of having higher power. However I want to make sure using this approach is valid.

Thank you.

Answer: Yes. Independent filtering is a methodology to improve F1 score by removing genes with low variance, mean expression, dispersion, or other noise surrogate. Genes that are removed by IF should be regarded as "not effectively profiled" and ignored in downstream analysis.
Question: I have seen cases that one PDB can have multiple BioUnit files, and I am not quite the sure the reason they make multiple biounits. Some biounits only contain the binding peptide. Could I just use the first biounit file? Or I need to be careful about which one to choose.

Because I want to do it in large scale, rules like "choosing the one with the largest number of residues" would be better.

Answer: There's 2 main reasons for having multiple biounits in the PDB:

1. When there are multiple equivalent copies of the biounit in the crystal, then they are annotated as several biounits. A very typical example would be 4re6: a dimeric protein with 4 chains in the asymmetric unit. In that example biounit 1 is the dimer between chains A and B, biounit2 the dimer between chains C and D. Both are equivalent, but as it appears twice in the crystal it is given as 2 biounits.
2. Another reason is that the PDB tries to accommodate different opinions about them. The annotations can come from the authors or from software predictions (mainly PISA). So many times you will encounter both of them: biounit1 authors and biounit2 PISA. It can also happen that a certain biounit is annotated by both PISA and authors, when both are agreeing. The same 4re6 example above is also showing this, the first 2 biounits are from both authors and PISA, whilst biounit3 is only PISA (predicting a tetramer). As you see even PISA can have multiple predictions and in some cases like 4re6 they are both added to the PDB.

An important point to make is that experimentally coming out with a correct biounit is not a simple task. From crystallographic data alone you can't really do much, you need to use other experimental methods like gel filtration, analytical ultra centrifugation, light scattering, mutagenesis etc in order to be sure of the oligomeric state in solution of the protein. Those methods can sometimes not be conclusive enough, thus making things complicated.

As to which of the biounits to use there's not a simple straight forward answer to that. One very accepted method is to simply use biounit1 and ignore the rest. That is counting on that the PDB annotators often will set biounit1 to their main preference (due to good authors data or software having confident predictions). Another method would be to use the first biounit annotation coming from authors, trusting that authors are doing their job correctly.

A last word of caution comes with errors in annotations: there are quite a few errors in biounit annotations in the PDB, see for instance our paper where we analyse the problem in some detail.
Question: I have paired end sequence data which I have mapped to a reference. I need to extract reads from the bam file that map either next to each other or within a defined interval from each other, for example, between 200 and 300 bp from each other. Is it possible to do this?
For instance, if my first read pair map between position 90,000 and 90,130 and a second read pair maps between 90,380 and 90,500. These 2 reads map 250 bp apart on the reference. Is it possible to extract all the reads that have this characteristic?
Thanks

Answer: Sure, you can do this with a few genomic set operations with [BEDOPS][1] tools:

1. Convert your reads to sorted BED with `bam2bed`
2. Use `bedmap --range` to find elements within the first, smaller range boundary (set `A`)
3. Use `bedmap --range` again to find elements within the second, larger range boundary (set `B`)
4. Use `bedops --not-element-of` to filter out elements from set `B` that overlap elements in set `A`

What's left at the end are your elements between the ranges specified for sets `A` and `B`.

For instance, to get reads that are between 200 and 300 bases from one another:

$ bam2bed < reads.bam > reads.bed
$ bedmap --range 200 --count --echo-map --delim '\t' reads.bed | awk '$1>1' | cut -f2- | tr ';' '\n' | sort-bed - > reads_within_200_bases.bed
$ bedmap --range 300 --count --echo-map --delim '\t' reads.bed | awk '$1>1' | cut -f2- | tr ';' '\n' | sort-bed - > reads_within_300_bases.bed
$ bedops --not-element-of 1 reads_within_300_bases.bed reads_within_200_bases.bed > answer.bed

The file `answer.bed` has reads between 200 and 300 bases of one another.

[1]: https://github.com/bedops/bedops
Question: Hi everyone. I would like to ask some doubts about the representation of data of FastQC. I have two cleaned files 1.fq and 2.fq of Illumina. When I make a analysis with FastQC. I saw that first file does not have any contigs with Phred score lower of 20. In contrast the second file shows a lot of contigs with lower Phred score of 20. However, when I merged both files into one. This unique file does not show a Phred score low of 20. I attached the images. Somebody can me explain the why of this? is it trustful?

Regards.

**File 1**

**File 2**

**File All merged**

[][1]

[1]: http://www.freeuploadimages.org/images/zzk2gchqlmqc463mfgza.jpg

Answer: simply, in these boxplots the lower whisker represents the 10th percentile, which means that 10% of the points are lower than that and are not displayed. when you merge two datasets the 10% will be higher, exactly as you see in your data. there are still quite a considerable number of sequences which have bases with a lower quality than 20 at the end.
Question: Hi,

I have a set of samples and log10 total read counts of them in a data frame like below

> head(a)
log10_total_counts
A1 6.468503
A10 6.565213
A11 6.752139
A12 5.078598
A2 6.277342
A3 6.473411
>

![enter image description here][1]

[1]: https://i.ibb.co/P6mmKc5/Picture1.png

For instance red dots have less than 1000000 counts and black dots have more thsn 1000000 counts.

How could I repreduce this plot please?

Answer: Something like the following would work with base R functions:

# make some data
a <- data.frame(sample=paste0("A", 1:50), log10_total_counts=rnorm(50, 6))

# get x points
x <- 1:length(a$log10_total_counts)

# create plot, drop x axis
plot(x, a$log10_total_counts, ylim=c(0,8), pch=19, xaxt="n", xlab="sample")
# add dashed line and X-axis labels with rotation
abline(h=6, lty="dashed")
axis(1,1:50, labels=a$sample, las=3)

# select points to color differently
LessThan6 <- a$log10_total_counts < 6
points(x[LessThan6], a$log10_total_counts[LessThan6], pch=19, col="red")
Question: Hello everyone!

I am learning to create a volcano plot by using EnhancedVolcano in R.
But I am struggling to color custom points (genes).

**Here is my input file:**

library('EnhancedVolcano')

res <- read.table("text.txt", header=TRUE)
res$gene <- as.character(res$gene)

res object contains:

head(res)
# gene log2foldchange pvalue padj
#1 Cx3cr1 -8.039239 1.29e-118 2.780000e-116
#2 Trem1 -5.258502 3.01e-44 1.770000e-42
#3 Serpina3f 3.202818 2.12e-09 2.400000e-08

Also, I have gene lists which I want to label(able to do it), and color (cannot do it) in "res$gene" object (overlapping):

gene_list <- scan("503-5ptarget.txt", what="", sep="\n")

[1] "gene" "Fam122A" "Ccnd2" "Usp2" "Arl2" "Rab9B" "Dcaf7" "Pom121" "Ccdc42B" "N4Bp1" "Tmem74B" "Akt3"

Now I created this volcano plot by injecting "gene_list" into "res$gene". But the problem was I can "label" overlapped gene, But I "cannot color that point".

aa2 <- EnhancedVolcano(res,
lab = res$gene,
x = 'log2foldchange',
y = 'pvalue',
title = 'test',
subtitle = "test",
pCutoff = FALSE,
FCcutoff = FALSE,
xlim = c(-10, 10),
pointSize = 2.0,
cutoffLineType = 'blank',
selectLab = gene_list,
labCol = 'black',
labFace = 'bold',
colAlpha = 1,
shade = gene_list,
shadeLabel = 'gene list 1',
shadeAlpha = 1/2,
shadeFill = 'red',
shadeSize = 1,
shadeBins = 5,
col = c('grey', 'grey','grey', 'grey'),
drawConnectors = TRUE,
gridlines.major = FALSE,
gridlines.minor = FALSE,
legendVisible = FALSE)

Here is the plot:

[][1]

How can I color these overlapped gene lists? thank you

[1]: http://%3Ca%20href=%22https://ibb.co/92snqsC%22%3E%3Cimg%20src=%22https://i.ibb.co/92snqsC/Rplot.png%22%20alt=%22Rplot%22%20border=%220%22%3E%3C/a%3E

Answer: Hey,

If you want to label specific genes, then you just need to use the `selectLab` parameter, as elaborated here:

- 4.7
Only label key variables

If you want to customise the colouring of the points, then you need to follow:

- 4.9
Over-ride colouring scheme with custom key-value pairs

A much easier way to colour-mark key genes is by following this section (but this is a relatively new feature and is not fully developed to my liking):

- 4.11
Shade certain variables

Kevin
Question: I have a transcriptome data of 10.8gb R1.fastq and R2.fastq each. I then compressed these pairs using `gzip R1.fastq` and `gzip R2.fast2`, and now the files are 2.2gb and 2.4gb. Is it possible for two compressed files to have different size when the uncompressed files are of same size?

Answer: Yes. It's perfectly possible, even if the reads are the same length. One might have sequences that are a little more repetitive, and therefore more compressible. If they have the same number of lines, that's all that matters.

It of course also possible to run gzip with different levels of compression, but you don't seem to have done that. in this case.
Question: Hi

I ask a question about term "downstream".

I know that downstream refers to the direction of 3'. That is all I have learned.

While reading below sentence, however, I am so confused that why the term "downstream" was used.

> After mapping reads to the reference genome, a three step post-alignment processing procedure is recommended to minimize the artifacts that may affect the quality of **downstream **variant calling. It consists of read duplicate removal, indel realignment, and base quality score recalibration.

> Base quality is a critical factor for accurate variant detection in the **downstream **analysis.

I fully understood above sentence except the downstream.

I cannot connect downstream I knew with the potential meaning in the sentence.

Could you help me why the term downstream was used for those sentence?

Answer: You might just want to look "downstream" and "upstream" up in a dictionary. These are generic English terms. But because I'm feeling nice this morning:

Imagine a river. Rivers flow in a single direction and if you placed something (that floats) in a river it would float with the current, passing by points along the coast in a set order (after all, it's not going to float against the current). If you were standing on the shore and the object were floating toward you, then it'd be upstream of you. After it passed you, it'd be downstream of you. In English, we generalize these terms to apply to any set of processes with a (roughly) linear ordering. So when you first get sequencing data, you'll probably perform some quality control steps and then alignment or assembly. These are the upstream-most steps. Once, that's done, the downstream steps (variant calling, counting alignments per gene in RNAseq, calling ChIP peaks, etc.) are performed, which generally is when the data starts getting more interesting.
Question: I'm using PhyML to compute phylogenetic trees. In principle, there's a parallel option (via MPI), but it doesn't work for me. Instead of spending lots of time debugging MPI, I was wondering if I could run the bootstraps independently (with a user-supplied tree) on my cluster. So the real question is: Given a list of trees in Newick format, how do I calculate the bootstrap support for the original tree?

Answer: Bootstrap replicates are independent, so to the first part of your question: yes, you can simply create however many bootstrapped matrices you require (e.g. 1000) and run a tree search on each of those on separate nodes of your cluster. You can use seqboot for this, or other utilities with similar functionality. As a plug for [Bio::Phylo](http://search.cpan.org/dist/Bio-Phylo/) you might, for example, do the following:

use Bio::Phylo::IO 'parse';

my ($matrix) = @{ parse(
-format => 'nexus', # or any of the other supported formats, e.g. 'phylip'
-file => 'myfile.nex', # or a string, url or handle
-as_project => 1,
)->get_matrices };

for ( 1 .. 1000 ) {
my $bootstrapped = $matrix->bootstrap;
open my $outfh, '>', "myfile{$_}.nex" or die $!;
print $outfh "#NEXUS\n", $bootstrapped->to_nexus;
}

...which gives you a thousand bootstrapped versions of 'myfile.nex', with names 'myfile[1..1000].nex'. Then, you run a tree search on each of those, concatenate the resulting trees into a list of newick strings (as per your original query) and do the following:

use Bio::Phylo::IO 'parse';

my $forest = parse(
-format => 'newick',
-file => 'mytrees.dnd',
);

my $consensus = $forest->make_consensus( -branches => 'frequency' );
open my $outfh, '>', 'consensus.dnd' or die $!;
print $outfh $consensus->to_newick;

This gives you a little more flexibility in terms of the file formats you can use beyond phylip and newick (and, consequently, the tree searching programs you can use) but other than that it is equivalent to using seqboot and consense - and, depending on the number of sequences and bootstrap replicates, there might be performance issues with using perl.