Design matrix after RUVg correction
1
0
Entering edit mode
@14ef1b09
Last seen 3 days ago
Egypt

I have 3 estimated W (factors of unwanted variation) after RUVg correction. I am trying to include W (estimated factors of unwanted variation) into the design matrix just as recommended by Prof. Davide Risso [2] but I am concerned about how this will affect the contrast matrix? They appear as 3 extra rows compared to the desired contrast matrix (to have 5 rows and 5 columns). Is there some wrong application here? W is a matrix of 3 columns (W_1, W_2, W_3)

design=model.matrix(~0+W+dge$samples$Stage)
colnames(design)[1] ="W_1"
colnames(design)[2] ="W_2"
colnames(design)[3] ="W_3"
colnames(design)[4] ="Control"
colnames(design)[5] ="Stage1"
colnames(design)[6] ="Stage2"
colnames(design)[7] ="Stage3"
colnames(design)[8] ="Stage4"

###Change columns order###
design <- design[, c(5,6,7,8,4,1,2,3)]

v = voomWithQualityWeights(dge, design = design, plot = TRUE)

vfit <- lmFit(v, design)

contrast.matrix <- makeContrasts(Stage1vsControl=Stage1-Control,
Stage2vsControl=Stage2-Control,
Stage3vsControl=Stage3-Control,
Stage4vsControl=Stage4-Control,
levels = colnames(design))

vfit <- contrasts.fit(vfit, contrasts=contrast.matrix)

efit <- eBayes(vfit)

plotSA(efit,main = "Final model: Mean-variance trend")

summary(decideTests(efit, method = "separate", adjust.method = "BH", p.value= 0.05,lfc = 1))

RUVSeq limma DifferentialExpression • 164 views
3
Entering edit mode
@gordon-smyth
Last seen 54 minutes ago
WEHI, Melbourne, Australia

Adding W to the design matrix makes no difference to the contrasts. If you were using a linear model like

design <- model.matrix(~ my.factors)


you just add W to the end of the linear model:

design <- model.matrix(~ my.factors + W)


All your contrasts remain exactly the same.

Yes, the W columns will naturally appear as rows in the contrast matrix, because all the coefficients of the linear model appear as rows of the contrast matrix, but that causes no complications.

By the way, I strongly advise against using logFC cutoffs such as lfc=1. I never do this myself in any of my own published analyses. The help page for decideTests says:

Note

Although this function enables users to set p-value and lfc cutoffs simultaneously, this combination criterion is not recommended. logFC cutoffs tend to favor low expressed genes and thereby reduce rather than increase biological significance. Unless the fold changes and p-values are very highly correlated, the addition of a fold change cutoff can increase the family-wise error rate or false discovery rate above the nominal level. Users wanting to use fold change thresholding are recommended to use treat instead of eBayes and to leave lfc at the default value when using decideTests.

0
Entering edit mode

Thanks for clarifying