Question: heatmap for differential expressed genes
4
gravatar for minni9234
5.4 years ago by
minni923450
United States
minni923450 wrote:

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)

ADD COMMENTlink modified 5.4 years ago by Sean Davis25k • written 5.4 years ago by minni923450
1

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.

ADD REPLYlink written 5.4 years ago by Devon Ryan93k

I tried that, and it was not working

ADD REPLYlink written 5.4 years ago by minni923450

Define "not working".

ADD REPLYlink written 5.4 years ago by Devon Ryan93k

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

ADD REPLYlink written 5.4 years ago by minni923450

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

ADD REPLYlink written 5.4 years ago by Devon Ryan93k

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

>heatmap.plus(top)

ADD REPLYlink written 5.4 years ago by minni923450
1

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.

ADD REPLYlink modified 5.4 years ago • written 5.4 years ago by Devon Ryan93k

as.numeric of the rownames? Stopped reading there.

ADD REPLYlink written 5.4 years ago by karl.stamm3.5k
2
gravatar for Sean Davis
5.4 years ago by
Sean Davis25k
National Institutes of Health, Bethesda, MD
Sean Davis25k wrote:
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')
ADD COMMENTlink written 5.4 years ago by Sean Davis25k

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"

ADD REPLYlink written 5.4 years ago by minni923450

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

ADD REPLYlink written 5.4 years ago by Devon Ryan93k

 

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

ADD REPLYlink written 5.4 years ago by minni923450

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.

ADD REPLYlink written 5.4 years ago by Devon Ryan93k

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

ADD REPLYlink written 5.4 years ago by minni923450

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

ADD REPLYlink written 5.4 years ago by Devon Ryan93k

Thanks a lot. I got it

ADD REPLYlink written 5.4 years ago by minni923450

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?

ADD REPLYlink written 3.7 years ago by Vasu410
1

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

ADD REPLYlink written 3.7 years ago by Devon Ryan93k

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 Im getting the Unit numbers. How to add annotation to get the gene names for this?

ADD REPLYlink modified 3.6 years ago • written 3.7 years ago by Vasu410

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

ADD REPLYlink written 3.7 years ago by Devon Ryan93k

Ok. Thank you very much !!

ADD REPLYlink modified 3.7 years ago • written 3.7 years ago by Vasu410

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.

ADD REPLYlink written 3.7 years ago by karl.stamm3.5k
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: 1600 users visited in the last hour