DESeq2 normFactors method: Error in fitBetaWrapper
1
0
Entering edit mode
Jakub ▴ 50
@jakub-9073
Last seen 9 weeks ago
United Kingdom

Dear all,

this is a follow up from my last question. When my run my DESeq2 analysis using the usual estimateSizeFactors, it runs fine, however if I provide a normalisation matrix, I get the following error on estimateDispersions:

Error in fitBetaWrapper(ySEXP = counts(object), xSEXP = modelMatrix, nfSEXP = normalizationFactors,  : in call to fitBeta, the following arguments contain NA: alpha_hatSEXP

I have generated the matrix & checked for the absence of NA and 0 & run DESeq2 as follows:

data <- read.table(file = "PATH.tsv", header=TRUE, row.name=1)
normFactors <- data.matrix(data)
normFactors <- normFactors / exp(rowMeans(log(normFactors)))
normFactors[is.na(normFactors)] <- 1
dds = DESeqDataSetFromMatrix(experiment, design, design = ~ genotype)
normalizationFactors(dds) <- normFactors
dds <- estimateDispersions(dds)

My design file shows that there are >3 replicates per group

   track group    batch genotype
1  A.12.R1 A-12   2nd       A
2  A.12.R2 A-12   2nd       A
3  A.12.R3 A-12   2nd       A
4  A.12.R4 A-12   2nd       A
5  B.12.R1 B-12   2nd       B
6  B.12.R2 B-12   2nd       B
7  B.12.R3 B-12   2nd       B 
8  B.12.R4 C-12   2nd       C  
9  C.12.R1 C-12   2nd       C  
10 C.12.R2 A-12   2nd       A 
11 C.12.R3 A-12   2nd       A  
12 C.12.R4 C-12   2nd       C   

Many thanks, Jakub

session.Info():

R version 3.2.3 (2015-12-10)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Scientific Linux release 6.7 (Carbon)

locale:
[1] C

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

other attached packages:
 [1] biomaRt_2.26.1                          
 [2] geneplotter_1.48.0                      
 [3] annotate_1.48.0                         
 [4] XML_3.98-1.4                            
 [5] lattice_0.20-33                         
 [6] MatchIt_2.4-21                          
 [7] MASS_7.3-45                             
 [8] TxDb.Mmusculus.UCSC.mm10.knownGene_3.2.2
 [9] TxDb.Mmusculus.UCSC.mm10.ensGene_3.2.2  
[10] GenomicFeatures_1.22.13                 
[11] AnnotationDbi_1.32.3                    
[12] VennDiagram_1.6.16                      
[13] futile.logger_1.4.1                     
[14] goseq_1.22.0                            
[15] RSQLite_1.0.0                           
[16] DBI_0.3.1                               
[17] geneLenDataBase_1.6.0                   
[18] BiasedUrn_1.07                          
[19] vsn_3.38.0                              
[20] dynamicTreeCut_1.63-1                   
[21] ggplot2_2.1.0                           
[22] gplots_3.0.1                            
[23] RColorBrewer_1.1-2                      
[24] DESeq2_1.10.1                           
[25] RcppArmadillo_0.6.600.4.0               
[26] Rcpp_0.12.4                             
[27] SummarizedExperiment_1.0.2              
[28] Biobase_2.30.0                          
[29] GenomicRanges_1.22.4                    
[30] GenomeInfoDb_1.6.3                      
[31] IRanges_2.4.8                           
[32] S4Vectors_0.8.11                        
[33] BiocGenerics_0.16.1                     

loaded via a namespace (and not attached):
 [1] locfit_1.5-9.1          GO.db_3.2.2             Rsamtools_1.22.0       
 [4] Biostrings_2.38.4       gtools_3.5.0            plyr_1.8.3             
 [7] futile.options_1.0.0    acepack_1.3-3.3         BiocInstaller_1.20.1   
[10] zlibbioc_1.16.0         gdata_2.17.0            Matrix_1.2-4           
[13] rpart_4.1-10            preprocessCore_1.32.0   splines_3.2.3          
[16] BiocParallel_1.4.3      foreign_0.8-66          RCurl_1.95-4.8         
[19] munsell_0.4.3           rtracklayer_1.30.2      mgcv_1.8-12            
[22] nnet_7.3-12             gridExtra_2.2.1         Hmisc_3.17-3           
[25] GenomicAlignments_1.6.3 bitops_1.0-6            nlme_3.1-126           
[28] xtable_1.8-2            gtable_0.2.0            affy_1.48.0            
[31] scales_0.4.0            KernSmooth_2.23-15      XVector_0.10.0         
[34] genefilter_1.52.1       affyio_1.40.0           limma_3.26.8           
[37] latticeExtra_0.6-28     Formula_1.2-1           lambda.r_1.1.7         
[40] survival_2.38-3         colorspace_1.2-6        cluster_2.0.3          
[43] caTools_1.17.1         
deseq2 normalization • 1.4k views
ADD COMMENT
0
Entering edit mode

It does indeed, sorry for not being explicit about this! Thanks again.

ADD REPLY
2
Entering edit mode
@mikelove
Last seen 6 hours ago
United States
What values are in normFactors right before you assign it to the slot in dds, what is the minimum, maximum, are all values finite?
ADD COMMENT
0
Entering edit mode

Many thanks!

Indeed, not all values were finite (lots of dividing by zero). I have added the line:

normFactors[!is.finite(normFactors)] <- 1

 A typical column now looks like this:

Min.   :0.1656  
1st Qu.:0.9380  
Median :1.0000  
Mean   :1.0208  
3rd Qu.:1.1189  
Max.   :4.2724
ADD REPLY
0
Entering edit mode
And does this resolve the issue?
ADD REPLY

Login before adding your answer.

Traffic: 836 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6