I have been using the Kallisto -> Sleuth pipeline on my RNAseq data consisting of a control and two treatments (in triplicates). From that I got 8 DETs (q-value < 0.05) for treatment 1 and 12 DETs for treatment 2. The highest beta ("equivalent" to log fold change) were around 6 and the lowest -4.
In order to validate my results from Sleuth, I wanted to use DESeq2 to see if I got similar results (+ make "prettier"/customizable figures). Tximport were used to load the count data from Kallisto into the DESeq2 pipeline. The first runs with DESeq2 were pretty consistent with the Sleuth results (considering Sleuth is more conservative than DESeq2) with 19 DETs for treatment1 and 34 for treatment 2. The highest/lowest log fold change was around 6.6 and -3.
BUT then... I updated by version of R and DESeq2 and suddenly I got completely different results. The exact same R script (with updated packages) have then 100 DETs for treatment1 and 143 for treatment2 with highest/lowest log fold changes around 21 and -25.5.
This is really weird - have others of you had similar experiences? What causes this? and how do I "fix" it (especially the log fold changes seems way off in the second analysis)
#Libraries library(tximport) library(rhdf5) library(DESeq2) #base directory: dir <- "~/Documents/RNASeq/kallisto/data/" #Metadata samples <- read.table(file.path(dir, "metadata.txt"), header=TRUE, stringsAsFactors = FALSE) #Kallisto-data files <- file.path(dir,"results", samples$path, "abundance.tsv") names(files)<- paste0(samples$Sample) txi <- tximport(files, type = "kallisto", txOut = TRUE) #DESeq2 sampleTable <- data.frame(condition = factor(rep(c("Treatment1", "Treatment2", "Control"), each = 3))) rownames(sampleTable) <- colnames(txi$counts) dds <- DESeqDataSetFromTximport(txi, sampleTable, ~condition) #Pre-filtering dds <- dds[rowSums(counts(dds))>1, ] #Control = first level dds$condition <- relevel(dds$condition, ref="Control") #Run DESeq2 dds <-DESeq(dds) #Call and inspect the results table res <- results(dds) # results for each treatment vs. control res_treat1 <- results(dds, contrast=c("condition", "Treatment1", "Control")) res_treat2 <- results(dds, contrast=c("condition", "Treatment2", "Control")) # Call information on the result object on which variables and tests were used: mcols(res_treat1, use.names=TRUE) mcols(res_treat2, use.names=TRUE) # Call number of significant results with BH-adjusted p-value < 0.05 sum(res_treat1$padj < 0.05, na.rm=TRUE) sum(res_treat2$padj < 0.05, na.rm=TRUE) #Subset the results and sort them by the log2 fold change down-regulation: resSig_treat1<- res[which(res_treat1$padj < 0.05), ] head(resSig_treat1[ order(resSig_treat1$log2FoldChange ), ]) resSig_treat2<- res[which(res_treat2$padj < 0.05), ] head(resSig_treat2[ order(resSig_treat2$log2FoldChange ), ]) #And sort by up-regulations: tail(resSig_treat1[order(resSig_treat1$log2FoldChange ), ] ) tail(resSig_treat2[order(resSig_treat2$log2FoldChange ), ] )
Session-info (for the updated run):
R version 3.4.0 (2017-04-21) Platform: x86_64-apple-darwin15.6.0 (64-bit) Running under: macOS Sierra 10.12.5 Matrix products: default BLAS: /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib LAPACK: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRlapack.dylib locale:  da_DK.UTF-8/da_DK.UTF-8/da_DK.UTF-8/C/da_DK.UTF-8/da_DK.UTF-8 attached base packages:  parallel stats4 stats graphics grDevices utils datasets methods base other attached packages:  DESeq2_1.17.2 SummarizedExperiment_1.7.4 DelayedArray_0.3.9 matrixStats_0.52.2  Biobase_2.37.2 GenomicRanges_1.29.4 GenomeInfoDb_1.13.4 IRanges_2.11.3  S4Vectors_0.15.3 BiocGenerics_0.23.0 rhdf5_2.21.1 tximport_1.5.0 loaded via a namespace (and not attached):  genefilter_1.59.0 locfit_1.5-9.1 splines_3.4.0 lattice_0.20-35  colorspace_1.3-2 htmltools_0.3.6 base64enc_0.1-3 survival_2.41-3  XML_3.98-1.7 rlang_0.1.1 DBI_0.6-1 foreign_0.8-68  BiocParallel_1.11.2 RColorBrewer_1.1-2 GenomeInfoDbData_0.99.1 plyr_1.8.4  stringr_1.2.0 zlibbioc_1.23.0 munsell_0.4.3 gtable_0.2.0  htmlwidgets_0.8 memoise_1.1.0 latticeExtra_0.6-28 knitr_1.16  geneplotter_1.55.0 AnnotationDbi_1.39.1 htmlTable_1.9 Rcpp_0.12.11  acepack_1.4.1 xtable_1.8-2 scales_0.4.1 backports_1.1.0  checkmate_1.8.2 Hmisc_4.0-3 annotate_1.55.0 XVector_0.17.0  gridExtra_2.2.1 ggplot2_2.2.1 digest_0.6.12 stringi_1.1.5  grid_3.4.0 tools_3.4.0 bitops_1.0-6 magrittr_1.5  RSQLite_1.1-2 lazyeval_0.2.0 RCurl_1.95-4.8 tibble_1.3.3  Formula_1.2-1 cluster_2.0.6 Matrix_1.2-10 data.table_1.10.4  rpart_4.1-11 nnet_7.3-12 compiler_3.4.0