[Michael Love]

note that the replicates are right on top of each other in the PCA plot. Are these technical or biological replicates?

The dispersion is calculated based on variance within conditions, so the dispersion is not necessarily large though you have large differences across conditions.

I'm not so familiar with microbial analysis. I'd guess, like others mentioned above, that you have many genes with counts for only one species. And there is not a clear group of genes which are not DE across the conditions. This makes normalization difficult, as the automatic methods within DESeq or edgeR are based on the assumption that there are enough genes that are not DE, such that robust measures like median or trimmed mean can find the center of the distribution of log ratios of samples.

Is there a set of genes that the biologists suspect might be equally expressed across the groups?

[alyamahmoud]

Hi Michael

These are biological replicates.

The dispersion is calculated based on variance within conditions, so the dispersion is not necessarily large though you have large differences across conditions.

I am not sure I get what you mean here.

There is a set of genes that the biologist know should be varying and these are non-metagenomics samples; single species per sample.

What would you suggest ?

[alyamahmoud]

If I use collapseReplicates the number of DEG decreases massively, however, this doesn't improve the ma plot (attached)!!

any help ?

[gringer]

genome

Just a general question related to this (based on the mapping IDs you have provided), is there anyone here on seqanswers who has successfully done a DESeq / DESeq2 run on E. coli?

[alyamahmoud]

hierarchial model

Is there any objection to applying a hierarchial model on the normalized counts ? I tried limma analysis on the normalized counts, the MA plot is also attached.

[Michael Love]

collapseReplicates() is for technical replicates only. We obviously do not recommend collapsing biological replicates, as you throw away information from the experiment.

The "problem" of too many p-values or a p-value distribution with a spike at 0 means that you have many large differences across the conditions.

Can you say more about the genes here? If you are sequencing multiple species, what is the relation of each species to the reference genome/transcriptome to which the reads were aligned?

[alyamahmoud]

They are not multiple species, only one species (same as reference) but under different environmental conditions (different pH ranges, anaerobic, water vs wt that is aerobic)

[alyamahmoud]

is it wrong to apply a hierarchal model ? it reduces the number of sig genes drastically but the the MA plot looks better I think ?

[Michael Love]

The last plot you posted, hm_norm_counts, was not an MA plot. it was logFC ~ adjusted p value.

An MA plot is logFC ~ log of mean counts, or mean of log counts.

Also, note the scale of the y axis is much larger than the previous plots.

Regarding the goal of reducing the number of sig genes: I don't know if you've reduced the number of significant genes for better or for worse. We can easily reduce the number of genes, either by reducing the FDR threshold or increasing the lfcThreshold argument.

We know for sure, from the PCA plot, that the differences between the samples are very large compared to the variation between biological replicates.

Could you send the dds object to me privately, so I can have a look?

My email is listed here:

maintainer("DESeq2")

[gringer]

That plot also doesn't look wonderful, presumably because you've got p-val on the X axis. MA plot is usually log fold change on Y, and average log expression on the X.

Are you able to do a scatter plot of the raw counts for each experiment, preferably log-transformed or using the VST from DESeq/DESeq2? If you're not getting a line that distributes around y=x with those plots, it's probably not a good idea trying to shoehorn in a differential expression analysis.

[Michael Love]

Yes, David's right.

Here's a pairs plot of your counts in the log scale

y <- log10(counts(dds)+1)
pairs(y, panel = function(...) smoothScatter(..., nrpoints = 0, add = TRUE),lower.panel=NULL)

For the samples other than 'water', we can see the diagonal line that would specify a log fold change of 0 between the two samples. This is the line that DESeq and edgeR use for defining a scaling factor for normalizing for sequencing depth.

However, for water vs others, a simple scaling factor automatically detected from the data will not work.

For the scatterplot of 1 vs 3 and 1 vs 10, there seems to be a faint line of genes on the diagaonal. Maybe you can investigate what is special about these genes. It is possible that nearly all the genes are differentially expressed (upregulated in the treated groups), but then the experiment really should use spike in controls for normalization.

I wonder if the experimental protocol might have been different for the water samples?

Another option for analysis would be to remove the water samples and use the 'contrast' argument to just compare the treatment groups against each other.