Dear Aaron, thank you for your quick confirmation !! So, i will need to rename the two specific column names in the design matrix, in order to not get any error produced by makeContrasts().

Regarding your extra suggestion: you meant by using duplicateCorrelation and avoid including the pairs factor in the design matrix, i can implement some direct comparisons, like directly comparing right-sided vs left-sided cancer samples right ? but this obviously will change the way variance is estimated, as this approach will model/estimate variability in absolute expression values between my patients, correct ? But is obviously will depend on the number of observations, that is my number of patients, which in this case in 30 paired observations(60 samples in total)..

Dear Aaron, sorry to return to this post, but i would like to ask you two very important "updated" questions and how i could interpret them !!

1) In detail, after i also used the duplicateCorrelation about the direct comparison of Right_sided vs Left_sided tumors(15 vs 15), in which (with also the implementation of arrayWeights & duplicateCorrelation), i resulted in a final list after topTable, with 27 DE genes with an adjusted p-val < 0.05. However, only 18 of these had an absolute value of logFC > 0.5, while the others had an absolute logFC from ~3.5 to 4.9. On the other hand, when i used these 28 genes to create a heatmap--also using only my cancer samples--there was a great separation of left and right cancer samples (with 1-2 mislabeled samples). Thus, what is your opinion on this matter ? That despite the fact that in this specific case there is not a great amount of differential expression (and possibly notable changes-but this could be attributed to also various reasons) between "right" and "left" cancer samples, could these 28 genes still considered "valid candidates" for follow up-studies or further research ? As they also show a relatively good separation in my heatmap plot ?

2) Finally,  i would like to ask you about the two lists i have also acquired from limma, regarding the comparison of each tumor anatomic location vs their adjucent controls:

So, my crusial question is: in order to identify DE genes that are unique to each Location (vs their paired controls), i should:

I) From topTable already acquired for specific Location -i.e. Right-sided tumors vs their adjucent controls:

get the DE genes. So for this specific purpose, i should consider only the adjusted p-value cutoff ? i.e. adjusted p.val < 0.05 ? Or i could also consider a further cutoff ? like Treat ? with something like lfc=0.5 ? But with this i could loose genes for the intersection below?

II) And then, in the other list (Left sided tumors vs controls), use in the topTable (after using confint=0.95) 3 simultaneous cutoffs:

top[which(top$CI.L >-1 & top$CI.R<1 & top$adj.P.Val > 0.05),] # the topTable without sorting and all returned probesets

And finally, with the intersection of these two above lists(probesetIDs, as they are in the sameExpressionSet), i will have identify the non-DE genes in the second list ? That is in my example, the Left_sided list?

Best,

Efstathios

Dear Aaron, one last comment regarding my yesterday questions---based also on a useful comment from a previous post (C: Pathway analysis of differentially methylated CpGs) regarding the use of argument robust---i tried it also in my direct comparison between tumors, to see if a get a bigger number of DE genes(from the 27 genes i reported). And then, i saw when i used robust=TRUE, the number of DE genes was reduced slightly from 27 to 24. I acknowledge that this change is possibly negligible, but essentially robust=TRUE should not increase the power for "more DE" ? or my notion is incorrect ? The code used is below:

pairs <- factor(rep(1:30, each=2))
g <- factor(paste0(eset.filtered$Disease, eset.filtered$Location))
design <- model.matrix(~0 + g)
colnames(design) <- levels(g)
aw <- arrayWeights(eset.filtered, design)
barplot(aw, xlab="Array", ylab="Weight", col="blue", las=2) # to ispect indeed the quality of the samples
abline(h=1, lwd=1, lty=2, col="red") 
w <- asMatrixWeights(aw, dim(eset.filtered))
dupcor <- duplicateCorrelation(eset.filtered, design, block=pairs, weights=w)
fit <- lmFit(eset.filtered, design, block=pairs, correlation=dupcor$consensus, weights=w)
cm <- makeContrasts(Direct_loc_comp=Cancerright_sided-Cancerleft_sided,levels=design)
fit2 <- contrasts.fit(fit, cm)
fit3 <- eBayes(fit2, trend=TRUE, robust=TRUE) # the only difference from running the same chunk code is the usage of robust=TRUE
top2_RvsL <- limma::topTable(fit3, coef=1, number=nrow(fit3), adjust.method="fdr", sort.by="none", confint=0.95)

Best,

Efstathios