non-typical p values distribution running DESeq

[alyamahmoud]

Hi All

I have 2 reps*5 conditions (4 + control). I ran followed the nature protocol to come to differentially expressed genes between each of the 4 conditions relative to the control.

The distribution of p values looks is attached for one of the conditions as well as the corresponding MA plot. I am getting too many differentially expressed genes (below) and the p values distribution doesn't look like expected.

Would you please advise what I could be missing ? or is such a distribution of p values expected in some cases and why ?

Thank you very much
Alyaa

> table (res$padj < 0.1)
FALSE  TRUE 
 3245  3578 

> table (res$padj < 0.01)
FALSE  TRUE 
 4813  2010 

> table (res$padj < 0.05)
FALSE  TRUE 
 3862  2961

[gringer]

Based on the MA plot, I'm guessing you're using DESeq2, rather than DESeq, which is good.

However, your MA plot looks crazy. It should be distributed around the y axis (0 log2 fold change).

I've got no idea what would do that, but it certainly indicates something screwy is going on. Are you sub-sampling genes prior to running them through DESeq2? Are you using normalised counts as input, instead of raw counts? What command did you run to produce this MA plot?

[alyamahmoud]

Hi gringer

Thanks for the prompt reply.

I am not subsampling, and I am using the raw counts as input to DESeq and not DESeq2.

I ran the following to generate the MA plot:

cds =  newCountDataSet(counts_table, conds)
cds <- estimateSizeFactors(cds)
cds <- estimateDispersions(cds)
res = nbinomTest(cds, "water", "aerobic") # one of the conditions vs control
plotMA(res)

[gringer]

I am not subsampling, and I am using the raw counts as input to DESeq and not DESeq2.

Oh, okay. Do you get the same results when using DESeq2?

The only reason I can think of off the top of my head why increased expression in both groups would result in increased log2 fold change is if one of the groups had 0 expression for all genes. Can you show the first few lines of your results, i.e. "head(res)"? I think DESeq (v1, not v2) should report estimated/normalised expression for each group in that result.

If that is the problem, then you may have chosen your condition names incorrectly in the nbinomTest command:

> conditions(cds)
[1] "water" "aerobic" # should be something like this

Or alternatively, the count table or condition list could be formatted incorrectly:

> colnames(counts_table)
[1] "water_r1" "water_r2" "aerobic_r1" "aerobic_r2" # something like this
> dim(counts_table)
[1] 30215 4 # should be something like this
> conds
[1] "water" "water" "aerobic" "aerobic"  # should be something like this

But please try DESeq2. It will complain a bit louder when you do things wrong, which will hopefully give you more insights into what went wrong.

[alyamahmoud]

I tried DESeq2 and the results are not very different. The p values distribution and the MA plot using DESeq2 are attached.


Here are the formats for the counts table and design:

>conditions(cds)
    water_1     water_2       ph5_1       ph5_2       ph9_1       ph9_2 
      water       water         ph5         ph5         ph9         ph9 
anaerobic_1 anaerobic_2   aerobic_1   aerobic_2 
  anaerobic   anaerobic     aerobic     aerobic 
Levels: aerobic anaerobic ph5 ph9 water

> colnames(counts_table)
 [1] "water_1"     "water_2"     "ph5_1"       "ph5_2"       "ph9_1"      
 [6] "ph9_2"       "anaerobic_1" "anaerobic_2" "aerobic_1"   "aerobic_2"

>conds
 [1] water     water     ph5       ph5       ph9       ph9       anaerobic
 [8] anaerobic aerobic   aerobic  
Levels: aerobic anaerobic ph5 ph9 water

[alyamahmoud]

Am I formatting the conditions file wrongly:

This is my exp_design file:

>exp_design
            condition
water_1         water
water_2         water
ph5_1             ph5
ph5_2             ph5
ph9_1             ph9
ph9_2             ph9
anaerobic_1 anaerobic
anaerobic_2 anaerobic
aerobic_1     aerobic
aerobic_2     aerobic

and these are the commands in DESeq and DEseq2, respectively.

conds = exp_design$condition
cds =  newCountDataSet(counts_table, conds)

ddsFullCountTable <- DESeqDataSetFromMatrix(countData = counts_table, colData = exp_design, design = ~condition)

[gringer]

Can you show the first few lines of results (and/or your count table)? The rest of what you've got looks fine.

[alyamahmoud]

> head (counts_table)
              water_1 water_2 ph5_1 ph5_2 ph9_1 ph9_2 anaerobic_1 anaerobic_2
FN649414.4579     500     243   133   647   141   114         222         106
FN649414.7957      23      20    10    91    12    13          13           6
FN649414.7767     135      50    55   321    52    43          96          53
p948.168            5       0     0     0     0     0          66          28

              aerobic_1 aerobic_2
FN649414.4579        50       113
FN649414.7957        16        34
FN649414.7767        33       101
p948.168              0         0

[alyamahmoud]

I tried

>cdsFilt = estimateDispersions(cdsFilt, method  = "blind", sharingMode="fit-only)

and the MA plot is attached, does it look any better ?

[gringer]

No. You're expecting points looking like a triangle (or diamond) shaped wedge elongated along the Y axis, centered on the Y axis. I've attached an example based on DESeq results. The DESeq2 plot should look similar, but narrows down to a point for low expression values. If you're not getting that (and all other steps check out), then it suggests that your experimental conditions aren't appropriately controlled.

[gringer]

You wouldn't happen to be trying to do a metagenomic analysis, would you? If you have a different sample population for each condition, then that might explain the differences that you're seeing.

[alyamahmoud]

I will check with the biologist who gave me the data, but I don't think they are coming from different population.

How reliable/ir-reliable would be the results of the DEG accordingly ?

[gringer]

Do you have mapping percentages for your genome? If they're wildly different, that might point to sample differences that could greatly influence the results.

What your MA plots seem to be showing is that one of the populations is producing normalised counts for genes that are effectively zero.

You can look at the count input data (or DESEq-normalised data) to see if this is the case -- just looking at total counts per experiment should show odd patterns. If the differences are as much as they look from the MA plot, you won't need DESeq at all to find differences.

[Michael Love]

Indeed, something looks strange about this experiment. Try running the PCA steps in the vignette to see how the samples are distributed. It looks like the different conditions are very different from each other for many rows. This can make the normalization difficult without spike-in controls.

Can you describe the experiment in more details? Is this standard RNA-Seq data?

[alyamahmoud]

Hi Michael

Thanks for your reply.

I attached the PCA plot.

Yes, these are very different conditions for some bacterial species, paired-end standard RNAseq, % of mapped reads >= 97%.

How can I handle this data without spike-in controls ? How much effect does this massive dispersion have on the final DEG ?