Hi there, I'm new to working with EdgeR's generalized linear models, but I have a data set where I'm trying to get complex comparisons between groups, times, and treatments, so here I am...

I've figured out that the best approach for me is to use the cell means model. Right now I have two groups and two treatments (at each time point, but ignore that for now).

Here is my code:

#countsDGE_d3 is my count matrix
>key_3d
                       sample    group condition
F.3dayInjury.IP1 F.3dayInjury forelimb    injury
F.3dayInjury.IP2 F.3dayInjury forelimb    injury 
F.3dayInjury.IP3 F.3dayInjury forelimb    injury
F.Intact.IP1         F.Intact forelimb no_lesion 
F.Intact.IP2         F.Intact forelimb no_lesion 
F.Intact.IP3         F.Intact forelimb no_lesion 
H.3dayInjury.IP1 H.3dayInjury hindlimb    injury 
H.3dayInjury.IP2 H.3dayInjury hindlimb    injury 
H.3dayInjury.IP3 H.3dayInjury hindlimb    injury 
H.Intact.IP1         H.Intact hindlimb no_lesion 
H.Intact.IP2         H.Intact hindlimb no_lesion 
H.Intact.IP3         H.Intact hindlimb no_lesion 
H.Intact.IP4         H.Intact hindlimb no_lesion

design3d <- model.matrix(~0+key_3d$sample)
colnames(design3d) <- levels(key_3d$sample)

countsDGE_d3 <- calcNormFactors(countsDGE_d3)
countsDGE_d3 <- estimateDisp(countsDGE_d3, design3d, robust = TRUE)
countsDGE_d3$common.dispersion

fit <- glmQLFit(countsDGE_d3, design3d)

cont.matrix <- makeContrasts(NvsINJinFL=F.3dayInjury-F.Intact, 
                             NvsINJinHL=H.3dayInjury-H.Intact,
                             Diff=(F.3dayInjury-F.Intact)-(H.3dayInjury-H.Intact),
                             levels=design3d)

qlf <- glmQLFTest(fit, contrast=cont.matrix[,"NvsINJinFL"])
topTags(qlf, n = 20)
table <- topTags(qlf, n = Inf)
write.csv( as.data.frame(table), file="All_FLDay3_InjuredvsIntact_results.csv")

I then repeat running glmQLFTest on all three contrasts and export the topTags table.

My question is, due to the data being pretty noisy, I used the RUVseq package to calculate the factors of unwanted variation for my samples. I would like to incorporate these into my design matrix and into my contrasts. How do I construct this model matrix and contrast matrix?

Does this make sense?

design3d <- model.matrix(~0+key_3d$sample+W_1, data=pData_3d)

I'm not sure what my contrasts should be. I tried:

cont.matrixW1 <- makeContrasts(NvsINJinFL=(F.3dayInjury-F.Intact)-W_1, 
                               NvsINJinHL=(H.3dayInjury-H.Intact)-W_1,
                               Diff=((F.3dayInjury-F.Intact)-(H.3dayInjury-H.Intact))-W_1,
                               levels=design3dW1)

Where pData_3d is a matrix of the factors of unwanted variation from:

seq_set <- newSeqExpressionSet(as.matrix(dge_counts, phenoData = key$sample))

design <- model.matrix(~0 + key$sample)
colnames(design) <- levels(key$sample)
dge <- DGEList(counts = dge_counts, group = key$sample)
dge <- calcNormFactors(dge, method="upperquartile", rubust = TRUE)
dge <- estimateGLMCommonDisp(dge, design)
dge <- estimateGLMTagwiseDisp(dge, design)

fit <- glmFit(dge, design)
res <- residuals(fit, type = "deviance")

#use all the genes to estimate the factors of unwanted variation.
seqUQ <- betweenLaneNormalization(seq_set, which="upper")
controls <- rownames(seq_set)
seqRUVr <- RUVr(seqUQ, controls, k=1, res)

Thank you!