Variance stabilizing transformation seems not working on NGS (bacterial 16srRNA based) data

0

Entering edit mode

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]]

0

Entering edit mode

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]]

0

Entering edit mode

Shucong Li
▴
10

@shucong-li-6267
Last seen 7.0 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="">

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
• link
7.8 years ago
Shucong Li
▴
60

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
• link
7.8 years ago
Michael Love
34k

0

Entering edit mode

minie
•
0

@minie-11306
Last seen 3.8 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

0

Entering edit mode

Wolfgang Huber
13k

@wolfgang-huber-3550
Last seen 7 weeks 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)) )

Similar Posts

Loading Similar Posts

Traffic: 340 users visited in the last hour

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