Hi,

I'm analysing RNA-seq data and have an experimental set up similar to this:

colData <- data.frame(sample=c(1,2,3,4,5,6,7,8),
                      group=c("A", "B", "C", "A", "B", "A", "C", "C"),
                      patient=c(1,1,1,2,2,3,3,4),
                      age=c(13.2, 15, 16.8, 13.4, 15.3, 12.8, 17, 16.5))

I'd like to do all pairwise group comparisons (AvB, AvC, BvC), whilst taking into account the ages of the patients and matched patient samples. However the issue is not all my patient samples are matched and in all groups.

I've decided to use limma, currently I understand duplicateCorrelation() can be used to account for missing matched patient samples?

I've generated the code below following the limma manual and other similar questions I've found online:

# reads is an expression matrix with columns as samples and rows as genes
dge <- DGEList(counts=reads)
# remove none expressed genes
cutoff <- 1
drop <- which(apply(cpm(dge), 1, max) < cutoff)
d <- dge[-drop,]
d <- calcNormFactors(d)

# fit model 
design = model.matrix(~0+colData$group+colData$age)
colnames(design) <- c("A", "B", "C", "age")

# normalise reads using voom
v <- voom(d, design=design)
# here the patient is treated as a random effect due to performing between and within patient comparisons
corfit <- duplicateCorrelation(v, design, block=colData$patient)

# repeat voom and duplicate correlation (https://support.bioconductor.org/p/59700/)
v <- voom(d, design=design, block=colData$patient, correlation=corfit$consensus)
corfit <- duplicateCorrelation(v, design, block=colData$patient)

fit <- lmFit(v, design=design, block=colData$patient, correlation=corfit$consensus)

contr <- makeContrasts(A - B, levels=design)

fit2 <- contrasts.fit(fit, contr)
fit2 <- eBayes(fit2)
top.table <- topTable(fit2, sort.by = "p", n = Inf)
head(top.table, 20)

The code runs fine, but I'm doubting my model is correct as the lowest adjusted p-value is 0.97 for any of the comparisons so no genes are coming up as significantly differentially expressed.

This may just be the data, but want to double check there is no issue with my model/code as I'm new to using limma. Any help appreciated, thanks!

The code seems fine but the analysis you are conducting has almost no statistical power to detect DE because the three patient groups are strongly correlated with age and you are adjusting for age. Groups A, B and C respectively consist of patients of increasing age, and therefore any gene that is DE between A and C and intermediate in B will be impossible to detect after adjusting for age effects.

I do not understand how you can have three different ages for the same patient, but that is another mattter.