Dear Bioconductor Community,

in a previous thread( https://support.bioconductor.org/p/64372/#64413) i have asked about the implementation of a paired analysis, and Mr. James W. MacDonald kindly provided me an suggestion of a paired limma analysis, which i used on the code below :

library(limma)
conditions <- data.trusted.eset$condition
condition <- factor(conditions, levels(condition)[c(2,1)])
pairs <- factor(rep(1:13, each = 2))
design <- model.matrix(~condition+pairs)
fit <- lmFit(data.trusted.eset, design)
fit2 <- eBayes(fit)

In the next step i wanted to implement a multifactorial analysis in limma to evaluate the interaction between two factors, condition(control & cancer tissue) & meta_factor(metastatic & non-metastatic cancer), in order to evaluate possible DE genes between metastatic cancer and non-metastatic cancer samples:

f <- paste(data.trusted.eset$condition, data.trusted.eset$Meta_factor, sep=".")
f <- factor(f)
pairs <- factor(rep(1:13, each = 2))
design <- model.matrix(~0 +f +pairs)

fit <- lmFit(data.trusted.eset, design)
head(names$coefficients)
contrast.matrix <- makeContrasts(M.CvsC="fCancer.1-fCancer.0", M.CvsControl="fCancer.1-fNormal.0", levels=design)
fitb <- contrasts.fit(fit, contrast.matrix)
fit2 <- eBayes(fitb)

My first question is if it is wrong for the statistical model that i didnt use an intercept term(because with the intercept term i couldnt use the coefficient "fCancer.0" for the contrast.matrix) ? And secondly, is it also important as with the first analysis with just one factor, to include the term pairs to get meaningful results ? My last question is that with topTable, i get only 4 genes with adj.p,Val < 0.05. Is that necessarily wrong or it might illustrate possible small differences between metastatic and non-metastatic tumors ?

Thank you for your patience(& please excuse me for any "begginer" questions, as i have only using R & Bioconductor for a few months).

You can no longer use a paired design, because Meta_factor is a between-subject factor and paired analyses only deal with within-subject comparisons. Instead you have to treat subject (pair) as a random effect, which is done in limma using the duplicateCorrelation function.

The limma User's Guide calls your experiment a "Multi-level experiment". Just look up Section 9.7 in the User's Guide and follow the instructions.

Your current design matrix is a bit tricky to interpret. I'd suggest doing something like this:

> f <- relevel(f, ref="Normal.0")
> design <- model.matrix(~0 + pairs + f)
> design <- design[,-ncol(design)]
> colnames(design)
 [1] "pairs1"    "pairs2"    "pairs3"    "pairs4"    "pairs5"    "pairs6"   
 [7] "pairs7"    "pairs8"    "pairs9"    "pairs10"   "pairs11"   "pairs12"  
[13] "pairs13"   "fCancer.0" "fCancer.1"

The first 13 coefficients describe the patient blocking effect. Coefficient 14 describes the fold change of Cancer.0 over Normal.0. Coefficient 15 describes the fold change of Cancer.1 over Normal.1.

Now, you can't directly compare between Cancer.1 and Cancer.0, because any differences will be absorbed into the patient blocking factors (this is also true of your original design). Instead, you can ask whether the log-fold change  between metastatic cancer and paired normal samples is equal to that between non-metastatic cancer and their normal samples, by using:

contrast.matrix <- makeContrasts(M.CvsC="fCancer.1-fCancer.0", level=design)

Of course, you can also compare between each cancer and its normal sample directly.

If you don't choose to include the pairing in design, you're going to be treating different patients as biological replicates. This will likely result in large dispersion estimates, which will confound detection of DE genes. On the other hand, it means that your dataset simplifies to a three-group design, so you can compare between metastatic and non-metastatic cancer samples directly.

As to the small number of DE genes; that depends on the variance of your data, the size of the fold-changes between normal and cancer samples for each patient, and the parametrization of your design and contrast matrices. For example, the contrast between Cancer.1 and Cancer.0 is unlikely to pick up any DE genes, as any changes will be absorbed by the blocking factors.