Question

mean log2CPM after model fit (lmFit)

0

Entering edit mode

anna • 0

@anna-23212

Last seen 4.7 years ago

Hi there, I am completely stumped in the analysis of my RNAseq data set using limma. But is there is a possibility to extract the mean expression levels for each sample (averaged among 3 biol. replicates and adjusted for block effects) from the lmFit output or can these be calculate with the help of the estimated coefficients? I have 4 genotypes and for each genotype 3 biological replicates. Prior to do any differential expression analysis, I would like extract the mean expression levels for each gene and each genotype adjusted for any block effects in my experimental design.

so far, I followed the limma pipeline:

design <- model.matrix(~0 + genotype + block) 
v <- voom(dge, design, plot=TRUE)
fit <- lmFit (v, design)

Thanks for your input and help.

limma lmFit mean expression level • 916 views

ADD COMMENT • link updated 4.7 years ago by James W. MacDonald 67k • written 4.7 years ago by anna • 0

score 2 · Accepted Answer · 2020-03-30

2

Entering edit mode

James W. MacDonald 67k

@james-w-macdonald-5106

Last seen 1 day ago

United States

You are asking about the coefficients, which are in

fit$coef

ADD COMMENT • link 4.7 years ago James W. MacDonald 67k

0

Entering edit mode

Thanks for your reply. So the fit$coef of lmFit are the adjusted mean log2CPM for each genotype? I guess I got confused with the coeficients used for DE determination and the lmFit function description, which says: "The coefficients of the fitted models describe the differences between the RNA sources hybridized to the arrays".

ADD REPLY • link 4.7 years ago anna • 0

1

Entering edit mode

What the coefficients describe is a function of the design matrix that you used. Your design matrix is a 'cell means' model, which computes the mean for each group. Since you also have block in the model, the computed means are adjusted for the block factor.

ADD REPLY • link 4.7 years ago James W. MacDonald 67k

1

Entering edit mode

You also want to specify

contrasts(block) <- contr.sum(levels(block))

to make sure that the block effects add to zero. Then the block effects won't bias the genotype coefficient, which will then be interpretable as genotype means.

ADD REPLY • link 4.7 years ago Gordon Smyth 51k

0

Entering edit mode

Thanks, James and Gordon, for your help and clarifications! My R script to obtain genotype means adjusted for block effects:

contrasts(block) <- contr.sum(levels(block))
design <- model.matrix(~0 + genotype + block) 
v <- voom(dge, design, plot=TRUE)
fit <- lmFit (v, design)
means <- fit$coefficients

ADD REPLY • link 4.7 years ago anna • 0

0

Entering edit mode

You actually need

means <- fit$coefficients[1:4]

ADD REPLY • link 4.7 years ago Gordon Smyth 51k