Hi all,
I am trying to identify genes that are specific to a cell type of interest by comparison with a large number of other cell types. To this end, I performed differential gene expression using DESeq2 using the following:
dds = DESeqDataSetFromMatrix(countData = as.matrix(dtCountsDDS),
colData = DataFrame(dtSampleDDS),
design = ~group)
where "group" goes from 0-48 with 0 is the group with the cell type of interest and 1-48 gives 48 other cell types which I am comparing to. Thus, the coldata looks like this:
description group
cell_type_of_interest_rep1 0
cell_type_of_interest_rep2 0
cell_type_of_interest_rep3 0
other_cell_type01_rep1 1
other_cell_type01_rep2 1
other_cell_type01_rep3 1
other_cell_type02_rep1 2
other_cell_type02_rep2 1
...
other_cell_type48_rep1 48
other_cell_type48_rep2 48
other_cell_type48_rep3 48
I then performed DESeq2 as follows:
dds = DESeq(dds, quiet = FALSE, betaPrior = FALSE, parallel = TRUE)
After which, I retrived the pairwise Wald comparison results for group0-vs-group1, group0-vs-group2, ..., group0-vs-group48 to find genes that are differentially expressed in the cell type of interest by the following:
res = results(dds, alpha = 0.05, independentFiltering = TRUE, contrast = c("group", 0, i))
where "i" goes from 1 to 48.
From the plotCounts, I am expecting the gene NANOG to be differentially expressed against most cell types, which is evident from the plotCounts.
However, the DESeq2 results only gives an FDR of 0.5-0.6 in all cases.
groupID geneID baseMean LFC lfcSE stat pvalue FDR
1 ENSG00000111704 9.399662 10.6784175 14.00173 0.76265003 0.4456721 0.5687196
2 ENSG00000111704 9.399662 8.7716355 14.00027 0.62653347 0.5309651 0.6534532
3 ENSG00000111704 9.399662 9.4506731 14.00173 0.67496482 0.4996981 0.633392
4 ENSG00000111704 9.399662 9.3446924 14.00173 0.66739571 0.5045194 0.6249148
5 ENSG00000111704 9.399662 10.3797534 14.00173 0.74131952 0.4584997 0.592654
6 ENSG00000111704 9.399662 10.6885146 14.00173 0.76337116 0.4452421 0.5774063
7 ENSG00000111704 9.399662 11.5510343 14.00173 0.82497212 0.4093874 0.565289
8 ENSG00000111704 9.399662 10.5822928 14.00173 0.75578484 0.4497782 0.5608032
9 ENSG00000111704 9.399662 9.5290049 14.00173 0.68055926 0.4961504 0.6209838
10 ENSG00000111704 9.399662 9.0412122 15.67196 0.57690359 0.5640046 0.6806634
11 ENSG00000111704 9.399662 10.8099388 14.00173 0.77204326 0.4400888 0.5565791
12 ENSG00000111704 9.399662 10.5248998 14.00173 0.75168584 0.45224 0.5566365
13 ENSG00000111704 9.399662 10.1456759 14.00173 0.72460176 0.4686964 0.5728484
14 ENSG00000111704 9.399662 8.7660278 14.00173 0.62606762 0.5312706 0.6353354
15 ENSG00000111704 9.399662 10.7080225 14.00173 0.76476442 0.4444118 0.5620308
16 ENSG00000111704 9.399662 9.3813843 14.00173 0.67001624 0.5028474 0.6240379
17 ENSG00000111704 9.399662 10.4812734 15.67196 0.66879132 0.5036286 0.6294616
18 ENSG00000111704 9.399662 10.9714755 14.00173 0.78358017 0.4332865 0.5349624
19 ENSG00000111704 9.399662 9.6212037 13.99186 0.68762874 0.4916866 0.6054409
20 ENSG00000111704 9.399662 8.0895122 13.97886 0.57869604 0.5627943 0.670864
21 ENSG00000111704 9.399662 9.6592287 14.00173 0.68985981 0.4902824 0.5920907
22 ENSG00000111704 9.399662 10.9413107 14.00173 0.7814258 0.4345521 0.5801206
23 ENSG00000111704 9.399662 9.354374 13.99711 0.66830755 0.5039373 0.6267258
24 ENSG00000111704 9.399662 7.7500671 13.98033 0.55435526 0.5793357 0.6769946
25 ENSG00000111704 9.399662 9.9603945 14.00173 0.711369 0.4768556 0.6080437
26 ENSG00000111704 9.399662 10.2005031 14.00173 0.7285175 0.4662969 0.6046647
27 ENSG00000111704 9.399662 9.1613179 15.67196 0.58456731 0.5588387 0.6709067
28 ENSG00000111704 9.399662 10.3497801 14.00173 0.73917883 0.4597984 0.5887615
29 ENSG00000111704 9.399662 8.929098 14.00173 0.63771405 0.5236598 0.6456075
30 ENSG00000111704 9.399662 11.1910169 14.00173 0.79925976 0.4241398 0.5865799
31 ENSG00000111704 9.399662 9.8669264 13.99732 0.70491541 0.4808629 0.6312367
32 ENSG00000111704 9.399662 10.4773122 14.00173 0.74828715 0.454287 0.7951142
33 ENSG00000111704 9.399662 -0.1575527 13.95228 -0.01129226 0.9909903 0.9973648
34 ENSG00000111704 9.399662 8.293104 14.00173 0.59229151 0.5536554 0.6722364
35 ENSG00000111704 9.399662 8.7458279 13.99514 0.62491901 0.5320242 0.6486637
36 ENSG00000111704 9.399662 7.8522695 14.00099 0.56083671 0.5749089 0.6777484
37 ENSG00000111704 9.399662 10.8459836 15.67196 0.69206282 0.4888979 0.6175993
38 ENSG00000111704 9.399662 10.390056 14.00173 0.74205533 0.4580538 0.5674383
39 ENSG00000111704 9.399662 8.6865408 14.00173 0.62039068 0.5350006 0.6595885
40 ENSG00000111704 9.399662 10.4439862 14.00173 0.74590701 0.4557236 0.5772359
41 ENSG00000111704 9.399662 9.8352453 14.00173 0.70243088 0.4824105 0.6273037
42 ENSG00000111704 9.399662 9.3964659 14.00173 0.67109335 0.5021611 0.629528
43 ENSG00000111704 9.399662 10.3623496 14.00173 0.74007654 0.4592536 0.5941857
44 ENSG00000111704 9.399662 8.3683002 15.67196 0.53396627 0.5933649 0.7084019
45 ENSG00000111704 9.399662 9.8391147 15.67196 0.62781632 0.5301243 0.6496046
46 ENSG00000111704 9.399662 10.4257006 14.00173 0.74460106 0.4565129 0.7593986
47 ENSG00000111704 9.399662 6.2449039 13.95703 0.44743778 0.654559 0.7649179
48 ENSG00000111704 9.399662 10.3231834 14.00173 0.7372793 0.4609525 0.5941485
Throughout this entire DESeq2 run, there is only a warning on the mean-dispersion fitting (reproduced below) and I think that is causing the lfcSE to be very large and consequently giving an insignificant FDR.
-- note: fitType='parametric', but the dispersion trend was not well captured by the
function: y = a/x + b, and a local regression fit was automatically substituted.
specify fitType='local' or 'mean' to avoid this message next time.
Thus, I have two questions:
-
Is there anyway to remedy this apparently strange result?
-
Is it sensible to do all the comparisons in one DESeq2 pass? Or should I perform 48 separate individual DESeq2 runs to compare the cell type of interest with each of the other cell types?
Hi Michael,
Thanks for the prompt response. I ran my dataset using your code and NANOG is now DE with a FDR of 10^-8 to 10^-16.
However, I have two further questions:
The calculation took a few hours as it is ran in serial. Is it possible to run this code in parallel instead?
I took a close look at the results. For group 1-10, the counts of NANOG is 0. So I was expecting the comparison group0-vs-group1, group0-vs-group2, ..., group0-vs-group10 to yield similar results. However, the FDR still varies quite a bit. Is there any reason why?
The differences in p-values are due to differences in the size factor. A "0" is not exactly a "0" due to the offset from sequencing depth. You shouldn't really interpret a difference between FDR of 1e-10 or 1e-12, this isn't giving you extra information. The question is really just, for an 'alpha' of interest to you, is padj < alpha or not.
In the next version of DESeq2 (1.20), you can just specify minmu directly to DESeq().
For now you can use a simple wrapper that uses BiocParallel directly:
dds <- makeExampleDESeqDataSet() library(BiocParallel) nworkers <- 3 register(MulticoreParam(nworkers)) # this makes a partition of the genes to send to diff workers idx <- factor(sort(rep(seq_len(nworkers),length.out=nrow(dds)))) dds <- estimateSizeFactors(dds) minmu <- 1e-3 # split and then bind together along partition dds <- do.call(rbind, bplapply(levels(idx), function(l) { estimateDispersionsGeneEst(dds[idx == l,], minmu=minmu) })) # this cannot be run in parallel, needs to see all genes dds <- estimateDispersionsFit(dds) # these two steps can be combined into one split and bind operation dds <- do.call(rbind, bplapply(levels(idx), function(l) { d <- estimateDispersionsMAP(dds[idx == l,]) nbinomWaldTest(d) })) # this is fast, does not benefit from parallel, # if you only compare to reference level res <- results(dds)Thanks for the quick wrapper solution. And I am looking forward to the next version of DESeq2. :)