Dear Bioconductor Community,

i have returned with a previous problem i have faced from a previous post (Possible implementation of multifactorial contrasts in limma regarding a microarray dataset). To provide a small summary from the link, with the great assistance and crusial suggestions of Aaron Lun i implemented in my merged colorectal affymetrix dataset, comprised of 30 different patients with two samples each(cancer and adjucent control), a multifactorial design in limma, incoprorating except from the "paired comparison" the specific tumor location of each patient:

head(pData(eset.2))
                                   Disease      Location                  Meta_factor        Study
St_1_WL57.CEL           Normal       sigmoid_colon           0                      hgu133plus2
St_2_WL57.CEL           Cancer       sigmoid_colon           0                      hgu133plus2
St_N_EC59.CEL           Normal       sigmoid_colon           0                      hgu133plus2
St_T_EC59.CEL           Cancer       sigmoid_colon           0                      hgu133plus2
St_N_EJ58.CEL           Normal                cecum              0                      hgu133plus2
St_T_EJ58.CEL           Cancer               cecum               0                     hgu133plus2

with the following code: 

pairs <- factor(rep(1:30, each=2))
design <- model.matrix(~ 0 + pairs + Disease:Location, data=eset.2)
design <- design[,-grep("Normal", colnames(design))] # get full rank
colnames(design)
fit <- lmFit(eset.2, design)
fit2 <- eBayes(fit, trend=TRUE)

Moreover, as Aaron mentioned, except from the first 30 coefficients which represent to the patient-specific blocking factors, the coefficients 31:39 correspond to the location-specific log-fold change of cancer vs control samples.

Thus, for istanse top2 <- topTable(fit2, coef="DiseaseCancer:Locationascending_colon", number=nrow(fit2), adjust.method="fdr", sort.by="none")  represents the log-fold change of the 8 ascending colon tumor samples versus their adjucent controls

So, to adress and highlight my problem, although there are 9 coefficients which represent 9 different anatomic locations, only in the above(ascending colon) i get  from topTable DE genes with an FDR < 0.05.

table(pData(eset.2)$Location)

    ascending_colon               cecum    descending_colon  left_colic_flexure 
                 10                  12                   6                   2 
 rectosigmoid_colon              rectum right_colic_flexure       sigmoid_colon 
                  4                   8                   6                  10 
   transverse_colon 
                  2 

As you can see from the above, although some comparisons include a very low number of samples, in the other comparisons it is very weird (like cecum and sigmoid colon) not to return any DEG genes. I acknoledge that each study has its specific details and experimental design regarding the biology of each system, as also other various reasons related with the analysis itself, but as it is a remaining issue i also discussed with my lab members and remains under discussion, i would like to adress some more questions(as im mainly a molecular biologist and perform data analysis about a year) to get also a statistical view or suggestions on the matter:

1. As Aaron also suggested various possible explanations, i would like first to highlight one of them(and please correct me for any misunderstanding or misconception): how it is possible(from a statistical view) ,a high variability in a specific location(although patients are different) to mask the DE genes expression in the other locations ? Even if generally cancer vs control samples generally yield a lot of different gene expression patterns ? And furthermore, as this model Aaron suggested models the variability in the cancer/normal logFC between the patients ?

2. One other suggestion was to construct an MDS plot but based on the log-fold changes of each of the above nine coefficients. Thus, i have to use something like this:

q <- fit2$coefficients[,31:39] and then run: plotMDS(x=q, top=500, gene.selection="pairwise") ?

And if the above implementation is correct, i should inspect the plot for samples that "behave" as outliers or gather far from the others ?

3. To continue from the question 3, then i should remove the specific patient samples from my expression set, and then re-run the analysis ? 

I know also the function removeBatchEffect() but i have never use it and i dont know if an how fits for my current problem ?

4. Finally, there are any other diagnostic plots or ways to adress the problem, or alternative solutions ?

5. To mention, before performing limma, from my EDA plots like clustering the samples and PCA, one specific control sample seems to be a common outlier. Also, i have perfomed paired analysis with limma without the location information

(*Please excuse me for any naive questions and my long message, but as im generally new in R and i dont have many experience in these specific issues and statistical models, i would like to cover any important part of my analysis and also to gain as much as possible feedback !!)

1. You never 'control' variance. You use variance as a yardstick to say if a given fold change is big enough to reliably say the underlying population difference is different from zero. The larger your variance, the bigger the fold change has to be. But do note that we are using sample variance as our estimator, and sample variance is not a very efficient statistic. So it takes more data for the sample estimate to converge to the underlying population parameter we are trying to estimate. This is why an ANOVA model with multiple groups is more powerful than a t-test; you borrow strength from all the data in the model. But if one group is much more variable, then it may undermine the results. Or maybe that group is more reliably estimating the true variance, and the other groups have lower than expected variance.
2. Yes. You estimate weights and then use them in the lmFit() step.
3. Methods like WGCNA don't accept sample weights, so you have to do other things. Removing outliers is one option.
4. I'll let Aaron answer questions about suggestions he has made.

I applaud your interest in learning statistics, and wanting to do things for yourself. But do note that the primary purpose of this support site is not teaching people statistics, but helping people who have technical problems when running Bioconductor packages. I (and Aaron, Gordon, and some others) do sometimes veer off topic to help people with fundamental statistical questions, but you should not consider this to be an on demand Coursera teaching site.

If you really plan to continue analyzing data for your group, you should take it upon yourself to learn enough statistics that you can both competently analyze data, as well as explain to your colleagues what you have done, why you did it, and the tradeoffs you had to make. You cannot expect to get that knowledge here, by asking repeated questions about each analysis you undertake. Both Rafael Irizarray's group (on edX) and Jeff Leek/Roger Peng/Brian Caffo (on Coursera) have online classes that you can take. There is also Julian Faraway's R modeling book that you can get free on r-project.org. Those are more appropriate ways to get the knowledge you seek.

1. When you fit an ANOVA model and compute contrasts, the denominator of your statistic is based on (essentially) the average within group variability, using all groups. In one sense this is better, as you are using more data to estimate variance, and all else equal, more data is better than less data for estimating variance. However, if you have one group that has an unexpectedly high within-group variance, it may inflate the variance estimate 'too much'.
2. If you have highly variable samples, you might consider using a weighted model fit. Sometimes you have enough samples that you can afford to remove outliers, but that is pretty ad hoc, and usually if you remove one sample, then another starts looking like an outlier, so you remove that one too, but then there is another outlier...
3. Like I said, using sample weights is often a better idea than seeking and removing outliers. removeBatchEffect() is only intended to remove known batch effects, and in that case only for plotting purposes. If you have known batches you should either fit as part of the linear model or use sva.