DESeq2 plotMA: data appears incorrect, patterned

[BPedler]

Hi,
I'm trying to learn DESeq2 using a simplified data set with 3 controls ("ctr"), and 3 treatments ("koh"). I used the summary from an earlier thread (MDonlin;120724) to get started. I am not getting any error messages, but the output from plotMA does not appear as it does in the DESeq2 "airway" vignette. It looks patterned in a non-random way suggesting something is incorrect (plot attached).

If anyone has any suggestions, I'd much appreciate the help.

Best,
Byron

code:
> library("DESeq2")
>
> #generate count table from text file
> GeneCountTable <- read.table("KM272.d3.ctr.koh.txt", header=TRUE, row.names=1)
>
> head(GeneCountTable)
ctr1d3 ctr2d3 ctr3d3 koh1d3 koh2d3 koh3d3
gi|10000000001|loc|edl|EDL_NS211000002.1| 0 3 0 0 0 0
gi|10000000002|loc|edl|EDL_NS211000003.1| 0 0 0 0 0 0
gi|10000000003|loc|edl|EDL_NS211000004.1| 0 0 0 0 0 0
gi|10000000008|loc|edl|EDL_NS211000009.1| 0 0 0 0 0 0
gi|10000000011|loc|edl|EDL_NS211000012.1| 0 0 0 0 0 0
gi|10000000018|loc|edl|EDL_NS211000019.1| 0 0 0 0 0 0
>
> #define samples
> samples <- data.frame(row.names=c("ctr1d3","ctr2d3","ctr3d3","koh1d3","koh2d3","koh3d3"), condition=as.factor(c(rep("ctr",3),rep("koh",3))))
> samples
condition
ctr1d3 ctr
ctr2d3 ctr
ctr3d3 ctr
koh1d3 koh
koh2d3 koh
koh3d3 koh
>
> #generate DESeq dataset
> KM272dds <- DESeqDataSetFromMatrix(countData = GeneCountTable, colData=samples, design=~condition)
> KM272dds
class: DESeqDataSet
dim: 559312 6
exptData(0):
assays(1): counts
rownames(559312): gi|10000000001|loc|edl|EDL_NS211000002.1|
gi|10000000002|loc|edl|EDL_NS211000003.1| ... gi|9964628|ref|NP_064758.1|
gi|99878752|ref|YP_615055.1|
rowRanges metadata column names(0):
colnames(6): ctr1d3 ctr2d3 ... koh2d3 koh3d3
colData names(1): condition
>
> #run DESeq on dataset
> KM272dds_1 <- DESeq(KM272dds)
estimating size factors
estimating dispersions
gene-wise dispersion estimates
mean-dispersion relationship
final dispersion estimates
fitting model and testing
>
> #generate results table
> KM272_res <- results(KM272dds_1)
>
> #reorder results table by lowest adjusted P-value
> KM272_resOrdered <- KM272_res[order(KM272_res$padj),]
> head(KM272_resOrdered)
log2 fold change (MAP): condition koh vs ctr
Wald test p-value: condition koh vs ctr
DataFrame with 6 rows and 6 columns
baseMean log2FoldChange lfcSE stat pvalue padj
<numeric> <numeric> <numeric> <numeric> <numeric> <numeric>
gi|426411412|ref|YP_007031511.1| 1257.5658 8.121059 0.8434747 9.628101 6.084362e-22 1.696016e-17
gi|537453526|ref|YP_008487251.1| 967.6043 7.914203 0.9057204 8.738021 2.372297e-18 3.306389e-14
gi|152995336|ref|YP_001340171.1| 294.0563 11.018376 1.3738295 8.020192 1.055802e-15 7.357622e-12
gi|224584543|ref|YP_002638341.1| 236.0105 9.143821 1.1382139 8.033482 9.474453e-16 7.357622e-12
gi|333907837|ref|YP_004481423.1| 291.8945 11.013159 1.3848704 7.952484 1.828092e-15 1.019161e-11
gi|285019583|ref|YP_003377294.1| 147.3368 6.864467 0.8663821 7.923140 2.315875e-15 1.075917e-11
>
> #write CSV file
> write.csv(KM272_resOrdered,file="KM272_RNA_results.csv")
>
> #summarize results
> summary(KM272_res)

out of 139894 with nonzero total read count
adjusted p-value < 0.1
LFC > 0 (up) : 1927, 1.4%
LFC < 0 (down) : 17, 0.012%
outliers [1] : 272, 0.19%
low counts [2] : 111747, 80%
(mean count < 0.5)
[1] see 'cooksCutoff' argument of ?results
[2] see 'independentFiltering' argument of ?results

> #VISUALIZE
> #In DESeq2, the function plotMA shows the log2 fold changes attributable to a given variable over the mean of normalized counts. Points will be colored red if the adjusted p value is less than 0.1. Points which fall out of the window are plotted as open triangles pointing either up or down.
> plotMA(KM272_res, main="DESeq2", ylim=c(-15,15))


> sessionInfo()
R version 3.2.2 (2015-08-14)
Platform: x86_64-apple-darwin13.4.0 (64-bit)
Running under: OS X 10.10.5 (Yosemite)

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] parallel stats4 stats graphics grDevices utils datasets methods base

other attached packages:
[1] DESeq2_1.8.1 RcppArmadillo_0.5.600.2.0 Rcpp_0.12.1
[4] GenomicRanges_1.20.8 GenomeInfoDb_1.4.3 IRanges_2.2.7
[7] S4Vectors_0.6.6 BiocGenerics_0.14.0

loaded via a namespace (and not attached):
[1] RColorBrewer_1.1-2 futile.logger_1.4.1 plyr_1.8.3 XVector_0.8.0
[5] futile.options_1.0.0 tools_3.2.2 rpart_4.1-10 digest_0.6.8
[9] RSQLite_1.0.0 annotate_1.46.1 gtable_0.1.2 lattice_0.20-33
[13] DBI_0.3.1 proto_0.3-10 gridExtra_2.0.0 genefilter_1.50.0
[17] stringr_1.0.0 cluster_2.0.3 locfit_1.5-9.1 nnet_7.3-11
[21] grid_3.2.2 Biobase_2.28.0 AnnotationDbi_1.30.1 XML_3.98-1.3
[25] survival_2.38-3 BiocParallel_1.2.21 foreign_0.8-66 latticeExtra_0.6-26
[29] Formula_1.2-1 geneplotter_1.46.0 ggplot2_1.0.1 reshape2_1.4.1
[33] lambda.r_1.1.7 magrittr_1.5 scales_0.3.0 Hmisc_3.17-0
[37] MASS_7.3-44 splines_3.2.2 xtable_1.7-4 colorspace_1.2-6
[41] stringi_0.5-5 acepack_1.3-3.3 munsell_0.4.2
>

[Michael Love]

One note: this dataset has very low counts overall, between 1-10 reads for most genes. Compare with the MA plot for the airway dataset: http://www.bioconductor.org/help/wor...seqGene/#plots

In your plot, those lines are 1 read vs 0, 2 reads vs 0, 3 reads vs 0, etc.

You can maybe see that more clearly if you set DESeq(dds, betaPrior=FALSE) and then plotMA(), although the LFC shrinkage (i.e. betaPrior=TRUE) seems useful here so you can compare 10 vs 0 and 100 vs 0, which un-shrunken would both go to infinity.