Error in estimateSizeFactorsForMatrix(counts.byGenes)

[HeenaBioinfo]

jscs.noNovel <- runJunctionSeqAnalyses(sample.files = countFiles.noNovel,

•                                    sample.names = sample_small$sample,
  

•                                    condition=factor(sample_small$clinGroup_PILOT),
  

•                                    flat.gff.file = gff.file,
  

•                                    nCores = 1,
  

•                                    analysis.type = "junctionsAndExons")
  

STARTING runJunctionSeqAnalyses (v1.14.0) (Wed Jul 29 16:13:26 2020)
rJSA: sample.files: /space/Junction_data/rawCts/S000084.l.r.m.c.lib.g.k2.a/QC.spliceJunctionAndExonCounts.forJunctionSeq.txt.gz, /space/Junction_data/rawCts/S000556.l.r.m3.c.lib.g.k2.a/QC.spliceJunctionAndExonCounts.forJunctionSeq.txt.gz, /space/Junction_data/rawCts/S000556.l.r.m4.c.lib.g.k2.a/QC.spliceJunctionAndExonCounts.forJunctionSeq.txt.gz, /space/Junction_data/rawCts/S000556.l.r.m8.c.lib.g.k2.a/QC.spliceJunctionAndExonCounts.forJunctionSeq.txt.gz
rJSA: sample.names: S000084.l.r.m.c.lib.g.k2.a, S000556.l.r.m3.c.lib.g.k2.a, S000556.l.r.m4.c.lib.g.k2.a, S000556.l.r.m8.c.lib.g.k2.a
rJSA: condition: ERpHER2p, ERpHER2nLNn, ERpHER2nLNn, ERpHER2nLNn
rJSA: analysis.type: junctionsAndExons
rJSA: use.junctions: TRUE
rJSA: use.novel.junctions: TRUE
rJSA: use.exons: TRUE
rJSA: nCores: 1
rJSA: use.covars:
rJSA: test.formula0: ~ sample + countbin
rJSA: test.formula1: ~ sample + countbin + condition:countbin
rJSA: use.multigene.aggregates: FALSE
rJSA: Reading Count files... Wed Jul 29 16:13:26 2020.
-> STARTING readJunctionSeqCounts (Wed Jul 29 16:13:26 2020)
---> RJSC; (v1.14.0)
---> RJSC: samplenames: S000084.l.r.m.c.lib.g.k2.a,S000556.l.r.m3.c.lib.g.k2.a,S000556.l.r.m4.c.lib.g.k2.a,S000556.l.r.m8.c.lib.g.k2.a
---> RJSC: flat.gff.file: /space/Junction_data/flattened.gtf.gz
---> RJSC: use.exons:TRUE
---> RJSC: use.junctions:TRUE
---> RJSC: use.novel.junctions:TRUE
---> File read complete.
---> Extracted counts. Found 669130 features so far.
---> Extracted gene-level counts. Found: 47257 genes and aggregate-genes.
---> Removed gene features. Found: 621873 features to be included so far.
---> Note: 524335 counting bins from overlapping genes
---> There are 16295 multigene aggregates.
---> There are 73216 genes that are part of an aggregate.
---> Removed multigene-aggregate features. Found: 97538 features to be included so far.
---> Final feature count: 97538 features to be included in the analysis.
---> Extracted feature counts.
---> counts complete.
-----> reading annotation...
-----> formatting annotation...
-----> initial generation...
-----> creating jscs...
-----> generating count vectors... (Wed Jul 29 16:14:13 2020)
Using single-core execution.
getAllJunctionSeqCountVectors: dim(counts) = 97538,4 (Wed Jul 29 16:14:13 2020)
getAllJunctionSeqCountVectors: dim(gct) = 47257,4
getAllJunctionSeqCountVectors: out generated. dim = 97538,8 (Wed Jul 29 16:14:16 2020)
-----> count vectors generated (Wed Jul 29 16:14:16 2020)
-----> generating DESeqDataSet... (Wed Jul 29 16:14:16 2020)
-----> DESeqDataSet generated (Wed Jul 29 16:14:16 2020)
rJSA: Count files read. Wed Jul 29 16:14:16 2020.
rJSA: Estimating Size Factors... Wed Jul 29 16:14:16 2020.
Error in estimateSizeFactorsForMatrix(counts.byGenes) :
every gene contains at least one zero, cannot compute log geometric means

Hi Can you please guide me through this error?

[HeenaBioinfo]

I am not able to generate the size factors, please guide..

[hartleys]

This occurs commonly in large samples and/or single-cell data. It can also occur if the samples are extremely variable and very different from one another (like in single-cell data). The geometric size factor calculation used by DESeq simply fails in this situation.

I would recommend skipping the size factor calculation step and replacing it with your own. For large samples / single sample data the Upper Quartile method seems to be gaining traction as the industry standard.

I would use the step-by-step analysis pipeline from the manual, but replace the estimateJunctionSeqSizeFactors step with the following:

#replace the estimation function with a version that takes an external size factor. 
# I Should have added this long ago...
estimateJunctionSeqSizeFactors.forced <- function( jscs , sizeFactors, verbose = FALSE){
   stopifnot( is( jscs, "JunctionSeqCountSet") )
   sizeFactors(jscs) <- sizeFactors
   sizeFactors(jscs@DESeqDataSet) <- rep(sizeFactors(jscs),2)
   fData(jscs)$baseMean <- rowMeans(counts(jscs, normalized= T))
   fData(jscs)$baseVar <- rowVars(counts(jscs, normalized=T))
   jscs@altSizeFactors <- data.frame()
   return(jscs);                         
}

#Calculate size factors using your favorite method.
#I recommend the upper quartile method, shown below:
SF <- sapply(1:ncol(jscs@geneCountData),function(i){
   quantile(jscs@geneCountData[,i],p=0.75);
})
#Now copy in those size factors:
jscs <- estimateJunctionSeqSizeFactors.forced(jscs,sizeFactors=SF)

For the future I'll fold this into the codebase with the next update.