Question: Limma residuals matrix
0
gravatar for agustin.gonvi
12 weeks ago by
Cleveland, OH
agustin.gonvi10 wrote:

I am looking for sexual dimorphism in gene expression so I fit a model

design <- model.matrix(~Sex, data = pData(eset))
fit <- lmFit(eset, design)

I need to obtain a matrix with the residuals of the linear model.

I look at the object fit and I cannot interpret where the residuals are; according to the manual fit$df.residual contains the degrees of freedom which is not what I need.

Is it possible to obtain the residuals matrix?

Thanks!

limma • 117 views
ADD COMMENTlink modified 12 weeks ago by Gordon Smyth39k • written 12 weeks ago by agustin.gonvi10
Answer: Limma residuals matrix
3
gravatar for Gordon Smyth
12 weeks ago by
Gordon Smyth39k
Walter and Eliza Hall Institute of Medical Research, Melbourne, Australia
Gordon Smyth39k wrote:
ResidualsMatrix <- residuals(fit, eset)
ADD COMMENTlink modified 12 weeks ago • written 12 weeks ago by Gordon Smyth39k
Answer: Limma residuals matrix
1
gravatar for Steve Lianoglou
12 weeks ago by
Denali
Steve Lianoglou12k wrote:

If you're looking for the genes that are dimorphic based on sex, you've just found them. If you follow through with your "standard" limma dge pipeline, you'll get the genes differentially expressed across the sexes of your samples, ie.

design <- model.matrix(~Sex, data = pData(eset))
fit <- lmFit(eset, design)
fit <- eBayes(fit)
res <- topTable(fit, n = Inf)

The p-values and logFCs reported in the res data.frame are those for male vs female, the order of which is determined by the levels(eset$Sex). Just look at the 2nd column name of your design matrix. Positive logFC's will be those that are more highly expressed in that group vs the intercept (the other sex).

Still, if you really want the residuals matrix, you could do:

resids <- limma::removeBatchEffect(eset, batch = eset$Sex)

Caveat emptor: this should work "just fine" assuming the data in exprs(eset) is appropriate (ie. on log2 scale and more-or-less normally distributed)

ADD COMMENTlink written 12 weeks ago by Steve Lianoglou12k

Thanks! I did find the sex-DEG. And also need the residuals. I wasn't sure whether using removeBatchEffect to regress sex out, would use the same linear model as lmFit.

I was wondering if there is an equivalent of resid() applied on the lm() output.

ADD REPLYlink written 12 weeks ago by agustin.gonvi10
2

Well, I'll be.

It does look like there is a limma::residuals.MArrayLM function defined. Try residuals(fit, eset) and see how it compares to resids above ...

ADD REPLYlink written 12 weeks ago by Steve Lianoglou12k
1

Yes that's right. Using removeBatchEffect is not the same as it does not subtract out the intercept.

ADD REPLYlink modified 12 weeks ago • written 12 weeks ago by Gordon Smyth39k
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: 333 users visited in the last hour