Hi,

I am following section 9.5.2 - 9.5.4 of the limma user guide. I understand that each of the output of each of these three sections should yield identical results, but my results differ in p-value and adjusted p-value, so I suspect that there are issues with my code but can't seem to find out what it is. This is a pre/post design where every patient has received a treatment. Certain patients are responders, others are not. We're interested in this specific interaction effect: the effect of treatment in responders versus non-responders

When I follow section 9.5.2:

make sure treatment is a factor variables with proper reference

>phenotype$treatment <- as.factor(phenotype$treatment)
>phenotype$treatment <- relevel(phenotype$treatment, ref = "0")

make sure responder is a factor variables with proper reference

>phenotype$responder <- as.factor(phenotype$responder)
>phenotype$responder <- relevel(phenotype$responder, ref = "non")

convert to cpm --> between sample, controlling for sequencing depth

>countsPerMillion <- cpm(dgList) ##divide each count per million for comparison

we retain only those genes that are represented at least 1cpm reads in at least two samples (cpm=counts per million)

>summary(countsPerMillion) 
>countCheck <- countsPerMillion > 1 
>head(countCheck)
>keep <- which(rowSums(countCheck) >= 2)
>dgList <- dgList[keep,]
>summary(cpm(dgList))

TMM normalization controls for library size

>dge <- calcNormFactors(dgList, method="TMM")

design matrix

>phenotype$treat <- factor(paste(phenotype$treatment,phenotype$responder,sep="."))
>treat <- paste(phenotype$responder,phenotype$treatment,sep=".")
>treat <- factor(treat, levels=c("non.0","non.1","responder.0","responder.1"))
>design <-model.matrix(~0 + treat + cov1 + cov2 + cov3 + cov4 + cov5 + cov6 + cov7,data=phenotype)
>colnames(design)[1:4] <- levels(treat)

estimate within subject correlation

>counts <- voom(dge,design)
>corfit <- duplicateCorrelation(counts,design,block=phenotype$study_id)
>corfit$consensus
>counts_voom <- voom(dge,design,block=phenotype$study_id,correlation=corfit$consensus)

account for correlation in a model with random effect

>fit <- lmFit(counts_voom,design,block=phenotype$study_id,correlation=corfit$consensus)
>fit <- eBayes(fit)

make contrasts

>cm <- makeContrasts(
  DOD = (responder.1 - responder.0) - (non.1 - non.0),
  RTE = responder.1 - responder.0,
  NRTE = non.1 - non.0,
  levels=design
)

fit model

>fit2 <- contrasts.fit(fit,cm)
>fit2 <- eBayes(fit2)
>results <- topTable(fit2,coef="DOD")

when I follow section 9.5.3:

create design matrix

>design <-model.matrix(~responder + responder:treatment + cov1 + cov2 + cov3 + cov4 + cov5 + cov6 + cov7,data=phenotype)

estimate within subject correlation

>counts <- voom(dge,design)
>corfit <- duplicateCorrelation(counts,design,block=phenotype$study_id)
>corfit$consensus
>counts_voom <- voom(dge,design,block=phenotype$study_id,correlation=corfit$consensus)

account for correlation in a model with random effect

>fit <- lmFit(counts_voom,design,block=phenotype$study_id,correlation=corfit$consensus)
>fit <- eBayes(fit)

make contrasts corresponding to the respondernon:treatment1 and responderresponder:treatment1 to extract interaction of the effect of treatment in responders versus nonresponders

>fit2 <- contrasts.fit(fit,c(0,0,0,0,0,0,0,0,0,-1,1)
>fit2 <- eBayes(fit2)
>results <- topTable(fit2)

The output from colnames(fit$design) is:

> [1] "Intercept"  "responderresponder"
> [3]  "cov1" "cov2"
> [5]  "cov3" "cov4"
> [7]  "cov5" "cov6"
> [9]  "cov7" "respondernon:treatment1"
> [11] "responderresponder:treatment1"

Based on this, and info gathered from my basic LDA class, I interpreted the coefficients as such: Intercept - The treatment0 and cov1 thru cov7 at 0 effect in nonresponder responderresponder - the difference between the non-responders and responders at treatment0 and cov1 thru cov7 at 0 respondernon:treatment1 - the difference between nonresponders at treatment 0 and treatment 1 adjusted for cov1-cov7 responderresponder:treatment1 - the difference between responders at treatment 0 and treatment 1 adjusted for cov1-cov7

In this case, I would assume that the difference between responderresponder:treatment1 and respondernon:treatment1 would be the interaction of interest. However, I don't believe that this is the case, given that the results I run with the 9.5.2 results don't match these (a la 9.5.3) results.