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

