heatmap for differential expressed genes

Question

[minni9234]

how to generate a heatmap after the limma analysis (for differential expressed genes)?

fit<-lmFit(selExpr,design)
contrast.matrix<-makeContrasts(BT20_1_day-BT20_0, levels=design)
fits<-contrasts.fit(fit,contrast.matrix)
ebfit<-eBayes(fits)
topTable(ebfit,number=10000,coef=1,adjust="fdr",p.value=0.05)

[Devon Ryan]

Have you checked the posts on the right of the screen? Just save the data.frame produced by topTable without sorting, subset the expressionSet accordingly, and make a heatmap of that.

[minni9234]

I tried that, and it was not working

[Devon Ryan]

Define "not working".

[minni9234]

It keeps showing me an error even I have already set it to be a numeric matrix

Error in heatmap.plus(top) : 'x' must be a numeric matrix

[Devon Ryan]

It would be rather more useful if you also showed the entirety of the relevant code...

[minni9234]

> top<-as.numeric(rownames(topTable(ebfit,n=10000,coef=1,adjust="fdr",p.value=0.05)))

>heatmap.plus(top)

[Devon Ryan]

I'm not sure how you could hope for that to ever work. You're trying to make a heatmap of probe names, which doesn't make any sense.

tt <- topTable(ebfit,n=Inf,coef=1,adjust="fdr",p.value=0.05, sort.by="none")
idx <- which(tt$fdr <0.1)
subEset <- selExpr[idx]
heatmap(assayData(subEset))

or something like that.

[karl.stamm]

as.numeric of the rownames? Stopped reading there.

Answer 1

[Sean Davis]

fit<-lmFit(selExpr,design)
contrast.matrix<-makeContrasts(BT20_1_day-BT20_0, levels=design)
fits<-contrasts.fit(fit,contrast.matrix)
ebfit<-eBayes(fits)
# use top 200 probes
tt = topTable(ebfit,number=200,coef=1,adjust="fdr",p.value=0.05)
library(gplots)
idx = rownames(tt)
# assume that selExpr is an ExpressionSet
heatmap.2(exprs(selExpr)[idx,],trace='none',scale='row')
# OR, if selExpr is a matrix
heatmap.2(selExpr[idx,],trace='none',scale='row')

[minni9234]

Hi, Sean Davis. Thanks for your useful explanation.

However, I have this same errors again.

> heatmap.2(exprs(selExpr)[index_top,],trace='none',scale='row')
Error in function (classes, fdef, mtable)  : 
  unable to find an inherited method for function "exprs", for signature "matrix"

[Devon Ryan]

Note the "assuming that selExpr is an ExpressionSet" comment... Try the next non-comment line.

[minni9234]

Sorry, now I have another error.

> heatmap.2(selExpr[idx,],trace='none',scale='row')

Error in heatmap.2(selExpr[idx, ], trace = "none", scale = "row") : 
  subscript out of bounds

[Devon Ryan]

Change the tt=... and idx=... lines to these:

tt = topTable(ebfit,n=Inf,coef=1,adjust="fdr", sort.by="none")
idx = which(tt$adj.P.Value <0.1)

The original version from Sean Davis won't work if the rownames are probe ID numbers.

[minni9234]

Hi Devon, this doesn't work either.

> heatmap.2(selExpr[idx,],trace='none',scale='row')
Error in heatmap.2(selExpr[idx, ], trace = "none", scale = "row") : 
  `x' must have at least 2 rows and 2 columns

Here is the result I got from topTable:

     row.names    ID    logFC    AveExpr    t    P.Value    adj.P.Val    B
1    10603    HSD3B1    -7.0675833    5.68016    -48.837013    3.100100e-22    1.341568e-17    38.539428
2    40338    LOC100505633    7.2969167    3.78948    45.883284    1.064792e-21    1.616203e-17    37.668604
3    14880    DUSP6    -6.0338333    7.95048    -45.765236    1.120418e-21    1.616203e-17    37.631914
4    11308    HCAR3    5.3886667    3.31332    38.487688    3.416048e-20    3.695737e-16    35.041866
5    9566    COL15A1    -5.3935833    5.30452    -33.774838    4.453433e-19    3.854446e-15    32.944817

[Devon Ryan]

Well, adj.P.Value should just be adj.P.Val, then.

[minni9234]

Thanks a lot. I got it

[Vasu]

Hi Devon,

I have used the same to get heatmap for differential expressed genes. But, got an error.

design <- model.matrix(~ -1+factor(c(1,1,2,2,3,3)))
colnames(design) <- c("group1", "group2", "group3")
fit <- lmFit(celfiles.rma, design)
contrast.matrix <- makeContrasts(group2-group1, group3-group2, group3-group1, levels=design)
fit2 <- contrasts.fit(fit, contrast.matrix)
fit2 <- eBayes(fit2)
tab <- topTable(fit2, coef=1, adjust="fdr", sort.by="B", number=Inf)
head(tab)
             logFC   AveExpr          t      P.Value adj.P.Val         B
17100820  3.813465 10.084444  13.990013 1.477686e-05 0.3544622 -4.142070
17100824  3.813465 10.084444  13.990013 1.477686e-05 0.3544622 -4.142070
16744689 -3.037688  1.577746 -13.252989 1.983301e-05 0.3544622 -4.143913
16842465 -6.954801  4.096664 -11.968499 3.444300e-05 0.4616826 -4.147931
16843156  3.222599  2.013231  11.375102 4.530001e-05 0.4857702 -4.150236
16706416 -2.014341  2.072348  -8.395743 2.270946e-04 1.0000000 -4.169346

idx = which(tab$P.Val < 0.9)
heatmap.2(tab[idx,],trace='none',scale='row')

Error in heatmap.2(tab[idx, ], trace = "none", scale = "row") : 
  `x' must be a numeric matrix

Can you help me in this?

[Devon Ryan]

Make a heatmap of the normalized counts, not the output of topTable().

[Vasu]

Thanks for the reply Devon. I have used oligo for my Affymetrix Human gene 2.0 ST array.

head(tab)
               logFC   AveExpr           t    P.Value adj.P.Val         B
16650001 -0.67713091 0.7013795 -2.11246203 0.07463104 0.8414740 -3.923065
16650003  0.01775323 0.8173575  0.05752413 0.95581904 0.9922585 -5.247542
16650005  0.20727356 1.2926154  0.28654513 0.78318782 0.9636612 -5.214740
16650007  0.13198594 0.8028477  0.37395467 0.72008320 0.9502862 -5.191016
16650009 -0.31117759 0.7187522 -0.92432176 0.38765340 0.8815813 -4.917247
16650011 -0.18142907 0.7450577 -0.46483921 0.65689248 0.9382330 -5.160082

results <- tab[which(tab$logFC >= 1.5 & tab$P.Value <= 0.05),]

idx = which(tab$P.Value < 0.05 & tab$logFC > 1.5)

heatmap(exprs(celfiles.rma[idx,]),trace='none',scale='row')

But I'm getting the Unit numbers. How to add annotation to get the gene names for this?

[Devon Ryan]

No clue on that part, I haven't used microarrays in a number of years.

[Vasu]

Ok. Thank you very much !!

[karl.stamm]

Do you have a question that we can answer? When you dump code on us and ask "how to do" there is no specific thing we can give you. You are asking the internet to do your job. Search "making heatmaps in R" on the google or on biostars and you'll see many tutorials. You have then to read them and try it. When you get a specific question, search that, and we may be able to help.