estimateEffectSizes: Error in if (any(dispersion < 0)) stop("dispersion values must be non-negative")

[ashish674675]

Hi
I am trying to run JunctionSeq on 15 samples with 5 conditions and three replicates each. Conditions are "wild type" and four mutants. I want to make pairwise compaisons between wild type and each mutant. My data comes from a single ended, stranded (fr-secondstrand) library. I keep running into this error in the last step done by estimateEffectSizes . This error comes for a lot of genes.

Error in if (any(dispersion < 0)) stop("dispersion values must be non-negative") : 
  missing value where TRUE/FALSE needed
In addition: Warning message:
In FUN(X[[i]], ...) : glmnb.fit failed for AT5G09885:

I first tried running runJunctionSeqAnalyses on all samples together and on pairs separately but the same error comes. Using section 5.2 of the JunctionSeq vignette I ran it step by step and this error comes in the last step for all samples.

I will really appreciate if you can help me troubleshoot this problem. I am using version 1.17.0 which I installed from bioconductor.
Thank You

decoder <- read.table(
  "decoder_file",
  header=T,stringsAsFactors=F)

wt_m6_dec <- decoder[1:6, ]
design6 <- data.frame(condition = factor(wt_m6_dec$group.ID));

countFiles6 <- paste0(
  "novel_and_known_splice_sites/",
  wt_m6_dec$sample.ID,
  "/QC.spliceJunctionAndExonCounts.withNovel.forJunctionSeq.txt.gz"
);
jscs6 = readJunctionSeqCounts(countfiles = countFiles6,
                             samplenames = wt_m6_dec$sample.ID,
                             design = design6,
                             flat.gff.file = "novel_and_known_splice_sites/withNovel.forJunctionSeq.gff.gz"
);
jscs6 <- estimateJunctionSeqSizeFactors(jscs6);
jscs6 <- estimateJunctionSeqDispersions(jscs6, nCores = 1);

# not using method.dispFit = local gives "unable to fit model" error
jscs6 <- fitJunctionSeqDispersionFunction(jscs6, method.dispFit = "local");
jscs6 <- testForDiffUsage(jscs6, nCores = 1);
jscs6 <- estimateEffectSizes(jscs6, nCores = 1)

[hartleys]

This will occur if ALL of the genes list a dispersion of NA. ie: the model dispersion fit always fails.

Have you tried running the analysis on a subset? For example: just the WT vs one of the mutants?

[ashish674675]

Yes I ran with only WT vs one mutant but got the same error. I have 4 mutants and I get the same error in every possible pairwise comparison. If I ignore this error and move to the next step of writing and plotting results, I do not get any errors and JunctionSeq gives me around 700 significant genes. I was expecting many of these genes to be significant.
I also tried subsetting the jscs considering it an ExpressionSet object but then I get unmatched dimensions error,