Complex contrasts in limma
1
0
Entering edit mode
tcalvo ▴ 70
@tcalvo-12466
Last seen 12 weeks ago
Brazil

Dear colleagues,

I'm analyzing a dataset with both within- and between-individuals comparisons. I'm insecure about some of contrasts and wondered if someone more experienced can take a look.

disease = between-subject (disease, control)
treatment = within-subject (after, before)
clin = between-subject (symptomatic, asymptomatic)

To reproduce my pdata.

Below there's an excerpt of the code.

Any input is appreciated. Thank you!


y <- DGEList(counts = txi$counts, samples = pdata) # collapse 3 factor vars. into one pdata$groups <- with(pdata, paste(disease, severity, treat, sep = "_"))

y <- y[filterByExpr(y, design = model.matrix(~groups, metadata)), ]
dge1 <- calcNormFactors(y)

# model matrix without intecept
design1 <- model.matrix(~0+groups+sex+age+bmi, dge1$samples) v <- voomWithQualityWeights(dge1, design1, plot = TRUE) # random effects for patient, 1st pass corfit <- duplicateCorrelation(v, design1, block = dge1$samples$pat_id) corfit$con

v <- voomWithQualityWeights(dge1, design1, block = dge1$samples$pat_id,
correlation = corfit$consensus) corfit <- duplicateCorrelation(v, design1, block = dge1$samples$pat_id) corfit$con

fit1 <- lmFit(v, design1, block = dge1$samples$pat_id,
correlation = corfit$consensus.correlation) # Contrasts cont1 <- makeContrasts( Disease_Symptom = groupsDISEASE_SYMPTOMATIC_AFTER - # Disease after vs before treat (symptomatic) groupsDISEASE_SYMPTOMATIC_BEFORE, Control_Asymptom = groupsCONTROL_ASYMPTOMATIC_AFTER - # Control after vs. before treat (asymptom) groupsCONTROL_ASYMPTOMATIC_BEFORE, Disease_Asymptom = groupsDISEASE_ASYMPTOMATIC_AFTER - # Disease after vs. before treat (asympto) groupsDISEASE_ASYMPTOMATIC_BEFORE, Treat_Diff_Disease_Symptom = # Difference in treatment response between disease symptomatic vs asymptomatic (diff of diff) (groupsDISEASE_SYMPTOMATIC_AFTER - groupsDISEASE_SYMPTOMATIC_BEFORE) - (groupsDISEASE_ASYMPTOMATIC_AFTER - groupsDISEASE_ASYMPTOMATIC_BEFORE), Disease_Symp_Control = ( # Difference in treatment response between disease symptom vs control asymptom (diff of diff) groupsDISEASE_SYMPTOMATIC_AFTER - groupsDISEASE_SYMPTOMATIC_BEFORE) - (groupsCONTROL_ASYMPTOMATIC_AFTER - groupsCONTROL_ASYMPTOMATIC_BEFORE), Eff_Treat = ( # Average efect of treatment regardless of symptom and disease groupsDISEASE_SYMPTOMATIC_AFTER + groupsCONTROL_ASYMPTOMATIC_AFTER)/2 - (groupsDISEASE_SYMPTOMATIC_BEFORE + groupsCONTROL_ASYMPTOMATIC_BEFORE)/2, levels = design1) fit2 <- contrasts.fit(fit1, cont1) fit2 <- eBayes(fit2, robust = TRUE) colnames(fit2$coefficients)

knitr::kable(summary(decideTests(fit2, adjust.method = "BH", p.value = 0.1)), format = "simple")

#           Disease_Symptom   Control_Asymptom   Disease_Asymptom   Treat_Diff_Disease_Symptom   Disease_Symp_Control   Eff_Treat
# -------  ----------------  -----------------  -----------------  ---------------------------  ---------------------  ----------
# Down                   27                151                  0                          178                      0         952
# NotSig              14134              13753              14417                        13980                  14417       11983
# Up                    256                513                  0                          259                      0        1482


The coefficients seem to agree with graphical EDA and descriptive summaries.

RNASeq limma DifferentialExpression • 207 views
2
Entering edit mode
@gordon-smyth
Last seen 4 hours ago
WEHI, Melbourne, Australia

0
Entering edit mode

Thank you very much.