Question: DESeq2 glibc error
0
gravatar for gregory.l.stone
2.8 years ago by
gregory.l.stone10 wrote:

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!

 

ADD COMMENTlink modified 2.8 years ago by Michael Love26k • written 2.8 years ago by gregory.l.stone10

> 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

ADD REPLYlink written 2.8 years ago by gregory.l.stone10

my colData looks like the following, where condition, sex, and nested are all appropriately leveled factors

sample condition sex nested
sample1 pre M 1
sample1 post M 1
sample2 pre M 2
sample2 post M 2
sample70 pre F 1
sample 70 post F 1
... ... ... ...

 

ADD REPLYlink written 2.8 years ago by gregory.l.stone10
Answer: DESeq2 glibc error
0
gravatar for Michael Love
2.8 years ago by
Michael Love26k
United States
Michael Love26k wrote:
Can you check that you aren't feeding full or reduced a matrix with a column of all zeros.
ADD COMMENTlink written 2.8 years ago by Michael Love26k

It seems that I have reached a positing limit (accurately reflecting my difficulties), so I am answering from another account. I am not feeding a full or reduced matrix with a column of all zeros. I remove these columns by:

idx <- which(colSums(model1 == 0) == nrow(model1))
model1 <- model1[,-idx]

and I double checked by doing idx_check <- which(colSums(model1 == 0) == nrow(model1)) and confirmed that idx_check is empty

ADD REPLYlink written 2.8 years ago by gstone20
And reduced=model0? You checked this as well?
ADD REPLYlink written 2.8 years ago by Michael Love26k

Correct. I did the same to that matrix, although I do not get to the command that takes that matrix as a parameter since R exits at the estimateDispersionsGeneEst() call

ADD REPLYlink written 2.8 years ago by gstone20

Oh gotcha. I see that now.

The only other thing left that seems strange in your code is why you've taken featureCounts output and used HTSeq constructor function. There could be something funny happening there. Why not use our recommended constructor function for a featureCount count matrix: DESeqDataSetFromMatrix?

ADD REPLYlink written 2.8 years ago by Michael Love26k

I used the HTSeq constructor function because it seemed much easier to use than that for a count matrix, and it was easy to format featureCount output appropriately. It is a pretty janky approach though, so I will go back and use the appropriate constructor function. Thank you for your help and for being so generous with your time!

ADD REPLYlink written 2.8 years ago by gstone20
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 16.09
Traffic: 223 users visited in the last hour