Question: DESeq2: How does the model of the trend function looks like?
0
gravatar for Fischer-philipp
6 months ago by
Fischer-philipp20 wrote:

Hello Community,

The paper of DESe1 says that the forumula (6) :
A parametric curve of the form (6) is fit by regressing the gene-wise dispersion estimates α_g^wi onto the means of the normalized counts, μ̄ i via a gamma-family GLM regression.

Wihtin the Source code of DESeq2 I found the following code:

# Estimate a parametric fit of dispersion to the mean intensity
parametricDispersionFit <- function( means, disps ) {
   coefs <- c( .1, 1 )
   iter <- 0
   while(TRUE) {
      residuals <- disps / ( coefs[1] + coefs[2] / means )

So am I right that one assumes something like that? Atr = a1 / m + a0 + error

Atr /( a1 / m + a0 ) = error /( a1 / m + a0 ) mit epsilon = error /( a1 / m + a0 )

Atr /( a1 / m + a0 ) = epsilon ~ Gamma(a0, a1)

Thank you in advance

ADD COMMENTlink modified 6 months ago by Michael Love25k • written 6 months ago by Fischer-philipp20
Answer: DESeq2: How does the model of the trend function looks like?
1
gravatar for Michael Love
6 months ago by
Michael Love25k
United States
Michael Love25k wrote:

By the way, if you want to see how it looks, you can run plotDispEsts, which plots the dispersion estimates.

https://bioconductor.org/packages/release/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#dispersion-plot-and-fitting-alternatives

The trend function is:

alpha_tr(x) = a1 / x + a0

No error, that's the function for the trend line.

Then the MLE dispersion estimates are assumed to follow a Gamma distribution around the trend, for the purposes of estimation:

alphai^MLE ~ Gamma(mean=alphatr(mu-bar_i), theta)

We use glm to perform the estimation of the coefficients a1 and a0.

ADD COMMENTlink written 6 months ago by Michael Love25k

Thank you very much for your very helpfull answers. Unfortunately I am stuck in the section of Dispersion Prior:

Additional file 1: Table S2 compares: - ψ1((m−p)/2) as an approximation of σ2lde with - the variance of logarithmic Cox–Reid adjusted dispersion. And says it is very similar.

As a result:

Therefore, the prior variance σ2d is obtained by subtracting the expected sampling variance from an estimate of the variance of the logarithmic residuals, s2lr:

σ2d=max{slr2−ψ1((m−p)/2),0.25}.

I do not really get why one substracts the squared logarithmic residuals from the trigamma approximated logarithmic dispersion estimators due to the fact that the logarithmic Cox–Reid adjusted dispersion are similar to the trigamma approximated ones. in order to obtain the prior variance of log(a_i).

I think I am missing something conceptionally - it would be great if you could give me a hint.

Thank you in advance.

ADD REPLYlink written 6 months ago by Fischer-philipp20
1

The subtraction is from what you would expect with a hierarchical Normal. The variances of the Normals add, so to get a reasonable prior from the observed distribution, we subtract the expected sampling variance.

ADD REPLYlink written 6 months ago by Michael Love25k

So I am wonderin what the hierarchical normal model looks like:

  • ) log(ai^gw) - log(atr(\bar mui) ~ N(aigw, Slr)

  • ) ai^gw ~ N(ai; sigma^i_lde)

  • ) log(ai) ~ N(log(atr(\bar mui); Sigma^2d) as prior

It is really interesting for me what you did - but I am really stuck.

Thank you again.

ADD REPLYlink modified 6 months ago • written 6 months ago by Fischer-philipp20
1

Yes, the second two lines are the assumptions that motivate subtracting the two estimates.

ADD REPLYlink written 6 months ago by Michael Love25k

Hey - Thank you again.

So I was trying to draw sthg like this: Hierarchical normal https://pasteboard.co/IdKFNpx.png Which is from: http://idiom.ucsd.edu/~rlevy/pmsltextbook/chapters/pmsl8.pdf

hierarchical normal model deseq2? https://pasteboard.co/IdKGGQF.jpg

Does it make sense?

Thank you in advance

ADD REPLYlink modified 5 months ago • written 5 months ago by Fischer-philipp20

Yes, that kind of diagram makes sense.

ADD REPLYlink written 5 months ago by Michael Love25k

=) nice thanks.

Further it would be great if you are so kind and comment my suggestion (the link) of how to derive the sigma^2d of the normal hierarchical.

model of hierarchical normal with full bayes https://ibb.co/xHcZ2MR

ADD REPLYlink written 5 months ago by Fischer-philipp20
1

This is beyond the scope of the support I can provide here. On the support site, I try to answer basic questions about how to use the software, or clarification questions like the above thread.

ADD REPLYlink written 5 months ago by Michael Love25k

I understand that - thank your for your help.

ADD REPLYlink written 5 months ago by Fischer-philipp20

Another question came to my mind.
Am I getting the cox-reid adjustment right. 1/2 log det xt w x One penalized values for alpha which have a lot of information on mu?

If so I do not get why this makes sense.
To emphasize on values of alpha coming more from the poisson distribution?

Thanks again

ADD REPLYlink modified 5 months ago • written 5 months ago by Fischer-philipp20

hi Philipp,

Again, I'm kind of pressed for time these days and need to reserve my support site allotment for pressing questions about usage of the software. I try to make myself very available for these questions, but I received for example 32 questions since Monday, so I have to be as efficient as possible. You can look over the Methods, and the open source code, and follow the references in the original paper to understand our use of the Cox-Reid adjustment, where we follow edgeR. I don't have time to check your notes here.

ADD REPLYlink modified 5 months ago • written 5 months ago by Michael Love25k

Hey Michael, Yes your support is great! I am sorry that I bothered you with my questions. All the best =).

ADD REPLYlink written 5 months ago by Fischer-philipp20
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: 471 users visited in the last hour