Microarray analysis beadarray with limma gives no differential expression

Question

[dharl]

I am fairly new to microarray analysis. I have 40 arrays from four phenotypes and would like to compare the gene expression among the four subtypes.

For this I used bead array (since the arrays are illumina based), preprocessed with BASH/HULK and summarised it and ran neqc for background correction and normalisation and finally have log2 expression data as an ExpressionSetillumina.

I then use this for carrying out differential expression analysis using Limma. I would like to compare the gene expression of one subtype with the others. For this I first use the design matrix and define the arrays corresponding to a subtype with one and the rest 0. I then compare the groups. Note that there are different number of repilcates in each subgroup.

Sadly I do not get any differential expression which is very unlikely and I am sure there is something wrong.

Any comments?

[David Westergaard]

To see if anything is wrong, you will have to post the code used. From your description, it is not obvious just how you made your design matrix, and which coefficient you are checking against.

[dharl]

Hi, I agree... please see the comments I have posted it there.

[dharl]

Here is what I have done: Phenos is the phenotypes whose RNA is measured and Eset is the expression set. PLease dp let me know if any relevant information is missing.. Here I compare the Pheotypes in pairs.. I also tried comparing one Phenotype to the other three. That also did not work. I also tried different normalizations apart from neqc.

---

table(Phenos)
Phenos
1  2  3  4
12  8 10 10 

Design<-model.matrix(~0+Phenos)

 Phenos1 Phenos2 Phenos3 Phenos4
1       0       1       0       0
2       1       0       0       0
3       0       0       0       1
4       0       1       0       0
5       0       0       0       1
6       1       0       0       0

Eset<-exprs(BSData.lumi.filt)[,rownames(Design)]

contrast.matrix <- makeContrasts(contrasts = c("Phenos1-Phenos2","Phenos1-Phenos3","Phenos1-Phenos4","Phenos2-Phenos3","Phenos2-Phenos4","Phenos3-Phenos4"), levels = Design)
aw<-arrayWeights(Eset,Design)
fit <- lmFit(Eset,Design,weights=aw)
fit2 <- contrasts.fit(fit, contrast.matrix)
Fit <- eBayes(fit2)

topTable(Fit,adjust='BH')
             Phenos1.Phenos2 Phenos1.Phenos3 Phenos1.Phenos4 Phenos2.Phenos3 Phenos2.Phenos4 Phenos3.Phenos4  AveExpr        F      P.Value adj.P.Val
ILMN_1718042     0.224781478     -0.08038455    -0.024472954     -0.30516603     -0.24925443      0.05591160 5.084373 9.621334 7.559494e-05  0.955246
ILMN_2375879    -0.197786684     -0.03969275     0.156960307      0.15809394      0.35474699      0.19665306 4.765450 8.860060 1.428595e-04  0.955246

[David Westergaard]

I am not familiar with Illumina beadchip arrays, but have you checked if any of the arrays are outliers or weird, and how the distributions look after normalisation?

Also, the results from topTable looks very odd. By default, it should only check coef=1, i.e. Group1 vs Group2. Maybe you should also try asking on the Bioconductor mailing list.

[dharl]

Ya I checked for outliers... I think I will ask the mailing group as you have suggested.

[dharl]

Well, turns out, like Sean said, it was lack of differential expression within my samples. Thanks a lot for all your help.

Answer 1

[Sean Davis]

You'll need to use four groups, not just 2. Then, to define differential expression, you'll need to define a contrast matrix. See the limma user guide for an example of single-color arrays with multiple groups.

[dharl]

Hi

I did as you have said. I have posted the code in the comments above. Thanks again.

[Sean Davis]

The code looks good. There are many reasons why you might not find significance, including lack of biological difference and lack of power (not enough samples). You can always try ranking your genes and then use another technology to try to validate some differences if sample size is limited.

[dharl]

Ya I was thinking of looking more into the data and see if any transformation can be applied.. Can you suggest any method which does this ranking?

[Sean Davis]

Just rank the genes by p-value. There is no "method", in that sense. No transformation is needed.