linear mixed models using edgeR+limma
1
1
Entering edit mode
weichengz ▴ 10
@weichengz-23557
Last seen 3.4 years ago
Melbourne, Australia

Hi there,

I am analysing my sequencing data with edgeR and limma packages. I was looking for the DEGs in fetal muscle as a result of a maternal stress treatment during gestation in pigs (HS vs. Control). The individual mother (n=15) was the biological replicate as they received the treatment (and control) while multiple fetuses (n=3 per litter) from each mother/litter were used for RNA-seq.

I hope to use a linear mixed model and to include the main effects (fixed effects) of treatment, fetal sex and their interactions with 'Litter' being the random effect. Also, I used voomQualityWeights as I've got outliers.

I've known that limma package allows us to use linear mixed models, but I am not sure if my code below was correct?

Your correction and suggestions are very much appreciated!

Cheers, Weicheng

sample matrix

# RNA-seq data fetal muscle 

#import table
countall<- read.table(file.choose(),header = F, row.names = 1)
countall

#convert to matrix
x <- as.matrix(countall)

#Create a DGEList object
Treatment <- factor(c(rep("Control",12),rep("HS",12),rep("Control",9),rep("HS",12)))
y <- DGEList(counts = x, group = Treatment)

#include Sex and Litter in the factor model

Sex <- factor(c("male","male","male","female","male","female","male","female","male","female","male","male","female","male","male","female","female","female","female","female","male","female","female","female","female","female","female","female","female","female","female","male","female","female","male","female","female","female","male","female","male","female", "female","male","female"))

Litter < factor(c("1","1","1","2","2","2","3","3","3","4","4","4","5","5","5","6","6","6","7","7","7","8","8","8","9","9","9","10","10","10","11","11","11","12","12","12","13","13","13","14","14","14", "15","15","15"))

y$samples$Sex <- Sex
y$samples$Litter <- Litter
y$samples

#Filter and convert to logCPM

keep <- filterByExpr(y)
summary(keep)
y <- y[keep, , keep.lib.sizes=FALSE]

#Apply TMM (trimmed mean of M-values) normalization 
y <- calcNormFactors(y,method = "TMM")

#Transformations from the raw-scale
logCPM <- cpm(y, log=TRUE, prior.count=2)

#Random effect correlation with Litter as a random effect

design <- model.matrix(~0+Treatment+Sex,data=y$samples)
dupC <- duplicateCorrelation(logCPM, design, block=y$samples$Litter)
dupC$consensus

# Analysis with combined voom and sample quality weights

vwts <- voomWithQualityWeights(logCPM, design, plot=TRUE)

#Then this inter-Litter correlation is input into the linear model fit:

vfit2 <- lmFit(vwts,design,block=y$samples$Litter,correlation=dupC$consensus)
efit3 <- eBayes(vfit2)
summary(decideTests(efit3))
linear mixed model edger limma voom RNA-seq • 2.5k views
ADD COMMENT
3
Entering edit mode
@gordon-smyth
Last seen 3 hours ago
WEHI, Melbourne, Australia

Your code is a bit inconsistent because you are mixing code from a few different limma pipelines. Assuming you're running R 4.0.0 and you want to use voom, the easiest and best would be:

y <- calcNormFactors(y)
design <- model.matrix(~Treatment+Sex)
fit <- voomLmFit(y, design, block=Litter, plot=TRUE)
fit <- eBayes(fit)
topTable(fit, coef="TreatmentHS")

You will need R 4.0.0 for this code because voomLmFit is new in the latest release of edgeR. The new function combines voom, duplicateCorrelation and lmFit into one, and also gives better protection against data with lots of zeros than the original functions would if run separately.

ADD COMMENT
0
Entering edit mode

Thanks so much Gordon.

I meant to follow the paper RNA-seq analysis is easy as 1-2-3 with limma,Glimma and edgeR

I used R3.6.1, does voomLmFit function also combines the voomQualityWeights which deals with outliers ?

Cheers,Weicheng

ADD REPLY
0
Entering edit mode

Yes it does, just add sample.weights=TRUE to the call:

fit <- voomLmFit(y, design, block=Litter, sample.weights=TRUE, plot=TRUE)
ADD REPLY
0
Entering edit mode

Hi Gordon,

Can I please further ask how to call interactions if I use design <- model.matrix(~Treatment+Sex)? For example, HS vs. Control in females (or males). And Female vs. Male in HS group?

Cheers, W

ADD REPLY
0
Entering edit mode

The model you are using assumes no interactions. If you want to make these extra comparisons, you would need to fit a design with four groups (Control.Female, Control.Male, HS.Female and HS.Male).

ADD REPLY
0
Entering edit mode

Hi Gordan,

Just to follow on from this post. I fit a design with four groups using the following codes. Would this be correct for each comparison of interest? What's the difference between design2 <- model.matrix(~0+combined) and design <- model.matrix(~Sex * Treatment) here?

If I would like to know the overall treatment effect (HS vs. Control) across Sex, should I go back and use design <- model.matrix(~Treatment+Sex) and call topTable(fit, coef="TreatmentHS")? Thank you.

y <- calcNormFactors(y)

combined <- paste(y$samples$group,y$samples$Sex,sep = ".")

combined <- factor(combined,levels=c("Control.male","Control.female","HS.male","HS.female"))

combined

design2 <- model.matrix(~0+combined)
colnames(design2) <- levels(combined)
colnames(design2)

fit2 <- voomLmFit(y, design2, block=Litter, sample.weights = TRUE, plot=TRUE)

my.contrast <- makeContrasts(HSvsControlinFemale=HS.female-Control.female,
                             HSvsControlinMale=HS.male-Control.male,
                             MalevsFemaleinHS=HS.male-HS.female,
                             MalevsFemaleinControl=Control.male-Control.female,levels = design2)
fit3 <- contrasts.fit(fit2, my.contrast)
fit3 <- eBayes(fit3)
summary(decideTests(fit3))

#HSvsConntrol in female
topTable(fit3, coef = "HSvsControlinFemale", sort.by="P",n=Inf)

#HSvsControl in male
topTable(fit3, coef = "HSvsControlinMale",sort.by="P",n=Inf)

#MalevsFemale in Control
topTable(fit3, coef = "MalevsFemaleinControl",sort.by="P",n=Inf)

#MalevsFemale in HS
topTable(fit3, coef = "MalevsFemaleinHS",sort.by="P",n=Inf)
ADD REPLY

Login before adding your answer.

Traffic: 837 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

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

Powered by the version 2.3.6