Question: Microarray analysis beadarray with limma gives no differential expression
2
gravatar for dharl
4.5 years ago by
dharl70
India
dharl70 wrote:

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 donot get any differential expression which is very unlikely and I am sure there is something wrong. 

Any comments?

illumina limma beadarray R • 2.1k views
ADD COMMENTlink modified 4.5 years ago • written 4.5 years ago by dharl70
1

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.

ADD REPLYlink written 4.5 years ago by David Westergaard1.4k

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

ADD REPLYlink written 4.5 years ago by dharl70

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
ADD REPLYlink written 4.5 years ago by dharl70

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. 

ADD REPLYlink written 4.5 years ago by David Westergaard1.4k

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

ADD REPLYlink written 4.5 years ago by dharl70
1
gravatar for Sean Davis
4.5 years ago by
Sean Davis25k
National Institutes of Health, Bethesda, MD
Sean Davis25k wrote:

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.

ADD COMMENTlink written 4.5 years ago by Sean Davis25k

Hi 

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

 

ADD REPLYlink written 4.5 years ago by dharl70

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.

ADD REPLYlink written 4.5 years ago by Sean Davis25k
0
gravatar for dharl
4.5 years ago by
dharl70
India
dharl70 wrote:

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? 

ADD COMMENTlink written 4.5 years ago by dharl70

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

ADD REPLYlink written 4.5 years ago by Sean Davis25k
0
gravatar for dharl
4.5 years ago by
dharl70
India
dharl70 wrote:

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

ADD COMMENTlink written 4.5 years ago by dharl70
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1097 users visited in the last hour