I am using DESeq2 to run a multifactor paired analysis to determine the change in gene expression between men and women due to a treatment on 70 males and 45 females. Because I have a differing number of samples for each condition, I'm using the following steps to create a custom design matrix: C: DESeq2 paired and multi factor comparison
Everything is fine until estimateDispersionsGeneEst(dds, modelMatrix=mm1) gets called. The command runs for a while then exits and terminates R. Before termination, the following gets called:
using supplied model matrix
*** glibc detected *** /opt/R-3.3.1/lib64/R/bin/exec/R: double free or corruption (!prev): 0x0000000042e3d4b0 ***
======= Backtrace: =========
(a bunch of calls through the R library)
I looked up what "glibc detected" indicates, and from what I could understand, it refers to erroneous memory allocation.
My code is:
samples <- read.csv(sampleTable.csv)
sampleTable <- data.frame(sample = samples$PID, fileName = samples$featureCounts_PID, sex = samples$Sex, condition = samples$Condition, nested = samples$Nested)
#nested is a factor
sampleTable$sex <- relevel(sampleTable$sex, ref = "Male")
sampleTable$condition <- relevel(sampleTable$condition, ref = "pre")
#even though counts are from featureCounts, formatted results to be same as HTSeq-count results to use following import method
dds <- DESeqDataSetFromHTSeqCount(sampleTable = sampleTable, directory = directory, design = ~ 1)
model1 = model.matrix(~ sex + sex:nested + sex:condition, colData(dds))
idx <- which(colSums(model1 == 0) == nrow(model1))
model1 <- model1[,-idx] # removing the columns with all 0's
model0 <- model.matrix(~ sex + sex:nested, colData(dds)) # this matrix is necessary to run nbinomLRT
design(dds) <- ~ sex + sex:nested + sex:condition
dds <- estimateSizeFactors(dds)
#stops working once I call next method
dds <- estimateDispersionsGeneEst(dds, modelMatrix=model1)
#glibc error, ever get to following methods
dds <- estimateDispersionsFit(dds)
dds <- estimateDispersionsMAP(dds, modelMatrix=model1)
dds <- nbinomLRT(dds, full=model1, reduced=model0)
Any help would be greatly appreciated!
> sessionInfo()
R version 3.3.1 (2016-06-21)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: CentOS release 6.8 (Final)
locale:
[1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
[3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
[5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
[7] LC_PAPER=en_US.UTF-8 LC_NAME=C
[9] LC_ADDRESS=C LC_TELEPHONE=C
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
attached base packages:
[1] parallel stats4 stats graphics grDevices utils datasets
[8] methods base
other attached packages:
[1] IHW_1.2.0 DESeq2_1.14.1
[3] SummarizedExperiment_1.4.0 Biobase_2.34.0
[5] GenomicRanges_1.26.2 GenomeInfoDb_1.10.2
[7] IRanges_2.8.1 S4Vectors_0.12.1
[9] BiocGenerics_0.20.0
loaded via a namespace (and not attached):
[1] genefilter_1.56.0 locfit_1.5-9.1 slam_0.1-40
[4] splines_3.3.1 lattice_0.20-33 colorspace_1.3-2
[7] htmltools_0.3.5 base64enc_0.1-3 survival_2.40-1
[10] XML_3.98-1.5 foreign_0.8-66 DBI_0.5-1
[13] BiocParallel_1.8.1 RColorBrewer_1.1-2 plyr_1.8.4
[16] stringr_1.1.0 zlibbioc_1.20.0 munsell_0.4.3
[19] gtable_0.2.0 htmlwidgets_0.8 memoise_1.0.0
[22] latticeExtra_0.6-28 knitr_1.15.1 geneplotter_1.52.0
[25] fdrtool_1.2.15 AnnotationDbi_1.36.2 htmlTable_1.9
[28] Rcpp_0.12.9 acepack_1.4.1 xtable_1.8-2
[31] backports_1.0.5 scales_0.4.1 checkmate_1.8.2
[34] Hmisc_4.0-2 annotate_1.52.1 XVector_0.14.0
[37] lpsymphony_1.2.0 gridExtra_2.2.1 ggplot2_2.2.1
[40] digest_0.6.12 stringi_1.1.2 grid_3.3.1
[43] tools_3.3.1 bitops_1.0-6 magrittr_1.5
[46] lazyeval_0.2.0 RCurl_1.95-4.8 tibble_1.2
[49] RSQLite_1.1-2 Formula_1.2-1 cluster_2.0.4
[52] Matrix_1.2-6 data.table_1.10.4 assertthat_0.1
[55] rpart_4.1-10 nnet_7.3-12
my colData looks like the following, where condition, sex, and nested are all appropriately leveled factors