Complex contrasts in limma
1
0
Entering edit mode
tcalvo ▴ 100
@tcalvo-12466
Last seen 18 months 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.

About the target variables:

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 • 1.5k views
ADD COMMENT
2
Entering edit mode
@gordon-smyth
Last seen 9 hours ago
WEHI, Melbourne, Australia

Your interpretations seem to match your contrast definitions.

ADD COMMENT
0
Entering edit mode

Thank you very much.

ADD REPLY

Login before adding your answer.

Traffic: 594 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