Question

Error message from updated DESeq2

0

Entering edit mode

asaferali • 0

@asaferali-6904

Last seen 10.0 years ago

Canada

Hello,

I recently updated to the latest version of DESeq2. When I try to run the exact same script that I had used in the previous version of DESeq2 the dds<-DESeq(dds) command gives me the following error:

Error in averagePriorsOverLevels(objectNZ, betaPriorVar) :
beta prior for Subject.T.1.,Condition.T.Control. is not greater than 0

This is my sessionInfo:

R version 3.1.3 (2015-03-09)
Platform: x86_64-apple-darwin13.4.0 (64-bit)
Running under: OS X 10.10.5 (Yosemite)

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

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

other attached packages:
[1] BiocInstaller_1.16.5 edgeR_3.8.6
[3] limma_3.22.7 DESeq2_1.6.3
[5] RcppArmadillo_0.6.200.2.0 Rcpp_0.12.1
[7] GenomicRanges_1.18.4 GenomeInfoDb_1.2.5
[9] IRanges_2.0.1 S4Vectors_0.4.0
[11] BiocGenerics_0.12.1 Rcmdr_2.1-7
[13] RcmdrMisc_1.0-3 sandwich_2.3-4
[15] car_2.0-25 FactoMineR_1.31.4

loaded via a namespace (and not attached):
[1] abind_1.4-3 acepack_1.3-3.3 annotate_1.44.0
[4] AnnotationDbi_1.28.2 base64enc_0.1-3 BatchJobs_1.6
[7] BBmisc_1.9 Biobase_2.26.0 BiocParallel_1.0.3
[10] brew_1.0-6 checkmate_1.6.3 class_7.3-14
[13] cluster_2.0.3 codetools_0.2-14 colorspace_1.2-6
[16] DBI_0.3.1 digest_0.6.8 e1071_1.6-7
[19] fail_1.3 flashClust_1.01-2 foreach_1.4.3
[22] foreign_0.8-66 Formula_1.2-1 genefilter_1.48.1
[25] geneplotter_1.44.0 ggplot2_1.0.1 grid_3.1.3
[28] gtable_0.1.2 Hmisc_3.17-0 iterators_1.0.8
[31] lattice_0.20-33 latticeExtra_0.6-26 leaps_2.9
[34] lme4_1.1-10 locfit_1.5-9.1 magrittr_1.5
[37] MASS_7.3-44 Matrix_1.2-2 MatrixModels_0.4-1
[40] mgcv_1.8-9 minqa_1.2.4 munsell_0.4.2
[43] nlme_3.1-122 nloptr_1.0.4 nnet_7.3-11
[46] pbkrtest_0.4-2 plyr_1.8.3 proto_0.3-10
[49] quantreg_5.19 RColorBrewer_1.1-2 readxl_0.1.0
[52] reshape2_1.4.1 rpart_4.1-10 RSQLite_1.0.0
[55] scales_0.3.0 scatterplot3d_0.3-36 sendmailR_1.2-1
[58] SparseM_1.7 stringi_1.0-1 stringr_1.0.0
[61] survival_2.38-3 tcltk_3.1.3 tcltk2_1.2-11
[64] tools_3.1.3 XML_3.98-1.3 xtable_1.8-0
[67] XVector_0.6.0 zoo_1.7-12

deseq2 beta prior error • 3.0k views

ADD COMMENT • link 10.0 years ago asaferali • 0

0

Entering edit mode

I'm also encountering the error:

> system.time(

+   dds <- DESeq(dds)
+ )
estimating size factors
estimating dispersions
gene-wise dispersion estimates
mean-dispersion relationship
final dispersion estimates
fitting model and testing
Error in averagePriorsOverLevels(objectNZ, betaPriorVar) : 
  beta prior for as.factor.HTN.0_No is not greater than 0
In addition: Warning messages:
1: In design(object) == formula(~1) :
  longer object length is not a multiple of shorter object length
2: In modelFormula == formula(~1) :
  longer object length is not a multiple of shorter object length
Timing stopped at: 11698 43893 5418

This occurs for DESeq2 versions 1.8.1 and 1.10.1.

ADD REPLY • link 9.8 years ago Gregory Warnes ▴ 50

