Dear Bioconductor community,

I am currently struggling with fitting a linear model to Illumina BeadChp HT12.v4 data. I have different cell generations (G1, G2, G3 etc) taken form the same clone. I assume these as biological replicates. For each generation I have different measurements. I assume these as technical replicates. Microarray data were measured in two different batches. Finally, I have a specific mutation that I am looking into: T790M: this is present in some of the samples and not present in others. Part of my targets file looks like this:

 I want to find differentially expressed genes between T790M positive (pos) and T790M negative (neg) samples. I am following limma's user guide for analysis of BeadArrays (section 17.3). I have reached the state where I need to factorize by my replicates. The guide says:"The study involves multiple cell types for the same patient..."

ct <- factor(targets$CellType)

and then the DE genes between these cell types are detected.

In my case it is different. I tried to factorize based on the Sample and on the CellClone and include the batch in my design matrix:

ct<-factor(targets$Sample)  
design <- model.matrix(~0+ct+targets$Batch)

but I realized that if I follow the example given by limma, I then have to detect the differences of the variables that rcorrespond to the colnames of my design matrix, which would be CL75_G.2, CL86_G5 etc..This is not the comparison that I want.

I then tried

ct<-factor(targets$T790M)  
design <- model.matrix(~0+ct+targets$Batch)
dupcor <- duplicateCorrelation(gefR, design, block=targets$Sample) #cor between samples of the same clone
dupcor$consensus.correlation
[1] 0.0723

then give the pos, neg and batch as column names of my design matrix and try to detect the DE genes between T790M pos and T790M neg, using the suggested procedure.

fit <- lmFit(myMatrix, design, block=targets$Sample, correlation=dupcor$consensus.correlation)
# the block indicates the technical replicates.
contrasts <- makeContrasts('pos-neg',levels=design)
fit2 <- contrasts.fit(fit, contrasts)
fit2 <- eBayes(fit2)
summary(decideTests(fit2, method="global"))

The outcome was 0 DE genes (4597 are the total genes that are truly expressed on my array, according to their  Detection p-value)

   pos-neg
-1       0
0     4597
1        0

I cannot imagine this result is true. I am sure I am missing something but I cannot understand what... The level of complexity of my experiment is grater than the complexity of the limma's example...I would really appreciate any comments/answers/insights from you.

Thank you very much!

Eleni

Your second approach seems to be correct. The only thing that's odd is that you run duplicateCorrelation on the gefR matrix, but lmFit on the myMatrix object. These had better hold the same data.

Anyway, there's a couple of reasons for why you don't have any DE genes. The first is that, in the subset of the targets file you've shown, there are only 4 positive samples and all of these belong in one batch. If there are no positive samples anywhere else, this means that the DE comparison is only being performed using the positive/negative samples in the first batch (as the Batch blocking factor will absorb the rest). This will reduce power even if you have lots of negative samples in the other batches. It's worth checking on a MDS plot whether the batch effect really exists; if it doesn't, you could improve power by dropping it from the model.

Another possible cause is outlier effects, either in terms of samples or gene-specific variances. The former you can identify by looking at a MDS plot (probably after running removeBatchEffect to get rid of any Batch effect, if one is present), and can be dealt with using arrayWeights. The latter you can protect against with robust=TRUE in eBayes.

If none of these suggestions result in DE genes, then either the mutation doesn't result in any DE, or your data set simply does not have enough statistical power to detect it. In such cases, it's often nice to have some positive controls (i.e., genes you know to be DE upon perturbation) and see how they're behaving in order to diagnose potential problems.