i would like to adress some specific design implementations in limma and ask your opinion about "valid" or appropriate they are for my comparison of interest. It refers to a cancer microarray dataset, and the phenoData looks like:
Meta_factor Disease Location Study
St_1_WL57A.CEL 0 Normal left_sided hgu133plus2
St_2_WL57A.CEL 0 Cancer left_sided hgu133plus2
St_N_EC59A.CEL 0 Normal left_sided hgu133plus2
St_T_EC59A.CEL 0 Cancer left_sided hgu133plus2
St_N_EJ31A.CEL 0 Normal right_sided hgu133plus2
St_T_EJ31A.CEL 0 Cancer right_sided hgu133plus2
As i have adressed in some previous posts, my initial goal was to identity DE genes regarding different anatomic tumor location, versus their adjucent controls. Thus, i used an interaction model, as my samples are paired-belong to each patient:
pairs <- factor(rep(1:30, each=2))
design <- model.matrix(~ 0 + pairs + Disease:Location, data=eset.filtered)
design1 <- design[,-grep("Normal", colnames(design))] # get full rank
arrayw <- arrayWeights(eset.filtered, design=design1)
fit <- lmFit(eset.filtered, design1, weights=arrayw)
fit2 <- eBayes(fit, trend=TRUE)..
But now, my collaborators asked me if i can through limma to test directly for genes differentially expressed between right_sided and left_sided tumors, and not just comparing the two DE resulted lists for common genuine or diffenent DE genes (which is also of the same importance):
So, from the above design matrix, i have two coefficients DiseaseCancer:Locationleft_sided & DiseaseCancer:Locationright_sided, which represent the average log2-fold change between the left sided tumors vs their adjucent controls, and also the same for the right-sided tumors respectively. I understand, that trying to compare directly the two above coefficients, would probably give no results, as any differences would mostly absorved in the pairs factor. Thus, for my above question, after the above step of fit, by implementing:
cm <- makeContrasts(RvsL= DiseaseCancer:Locationright_sided-DiseaseCancer:Locationleft_sided, levels=design)
fit2 <- contrasts.fit(fit, cm)
fit3 <- eBayes(fit2, trend=TRUE)...
and with this contrast i can ask in which genes the log2-fold change of the right-sided tumors vs their controls, is significantly different in the left_sided tumors ? And this could at "least" partially give some insights in my question ?
Or my approach is incorrect ?