Hi,

I have 3 replicates of 4 different samples sequenced by RNA-Seq (Untreated_GFP(-), Untreated_GFP(+), Treated_GFP(-) and Treated_GFP(+)) and I would like compare in 4 different ways for exon usage using DEXSeq;

(1) Untreated_GFP(+) vs Untreated_GFP(-) (untGFP-P vs untGFP-N)

(2) Treated_GFP(+) vs Treated_GFP(-) (trtGFP-P vs trtGFP-N)

(3) Treated_GFP(+) vs Untreated_GFP(-) (trtGFP-P vs untGFP-N)

(4) Treated_GFP(-) vs Untreated_GFP(-) (trtGFP-N vs untGFP-N).

Weeks ago I ran the script below and it worked fine. However, I upgraded Linux and then try the same code and I have met a problem while running estimateExonFoldChanges, with error:

dxd = estimateExonFoldChanges(dxd,fitExpToVar="condition", BPPARAM = mP)
Error in exp(alleffects) : non-numeric argument to mathematical function.

I tried DEXSeq Error when calculating Exon Fold-changes. But it did not help.

I want to first estimateSizeFactors and then do analyses for the subsets. If I do analyses for conditions one by one, it works.

thanks in advance.

ilyas.

my script is as below.

suppressPackageStartupMessages(library("DEXSeq"))
library(BiocParallel)
curDir = getwd()
setwd(curDir)

nWorkers <- 8
if Sys.info()[1][[1]] == "Windows"){
    mP <- SnowParam(workers = nWorkers)
}else{
    mP = MulticoreParam(workers = nWorkers)
}

flattenedFile = "Danio_rerio.GRCz10.82.gff"
allConditions <- c("trtGFP-N_R1","trtGFP-N_R2","trtGFP-N_R3","trtGFP-P_R1","trtGFP-P_R2","trtGFP-P_R3","untGFP-N_R1","untGFP-N_R2","untGFP-N_R3","untGFP-P_R1","untGFP-P_R2","untGFP-P_R3")
allComparisons <- c("trtGFP-P_vs_untGFP-P","trtGFP-N_vs_untGFP-N","trtGFP-P_vs_trtGFP-N","untGFP-P_vs_untGFP-N")
sampleTable= data.frame(row.names= allConditions, condition = c(rep("trtGFP-N",3),rep("trtGFP-P",3),rep("untGFP-N",3),rep("untGFP-P",3)),libType = c(rep("single-end",12)))
countFiles = list.files(paste(curDir,"/counts",sep = ""),pattern="_NH.sam.txt$",full.names=TRUE)
dxdFull = DEXSeqDataSetFromHTSeq(countFiles, sampleData = sampleTable, design= ~ sample + exon + condition:exon, flattenedfile = flattenedFile)

dxdFull = estimateSizeFactors(dxdFull) #normalise first and then do the further analyses!!! estimateSizeFactors first or later?

for (i in 1:length(allComparisons)){
  output_folder <- allComparisons[i]
  dir.create(output_folder ,showWarnings = FALSE)
  conditionSplit <- strsplit(output_folder,"_vs_")
  trt_sp = conditionSplit[[1]][1]
  unt_sp = conditionSplit[[1]][2]

  dxd <- dxdFull[,colData(dxdFull)$condition %in% c(trt_sp,unt_sp)]
  dxd$condition <- droplevels(dxd$condition)
  dxd$condition <- relevel(dxd$condition,unt_sp)

  dxd= estimateDispersions(dxd, BPPARAM = mP)
    
  tiff(paste(output_folder,"/",paste(output_folder,"_plotDispersions.tiff", sep = ""), sep = ""))
  plotDispEsts(dxd)
  dev.off()

  dxd=testForDEU(dxd, BPPARAM = mP)

  dxd = estimateExonFoldChanges(dxd,fitExpToVar="condition", BPPARAM = mP)

  dxr1 = DEXSeqResults(dxd)
  tiff(paste(output_folder,"/",paste(output_folder,"_plotMA.tiff", sep = ""), sep = ""))
  plotMA(dxr1,cex=0.8)
  dev.off()

  DEXSeqHTML( dxr1, FDR=0.1, color=c("#FF000080", "#0000FF80") )

}

> SessionInfo()

R version 3.2.2 (2015-08-14)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 15.10

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C               LC_TIME=de_DE.UTF-8        LC_COLLATE=en_US.UTF-8     LC_MONETARY=de_DE.UTF-8   
 [6] LC_MESSAGES=en_US.UTF-8    LC_PAPER=de_DE.UTF-8       LC_NAME=C                  LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=de_DE.UTF-8 LC_IDENTIFICATION=C       

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

other attached packages:
 [1] DEXSeq_1.16.4              DESeq2_1.10.0              RcppArmadillo_0.6.300.2.0  Rcpp_0.12.2                SummarizedExperiment_1.0.1
 [6] GenomicRanges_1.22.1       GenomeInfoDb_1.6.1         IRanges_2.4.4              S4Vectors_0.8.3            Biobase_2.30.0            
[11] BiocGenerics_0.16.1        BiocParallel_1.4.0        

loaded via a namespace (and not attached):
 [1] genefilter_1.52.0    statmod_1.4.22       locfit_1.5-9.1       reshape2_1.4.1       splines_3.2.2        lattice_0.20-33     
 [7] colorspace_1.2-6     survival_2.38-3      XML_3.98-1.3         foreign_0.8-66       DBI_0.3.1            RColorBrewer_1.1-2  
[13] lambda.r_1.1.7       plyr_1.8.3           stringr_1.0.0        zlibbioc_1.16.0      Biostrings_2.38.2    munsell_0.4.2       
[19] gtable_0.1.2         futile.logger_1.4.1  hwriter_1.3.2        latticeExtra_0.6-26  geneplotter_1.48.0   biomaRt_2.26.1      
[25] AnnotationDbi_1.32.0 proto_0.3-10         acepack_1.3-3.3      xtable_1.8-0         scales_0.3.0         Hmisc_3.17-0        
[31] annotate_1.48.0      XVector_0.10.0       Rsamtools_1.22.0     gridExtra_2.0.0      ggplot2_1.0.1        digest_0.6.8        
[37] stringi_1.0-1        grid_3.2.2           tools_3.2.2          bitops_1.0-6         magrittr_1.5         RCurl_1.95-4.7      
[43] RSQLite_1.0.0        Formula_1.2-1        cluster_2.0.3        futile.options_1.0.0 MASS_7.3-45          rpart_4.1-10        
[49] nnet_7.3-11 

Hi Mehmet, 

Thanks for reporting this. I introduced a bug while trying to fix an old bug, and this was the reason it was giving you that error. Sorry about that, it should work if you install the latest version of DEXSeq (1.16.5).

Alejandro