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

# 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))

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.