Question: Variance stabilizing transformation seems not working on NGS (bacterial 16srRNA based) data
0
gravatar for Michael Love
5.6 years ago by
Michael Love24k
United States
Michael Love24k wrote:
hi Shucong, I'm replying and including the Bioconductor list. It depends on what you mean by 'not working', I see that the plots of standard deviation over mean are not flat using the two transformations from DESeq2. For both transformations, you are specifying blind=FALSE, which means that the group means are taken into consideration when estimating the dispersion-mean relationship. Then for the plot, you are calculating the standard deviation over all the samples, regardless of groups. I would say that the plot does not imply that the transformation did not work, as the dispersion-mean relationship was not calculated with respect to inter-group differences because 'blind' was set to FALSE. I can't tell if the manuscript you reference used blind dispersion estimation or not when they were using the VST (I couldn't find a link to supplemental file 3). hope this helps, Mike On Wed, Nov 27, 2013 at 3:21 PM, Shucong Li <shucong.li@umanitoba.ca> wrote: > Hello Dr. Love, > > I tried to used DESeq2 to do VST transformation with my NGS (bacterial > 16srRNA based) data. I came to this approach after I read the manuscript > written by Profs McMurdie and Holmes titled "Waste Not, Want Not: Why > Rarefying Microbiome Data Is Inadmissible". However, the VST > transformation did work as what suggested in the Vignettes. I started > wondering this VST from DESeq2 may not work to my data. I emailed Prof. > Holmes this morning to get her opinion. By the mean time, it would be great > i can have your opinion regarding this issue. > > I also attached the R codes and data (R workspace file). > > Thank you in advance. > > Shucong > > > > > Shucong Li PhD Research Associate > ----------------------------------- > Department of Animal Science > University of Manitoba > 12 Dafoe Rd, Winnipeg, Manitoba > R3T 2N2, Canada > > Tel: (204) 474 9841 > Fax: +1(204) 474 7628 > Email: shucong.li@ad.umanitoba.ca > ----------------------------------- > > > > dds <- DESeq(dds, fitType="local”) > > > > # ddsClean <- replaceOutliersWithTrimmedMean(dds) > > # ddsClean <- DESeq(ddsClean,fitType="local") > > > > plotDispEsts(dds) > > > > rld <- rlogTransformation(dds, blind=FALSE) > > vsd <- varianceStabilizingTransformation(dds, blind=FALSE) > > > > > > library("vsn") > > par(mfrow=c(1,3)) > > notAllZero <- (rowSums(counts(dds))>0) > > meanSdPlot(log2(counts(dds,normalized=TRUE)[notAllZero,] + 1), > > ylim = c(0,2.5)) > > meanSdPlot(assay(rld[notAllZero,]), ylim = c(0,2.5)) > > meanSdPlot(assay(vsd[notAllZero,]), ylim = c(0,2.5)) > > > > > > par(mfrow=c(1,2)) > > dist <- dist(t(assay(rld))) > > fit <- hclust(dist, method="average") > > plot(fit, labels=dds@colData$treat) > > > > dists <- dist(t(assay(vsd))) > > fit <- hclust(dists, method="average") > > plot(fit, labels=dds@colData$treat) > > > > > [[alternative HTML version deleted]]
microbiome deseq2 • 1.7k views
ADD COMMENTlink modified 2.9 years ago by Wolfgang Huber13k • written 5.6 years ago by Michael Love24k
Answer: Variance stabilizing transformation seems not working on NGS (bacterial 16srRNA
0
gravatar for Shucong Li
5.6 years ago by
Shucong Li60
Canada
Shucong Li60 wrote:
Hi Michael, The supplemental file 3 in the manuscript is not available yet. I was told it will be available in a week or so. I will try to rerun the test with blind set as TRUE. Shucong On Nov 27, 2013, at 3:13 PM, Michael Love <michaelisaiahlove@gmail.com> wrote: > hi Shucong, > > I'm replying and including the Bioconductor list. > > It depends on what you mean by 'not working', I see that the plots of standard deviation over mean are not flat using the two transformations from DESeq2. For both transformations, you are specifying blind=FALSE, which means that the group means are taken into consideration when estimating the dispersion-mean relationship. Then for the plot, you are calculating the standard deviation over all the samples, regardless of groups. I would say that the plot does not imply that the transformation did not work, as the dispersion-mean relationship was not calculated with respect to inter-group differences because 'blind' was set to FALSE. I can't tell if the manuscript you reference used blind dispersion estimation or not when they were using the VST (I couldn't find a link to supplemental file 3). > > hope this helps, > > Mike > > On Wed, Nov 27, 2013 at 3:21 PM, Shucong Li <shucong.li@umanitoba.ca> wrote: > Hello Dr. Love, > > I tried to used DESeq2 to do VST transformation with my NGS (bacterial 16srRNA based) data. I came to this approach after I read the manuscript written by Profs McMurdie and Holmes titled "Waste Not, Want Not: Why Rarefying Microbiome Data Is Inadmissible". However, the VST transformation did work as what suggested in the Vignettes. I started wondering this VST from DESeq2 may not work to my data. I emailed Prof. Holmes this morning to get her opinion. By the mean time, it would be great i can have your opinion regarding this issue. > > I also attached the R codes and data (R workspace file). > > Thank you in advance. > > Shucong > > > > > Shucong Li PhD Research Associate > ----------------------------------- > Department of Animal Science > University of Manitoba > 12 Dafoe Rd, Winnipeg, Manitoba > R3T 2N2, Canada > > Tel: (204) 474 9841 > Fax: +1(204) 474 7628 > Email: shucong.li@ad.umanitoba.ca > ----------------------------------- > > > > dds <- DESeq(dds, fitType="local”) > > > > # ddsClean <- replaceOutliersWithTrimmedMean(dds) > > # ddsClean <- DESeq(ddsClean,fitType="local") > > > > plotDispEsts(dds) > > > > rld <- rlogTransformation(dds, blind=FALSE) > > vsd <- varianceStabilizingTransformation(dds, blind=FALSE) > > > > > > library("vsn") > > par(mfrow=c(1,3)) > > notAllZero <- (rowSums(counts(dds))>0) > > meanSdPlot(log2(counts(dds,normalized=TRUE)[notAllZero,] + 1), > > ylim = c(0,2.5)) > > meanSdPlot(assay(rld[notAllZero,]), ylim = c(0,2.5)) > > meanSdPlot(assay(vsd[notAllZero,]), ylim = c(0,2.5)) > > > > > > par(mfrow=c(1,2)) > > dist <- dist(t(assay(rld))) > > fit <- hclust(dist, method="average") > > plot(fit, labels=dds@colData$treat) > > > > dists <- dist(t(assay(vsd))) > > fit <- hclust(dists, method="average") > > plot(fit, labels=dds@colData$treat) [[alternative HTML version deleted]]
ADD COMMENTlink written 5.6 years ago by Shucong Li60
Answer: Variance stabilizing transformation seems not working on NGS (bacterial 16srRNA
0
gravatar for Shucong Li
5.6 years ago by
Shucong Li10
Shucong Li10 wrote:
An embedded and charset-unspecified text was scrubbed... Name: not available URL: <https: stat.ethz.ch="" pipermail="" bioconductor="" attachments="" 20131127="" c2602b76="" attachment.asc="">
ADD COMMENTlink written 5.6 years ago by Shucong Li10
Hi Michael, I have rerun the test with 'blind' set as TRUE. The graph didn't change. It is interesting that rlogTransformation seems to work very well. Please allow me ask one more question about differential analysis with DESeq(). I also run the differential analysis with the same set of data. Original counts was used as suggested. I'm wondering if any variance stabilizing transformation is involved in DESeq() procedure? Shucong On Nov 27, 2013, at 3:31 PM, Shucong Li <shucong.li at="" umanitoba.ca=""> wrote: > Hi Michael, > > The supplemental file 3 in the manuscript is not available yet. I was told it will be available in a week or so. I will try to rerun the test with blind set as TRUE. > > > Shucong > On Nov 27, 2013, at 3:13 PM, Michael Love <michaelisaiahlove at="" gmail.com=""> wrote: > >> hi Shucong, >> >> I'm replying and including the Bioconductor list. >> >> It depends on what you mean by 'not working', I see that the plots of standard deviation over mean are not flat using the two transformations from DESeq2. For both transformations, you are specifying blind=FALSE, which means that the group means are taken into consideration when estimating the dispersion-mean relationship. Then for the plot, you are calculating the standard deviation over all the samples, regardless of groups. I would say that the plot does not imply that the transformation did not work, as the dispersion-mean relationship was not calculated with respect to inter-group differences because 'blind' was set to FALSE. I can't tell if the manuscript you reference used blind dispersion estimation or not when they were using the VST (I couldn't find a link to supplemental file 3). >> >> hope this helps, >> >> Mike >> >> On Wed, Nov 27, 2013 at 3:21 PM, Shucong Li <shucong.li at="" umanitoba.ca=""> wrote: >> Hello Dr. Love, >> >> I tried to used DESeq2 to do VST transformation with my NGS (bacterial 16srRNA based) data. I came to this approach after I read the manuscript written by Profs McMurdie and Holmes titled "Waste Not, Want Not: Why Rarefying Microbiome Data Is Inadmissible". However, the VST transformation did work as what suggested in the Vignettes. I started wondering this VST from DESeq2 may not work to my data. I emailed Prof. Holmes this morning to get her opinion. By the mean time, it would be great i can have your opinion regarding this issue. >> >> I also attached the R codes and data (R workspace file). >> >> Thank you in advance. >> >> Shucong >> >> >> >> >> Shucong Li PhD Research Associate >> ----------------------------------- >> Department of Animal Science >> University of Manitoba >> 12 Dafoe Rd, Winnipeg, Manitoba >> R3T 2N2, Canada >> >> Tel: (204) 474 9841 >> Fax: +1(204) 474 7628 >> Email: shucong.li at ad.umanitoba.ca >> ----------------------------------- >> >> >> >> dds <- DESeq(dds, fitType="local?) >> >> >> >> # ddsClean <- replaceOutliersWithTrimmedMean(dds) >> >> # ddsClean <- DESeq(ddsClean,fitType="local") >> >> >> >> plotDispEsts(dds) >> >> >> >> rld <- rlogTransformation(dds, blind=FALSE) >> >> vsd <- varianceStabilizingTransformation(dds, blind=FALSE) >> >> >> >> >> >> library("vsn") >> >> par(mfrow=c(1,3)) >> >> notAllZero <- (rowSums(counts(dds))>0) >> >> meanSdPlot(log2(counts(dds,normalized=TRUE)[notAllZero,] + 1), >> >> ylim = c(0,2.5)) >> >> meanSdPlot(assay(rld[notAllZero,]), ylim = c(0,2.5)) >> >> meanSdPlot(assay(vsd[notAllZero,]), ylim = c(0,2.5)) >> >> >> >> >> >> par(mfrow=c(1,2)) >> >> dist <- dist(t(assay(rld))) >> >> fit <- hclust(dist, method="average") >> >> plot(fit, labels=dds at colData$treat) >> >> >> >> dists <- dist(t(assay(vsd))) >> >> fit <- hclust(dists, method="average") >> >> plot(fit, labels=dds at colData$treat) >
ADD REPLYlink written 5.6 years ago by Shucong Li60
Hi Shucong, On Nov 28, 2013 11:27 AM, "Shucong Li" <lishucongnn@gmail.com> wrote: > > Hi Michael, > > I have rerun the test with 'blind' set as TRUE. The graph didn't change. It is interesting that rlogTransformation seems to work very well. > > Please allow me ask one more question about differential analysis with DESeq(). I also run the differential analysis with the same set of data. Original counts was used as suggested. I'm wondering if any variance stabilizing transformation is involved in DESeq() procedure? > There are no transformations used in DESeq(). Instead the dispersion- mean relationship is used in estimating the dispersions for each gene, which is subsequently used in the generalized linear model of raw counts. It's an important but perhaps subtle distinction between transforming data and incorporating information into modeling of raw data. Mike > Shucong > > On Nov 27, 2013, at 3:31 PM, Shucong Li <shucong.li@umanitoba.ca> wrote: > > > Hi Michael, > > > > The supplemental file 3 in the manuscript is not available yet. I was told it will be available in a week or so. I will try to rerun the test with blind set as TRUE. > > > > > > Shucong > > On Nov 27, 2013, at 3:13 PM, Michael Love <michaelisaiahlove@gmail.com> wrote: > > > >> hi Shucong, > >> > >> I'm replying and including the Bioconductor list. > >> > >> It depends on what you mean by 'not working', I see that the plots of standard deviation over mean are not flat using the two transformations from DESeq2. For both transformations, you are specifying blind=FALSE, which means that the group means are taken into consideration when estimating the dispersion-mean relationship. Then for the plot, you are calculating the standard deviation over all the samples, regardless of groups. I would say that the plot does not imply that the transformation did not work, as the dispersion-mean relationship was not calculated with respect to inter-group differences because 'blind' was set to FALSE. I can't tell if the manuscript you reference used blind dispersion estimation or not when they were using the VST (I couldn't find a link to supplemental file 3). > >> > >> hope this helps, > >> > >> Mike > >> > >> On Wed, Nov 27, 2013 at 3:21 PM, Shucong Li <shucong.li@umanitoba.ca> wrote: > >> Hello Dr. Love, > >> > >> I tried to used DESeq2 to do VST transformation with my NGS (bacterial 16srRNA based) data. I came to this approach after I read the manuscript written by Profs McMurdie and Holmes titled "Waste Not, Want Not: Why Rarefying Microbiome Data Is Inadmissible". However, the VST transformation did work as what suggested in the Vignettes. I started wondering this VST from DESeq2 may not work to my data. I emailed Prof. Holmes this morning to get her opinion. By the mean time, it would be great i can have your opinion regarding this issue. > >> > >> I also attached the R codes and data (R workspace file). > >> > >> Thank you in advance. > >> > >> Shucong > >> > >> > >> > >> > >> Shucong Li PhD Research Associate > >> ----------------------------------- > >> Department of Animal Science > >> University of Manitoba > >> 12 Dafoe Rd, Winnipeg, Manitoba > >> R3T 2N2, Canada > >> > >> Tel: (204) 474 9841 > >> Fax: +1(204) 474 7628 > >> Email: shucong.li@ad.umanitoba.ca > >> ----------------------------------- > >> > >> > >> > >> dds <- DESeq(dds, fitType="local”) > >> > >> > >> > >> # ddsClean <- replaceOutliersWithTrimmedMean(dds) > >> > >> # ddsClean <- DESeq(ddsClean,fitType="local") > >> > >> > >> > >> plotDispEsts(dds) > >> > >> > >> > >> rld <- rlogTransformation(dds, blind=FALSE) > >> > >> vsd <- varianceStabilizingTransformation(dds, blind=FALSE) > >> > >> > >> > >> > >> > >> library("vsn") > >> > >> par(mfrow=c(1,3)) > >> > >> notAllZero <- (rowSums(counts(dds))>0) > >> > >> meanSdPlot(log2(counts(dds,normalized=TRUE)[notAllZero,] + 1), > >> > >> ylim = c(0,2.5)) > >> > >> meanSdPlot(assay(rld[notAllZero,]), ylim = c(0,2.5)) > >> > >> meanSdPlot(assay(vsd[notAllZero,]), ylim = c(0,2.5)) > >> > >> > >> > >> > >> > >> par(mfrow=c(1,2)) > >> > >> dist <- dist(t(assay(rld))) > >> > >> fit <- hclust(dist, method="average") > >> > >> plot(fit, labels=dds@colData$treat) > >> > >> > >> > >> dists <- dist(t(assay(vsd))) > >> > >> fit <- hclust(dists, method="average") > >> > >> plot(fit, labels=dds@colData$treat) > > > [[alternative HTML version deleted]]
ADD REPLYlink written 5.6 years ago by Michael Love24k
Answer: Variance stabilizing transformation seems not working on NGS (bacterial 16srRNA
0
gravatar for minie
2.9 years ago by
minie0
minie0 wrote:

Hello!

 

I am facing some problem while running the same code. I am using following code:

rld <- rlogTransformation(dds, blind=TRUE)
vsd <- varianceStabilizingTransformation(dds, blind=TRUE)

library("vsn")
pdf(paste('plots/vsn_variance',expname,'.pdf',sep=''), height=4,width=12)
par(mfrow=c(1,3))
notAllZero <- (rowSums(counts(dds))>0)
meanSdPlot(log2(counts(dds,normalized=TRUE)[notAllZero,] + 1),
ylim = c(0,2.5))
meanSdPlot(assay(rld[notAllZero,]), ylim = c(0,2.5))
meanSdPlot(assay(vsd[notAllZero,]), ylim = c(0,2.5))
dev.off()
rm(rld,vsd)

and the error I am getting is          :
Error: Unknown parameters: ylim
Execution halted

Can any one help me to rectify the error.

Thanks in advance

 

ADD COMMENTlink written 2.9 years ago by minie0
Answer: Variance stabilizing transformation seems not working on NGS (bacterial 16srRNA
0
gravatar for Wolfgang Huber
2.9 years ago by
EMBL European Molecular Biology Laboratory
Wolfgang Huber13k wrote:

Minie:

Have a look at the manual page of meanSdPlot, and in particular the example. Since Bioconductor 3.2, this function no longer takes R base graphics parameters, rather it returns a ggplot object. So try

msd = meanSdPlot(...)
print(
  msd$gg + scale_y_continuous(limits = c(0, 2.5)) 
)


 

ADD COMMENTlink modified 2.9 years ago • written 2.9 years ago by Wolfgang Huber13k

Dear Wolfganag thanks for your reply, but with the limits also I am facing same problem :

Error: Unknown parameters: limits
Execution halted

ADD REPLYlink written 2.9 years ago by minie0

Try instead

+ ylim(0,2.5)

 

ADD REPLYlink written 2.9 years ago by Michael Love24k
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 16.09
Traffic: 134 users visited in the last hour