**380**wrote:

Hi

So I am trying to work out why these two bits of code give different results. The overall aim is to look for genes that are changing differently between responders and non-responders (to therapy) between two different treatment groups (baseline only data i.e. a single timepoint). However, I am also interested if there are genes that are changing differently between the two treatment groups aside from responder and non-responder, kind of as a control to make sure the two groups are similar.

This code is (I believe) for answering the last question, which is to look for genes that are changing depending on medication assignment:

### code 1

design <- model.matrix(~0+meta$Current.Medication+meta$Age+meta$Gender+meta$Ethnicity)

colnames(design) <- c('Drug1','Drug2','Age','Gender','Ethnicity*A','Ethnicity*C','Ethnicity_O')

contrast.matrix <- makeContrasts('Drug1-Drug2', levels = design)

dge <- DGEList(counts=countsn)

keep <- filterByExpr(dge, design)

dge <- dge[keep,,keep.lib.sizes=FALSE]

dge <- calcNormFactors(dge)

v <- voom(dge, design, plot=TRUE)

fit <- lmFit(v,design)

fit2 <- contrasts.fit(fit, contrast.matrix)

fit2 <- eBayes(fit2)

top2 <- topTable(fit2,coef=1,number=Inf,sort.by="P")

qs <- qvalue::qvalue(p = top2$P.Value)

top2$q.Val <- qs$qvalues

nrow(subset(top2, top2$q.Val < 0.05))

[1] 0

So 0 DEGs are found.

This next code can test the interaction of medication*response, but also can (I think) compare these two medication groups Drug1-Drug2 like I have done above. However, the results for the medication only comparison below are different.

### code 2

design <- model.matrix(~meta$Current.Medication*meta$DAS28.CRP.EULARresp.V7+meta$Age+meta$Gender+meta$Ethnicity)

colnames(design) <- c('Intercept','Medication','Response','Response:Medication','Age','Gender','Ethnicity*A','Ethnicity*C','Ethnicity_O')

dge <- DGEList(counts=countsn)

keep <- filterByExpr(dge, design)

dge <- dge[keep,,keep.lib.sizes=FALSE]

dge <- calcNormFactors(dge)

v <- voom(dge, design, plot=TRUE)

fit <- lmFit(v,design)

cont.matrix <- cbind(Medication=c(0,1,0,0,0,0,0,0,0),Response=c(0,0,1,0,0,0,0,0,0),Diff=c(0,0,0,1,0,0,0,0,0))

fit2 <- contrasts.fit(fit, cont.matrix)

fit2 <- eBayes(fit2)

top2 <- topTable(fit2,coef=1,number=Inf,sort.by="P")

qs <- qvalue::qvalue(p = top2$P.Value)

top2$q.Val <- qs$qvalues

nrow(subset(top2, top2$q.Val < 0.05))

[1] 14451

So 14451 DEGs instead of 0, seems like a different test has been done, but I thought these would be the same.

Why do I see such dramatic differences? Can some one help?

Thanks,

Chris

**38k**• written 4 months ago by chris86 •

**380**