Variance stabilizing transformation seems not working on NGS (bacterial 16srRNA based) data
4
0
Entering edit mode
@mikelove
Last seen 10 hours ago
United States
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 Microbiome DESeq2 • 3.6k views
ADD COMMENT •
0
Entering edit mode
Shucong Li ▴ 60
@shucong-li-6266
Last seen 8.8 years ago
Canada
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 COMMENT •
0
Entering edit mode
Shucong Li ▴ 10
@shucong-li-6267
Last seen 10.3 years ago
An embedded and charset-unspecified text was scrubbed... Name: not available URL: <https: stat.ethz.ch="" pipermail="" bioconductor="" attachments="" 20131127="" c2602b76="" attachment.asc="">
ADD COMMENT •
0
Entering edit mode
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 REPLY •
0
Entering edit mode
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 REPLY •
0
Entering edit mode
minie • 0
@minie-11306
Last seen 7.1 years ago

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 COMMENT •
0
Entering edit mode
@wolfgang-huber-3550
Last seen 3 months ago
EMBL European Molecular Biology Laborat…

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 COMMENT •
0
Entering edit mode

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

Error: Unknown parameters: limits
Execution halted

ADD REPLY •
0
Entering edit mode

Try instead

+ ylim(0,2.5)

 

ADD REPLY •

Login before adding your answer.

Similar Posts
Loading Similar Posts
Traffic: 659 users visited in the last hour
Content Search
Users
Tags
Badges
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6