I have recently been trying to analyse 12 non-replicated metagenome samples using DESeq2 as an exploratory analysis to find genes with differing abundance to allow for inferences to be obtained with respect to GO terms and KEGG pathways associated with these genes. I have formatted a counts table from a Bowtie2 mapping into HTSeq gene region count extraction procedure which I can successfully read into DESeq2. However, when I attempt to run the program, I receive the following error (note that I'm running it in parallel, and my data set is very, very large):
mean-dispersion relationship final dispersion estimates, MLE betas: 12 workers fitting model and testing: 12 workers -- replacing outliers and refitting for 2151512 genes -- DESeq argument 'minReplicatesForReplace' = 7 =-- original counts are preserved in counts(dds) estimating dispersions Error in estimateDispersionsGeneEst(objectSub, quiet = quiet, modelMatrix = modelMatrix) : the number of samples and the number of model coefficients are equal, i.e., there are no replicates to estimate the dispersion. use an alternate design formula Calls: DESeq -> refitWithoutOutliers -> estimateDispersionsGeneEst In addition: Warning message: In checkForExperimentalReplicates(object, modelMatrix) : same number of samples and coefficients to fit, estimating dispersion by treating samples as replicates. read the ?DESeq section on 'Experiments without replicates' Execution halted
The R output including headers for my input files as well as the sessionInfo are provided below:
R version 3.2.5 (2016-04-14) -- "Very, Very Secure Dishes" Copyright (C) 2016 The R Foundation for Statistical Computing Platform: x86_64-pc-linux-gnu (64-bit)
R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details.
Natural language support but running in an English locale
R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R.
> ############## START > > # Load deseq2 and set working directory > library(DESeq2) > library(pheatmap) > library(RColorBrewer) > library(BiocParallel) > setwd("/home/workingdir") > > register(MulticoreParam(12)) > > # Load in data > countData = as.matrix(read.table(file='combined.htseq.counts', header = TRUE, sep = "\t", stringsAsFactors = FALSE, row.names="gene_id")) > head(countData) AR DW GH1 LS OF STP1 STP2 STP11A STP12B STP13A STP15C WS OF0002336 0 8 0 0 1744 0 0 0 0 0 0 0 OF0057878 0 0 0 0 1560 0 6 2 4 0 0 0 OF0008217 0 0 0 0 1004 0 0 0 0 0 0 0 STP11A0014293 0 0 2 0 33 2 378 697 1674 294 474 0 LS0032436 0 0 0 981 0 0 0 0 0 1 0 149 OF0055859 0 0 0 0 892 0 0 0 1 1 0 0 > > colData = read.csv('colData.csv') > head(colData) condition 1 AR 2 DW 3 GH1 4 LS 5 OF 6 STP1 > > # Create DESeq2 data set > dds <- DESeqDataSetFromMatrix(countData = countData, colData = colData, design = ~ condition) > dds <- dds[ rowSums(counts(dds)) > 1, ] # 1.3.6 Pre-filtering - removes genes for which only 0 or 1 reads are present > > # Reduce RAM usage > rm(countData) > rm(colData) > > # Session info > sessionInfo() R version 3.2.5 (2016-04-14) Platform: x86_64-pc-linux-gnu (64-bit) Running under: SUSE Linux Enterprise Server 11 SP4
locale:  LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C  LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8  LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8  LC_PAPER=en_US.UTF-8 LC_NAME=C  LC_ADDRESS=C LC_TELEPHONE=C  LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
attached base packages:  parallel stats4 stats graphics grDevices utils datasets  methods base
other attached packages:  BiocParallel_1.4.3 RColorBrewer_1.1-2  pheatmap_1.0.8 DESeq2_1.10.1  RcppArmadillo_0.7.500.0.0 Rcpp_0.12.8  SummarizedExperiment_1.0.2 Biobase_2.30.0  GenomicRanges_1.22.4 GenomeInfoDb_1.6.3  IRanges_2.4.8 S4Vectors_0.8.11  BiocGenerics_0.16.1
loaded via a namespace (and not attached):  genefilter_1.52.1 locfit_1.5-9.1 splines_3.2.5  lattice_0.20-33 colorspace_1.3-1 htmltools_0.3.5  chron_2.3-47 survival_2.40-1 XML_3.98-1.5  foreign_0.8-66 DBI_0.5-1 lambda.r_1.1.9  plyr_1.8.4 stringr_1.1.0 zlibbioc_1.16.0  munsell_0.4.3 gtable_0.2.0 futile.logger_1.4.3  latticeExtra_0.6-28 knitr_1.15.1 geneplotter_1.48.0  AnnotationDbi_1.32.3 htmlTable_1.7 acepack_1.4.1  xtable_1.8-2 scales_0.4.1 Hmisc_4.0-0  annotate_1.48.0 XVector_0.10.0 gridExtra_2.2.1  ggplot2_2.2.0 digest_0.6.10 stringi_1.1.2  grid_3.2.5 tools_3.2.5 magrittr_1.5  lazyeval_0.2.0 tibble_1.2 RSQLite_1.0.0  Formula_1.2-1 cluster_2.0.5 futile.options_1.0.0  Matrix_1.2-4 data.table_1.9.6 assertthat_0.1  rpart_4.1-10 nnet_7.3-12 > > # Perform the differential expression (1.4) > dds <- DESeq(dds, parallel=TRUE)
^^ This is where the code dies
If someone could help me figure out why DESeq2 is not allowing me to run without replicates (as I understand it is capable of despite encouraging you not to take any results as conclusive) that would be appreciated. Thanks.