DESeq2 error on analyzing microbiome data
2
0
Entering edit mode
zhigang.li • 0
@zhigangli-21311
Last seen 2.5 years ago

Hi There,

I am trying to do a differential abundance analysis with microbiome sequencing data suing DESeq2 package, but I keep getting errors after trying different approaches within the package. The data is a 50 by 501 matrix with each row being a sample and the first column being the group indicator and the other 500 columns are sequencing reads for 500 taxa. No missing data in the data set. I can upload the data, but I don't know how to do that. Below is the R script for the analysis. Please help. Thanks a lot!

class(data) [1] "data.frame"

# make data in the format for DESeq2

countsData=t(data[,-1])+1

rownames(countsData)=colnames(data)[-1]

class(countsData)

[1] "matrix"

print(dim(countsData))

[1] 500 50

print(countsData[1:5,1:5]

                  [,1] [,2] [,3] [,4] [,5]


rawCount1 8 8 10 11 9

rawCount2 8 9 8 7 8

rawCount3 36 34 30 36 33

rawCount4 7 8 6 5 8

rawCount5 10 7 63 13 64

# get predictor data and annotate the predictor data

xData=as.data.frame(data[,1])

colnames(xData)=colnames(data)[1]

xData[,"x"]=as.factor(xData[,"x"])

class(xData)

[1] "data.frame"

print(dim(xData))

[1] 50 1

print(xData[1:5,])

[1] 0 0 1 0 1

Levels: 0 1

# prepare for analysis

dds <- DESeqDataSetFromMatrix(countData=countsData,colData=xData,design=~x)

# run analysis

suppressWarnings(runDeseq<- try(DESeq(dds, quiet = TRUE), silent = TRUE))

if (inherits(runDeseq, "try-error")) {

# If the parametric fit failed, try the local.

suppressWarnings(runDeseq<- try(DESeq(dds, fitType = "local", quiet = TRUE),
silent = TRUE))

if (inherits(runDeseq, "try-error")) {

# If local fails, try the mean

suppressWarnings(runDeseq<- try(DESeq(dds, fitType = "mean", quiet = TRUE),
silent = TRUE))
}

if (inherits(runDeseq, "try-error")) {

# If still bad, quit with error.

stop("DESeq1 fail")
}


}

When I dig into the error message, it says "simpleError in estimateDispersionsFit(object, fitType = fitType, quiet = quiet): all gene-wise dispersion estimates are within 2 orders of magnitude from the minimum value, and so the standard curve fitting techniques will not work. One can instead use the gene-wise estimates as final estimates: dds <- estimateDispersionsGeneEst(dds) dispersions(dds) <- mcols(dds)$dispGeneEst ...then continue with testing using nbinomWaldTest or nbinomLRT". So I tried dds <- DESeqDataSetFromMatrix(countData=countsData,colData=xData,design=~x) dds <- estimateDispersionsGeneEst(dds) but it gave me another error message: "Error in .local(object, ...) : first calculate size factors, add normalizationFactors, or set normalized=FALSE". Here is some session info: sessionInfo() R version 3.6.1 (2019-07-05) Platform: x86_64-w64-mingw32/x64 (64-bit) Running under: Windows 10 x64 (build 16299) packageVersion("DESeq2") [1] ‘1.24.0’ deseq2 microbiome sequencing • 452 views ADD COMMENT 0 Entering edit mode @mikelove Last seen 9 hours ago United States The second error seems easy to overcome. Did you try following the advice printed there, that is “first estimate size factors...”? See estimateSizeFactors(). ADD COMMENT 0 Entering edit mode zhigang.li • 0 @zhigangli-21311 Last seen 2.5 years ago Thanks a lot for your help. I tried estimateSizeFactors() and changed the code a little bit. Now it is working. Here is my new code for the analysis. Does it look correct to you? Thanks! # run analysis suppressWarnings(runDeseq<- try(DESeq(dds, quiet = TRUE), silent = TRUE)) if (inherits(runDeseq, "try-error")) { # If the parametric fit failed, try the local. suppressWarnings(runDeseq<- try(DESeq(dds, fitType = "local", quiet = TRUE), silent = TRUE)) if (inherits(runDeseq, "try-error")) { # If local fails, try the mean suppressWarnings(runDeseq<- try(DESeq(dds, fitType = "mean", quiet = TRUE), silent = TRUE)) } if (inherits(runDeseq, "try-error")) { dds <- estimateSizeFactors(dds) dds <- estimateDispersionsGeneEst(dds) dispersions(dds) <- mcols(dds)$dispGeneEst

runDeseq <- try(nbinomWaldTest(dds, quiet = TRUE),silent = TRUE)
}


}