0

Entering edit mode

Relevant source code:

# Proposed model
```{r}
pData$GeneExpr <- pca[,'::PC1::']
pData$HTN <- factor(pData$HTN, ordered=FALSE)

fit <- lm(GeneExpr ~ 1
                    # Technical Covariates
                    + Site
                    + T_Batch
                    + RIN
                    + Bias35
                    #+ Pct_Usable_Bases
                    + Pct_Ribosomal
                    #+ N_Reads
                    #+ N_Ribosomal
                    # Demographics
                    + Gender
                    + Age
                    + Race
                    # Diagnosis Flags
                    + DM
                    + as.factor(HTN)
                    # Medical Covariates
                    #+ BP_systole   # only 185 observations
                    #+ BP_diastole  # only 185 observations
                    + log(BMI, 2)
                    # GFR
                    + GFR
          ,
          data=pData,
          subset=OMIT=="No" & GFR>20,
          na.action=na.omit)
```

```{r}
summary(fit)
```

--------------------------------

# Fit model to all genes

```{r}

modelRows <- rownames(model.matrix(fit))

pData.model <- pData[modelRows, ]
cData.model <- cData[, modelRows]

dds <- DESeqDataSetFromMatrix(countData = cData.model,
                              colData   = pData.model,
                              design    = formula(fit)
                              )

times <- system.time(
  dds <- DESeq(dds)
)
times[3]/3600
```

## Elapsed: 11698 43893 5418

```{r}
fit.deseq <- results(dds,
                     name="GFR",
                     parallel=TRUE
                     )
```

ADD REPLY • link 9.8 years ago Gregory Warnes ▴ 50

0

Entering edit mode

Hi Gregory, Can you try running without the beta prior? DESeq(dds, betaPrior=FALSE). This is a perfectly valid analysis, giving more standard LFCs as one would get from DESeq or other methods. You have many covariates in the design formula, and I'd guess the correlation of some of these causes some problem for the prior estimation.

ADD REPLY • link 9.8 years ago Michael Love 43k

score 0 · Answer 1 · 2015-11-13

0

Entering edit mode

Michael Love 43k

@mikelove

Last seen 21 hours ago

United States

I can't really diagnose without also seeing more of the code you used here. And what is the colData(dds) look like? Note that we changed the DESeq2 approach to interaction terms, they now do not use a prior. You may have to rebuild the 'dds' object using the latest version of DESeq2 to continue.

ADD COMMENT • link 10.0 years ago Michael Love 43k

score 0 · Answer 2 · 2015-11-13

Ok, here is the code I used until the error:

library("DESeq2")
library("edgeR")
samples = read.csv("samples.csv",header=TRUE)
samples$countf = paste(samples$LibraryName, "count", sep=".")
counts = readDGE(samples$countf)$counts
countData<-counts
coldata=with(samples, data.frame(LibraryName=I(LibraryName), Subject=Subject, Condition=Condition, Replicate=Replicate))
coldata$Subject<-factor(coldata$Subject)
coldata$Replicate<-factor(coldata$Replicate)

#Simple model
dds<-DESeqDataSetFromMatrix(countData=countData, colData=coldata, design=~Subject+Condition)

#relevel samples
dds$Condition<-relevel(dds$Condition, "Control")

#collapse technical replicates
ddscollapsed <- collapseReplicates ( dds, groupby = dds$Replicate, renameCols = TRUE)
dds<-ddscollapsed

#test for DE genes
dds<-DESeq(dds)
resultsNames(dds)

> colData(dds)
DataFrame with 20 rows and 4 columns
LibraryName Subject Condition Replicate
<factor> <factor> <factor> <factor>
1 HAP308T_Control_1 1 Control 1
2 HAP307T_Control_1 2 Control 2
3 HAP305T_Control_1 3 Control 3
4 HAP304T_Control_1 4 Control 4
5 HAP303T_Control_1 5 Control 5
... ... ... ... ...
16 HAP302T_HRV1B_1 6 Virus 16
17 CCF1T_HRV1B_1 1 Virus 17
18 CF745T_HRV1B_1 2 Virus 18
19 CF756T_HRV1B_1 3 Virus 19
20 CF757T_HRV1B_1 4 Virus 20