DEGs of microarray Data and hierarchical clustering and visualizing using heatmap
1
0
Entering edit mode
@lawardeankita1-14844
Last seen 2.5 years ago

Hello Sir/madam,

I have microarray data from the GEO database  (GSE51401) for which i have to find the DEGs and to generate a heatmap of those DEGs.

I am using .CEL files (raw data) for the analysis using limma package. and the data has NA values, so can anyone tell me how to remove the NA using limma and then the R code for the DGEs.

limma heatmap DEGs • 980 views
ADD COMMENT
1
Entering edit mode

Have you read the limma user guide? If so, post the code you have tried. If not, that is the place to start.

ADD REPLY
1
Entering edit mode

Hello Sir,

using limma package the following R code i have used

library(affy)
library(limma)
library(gplots)

pd <- read.AnnotatedDataFrame("tnt.txt",header = T,sep = "\t")
pData(pd)
my_expressionset <- ReadAffy(filenames =rownames(pData(pd)))
my_expressionset
eset <- rma(my_expressionset)
pData(eset)

# design matrix
design <- model.matrix(~Substance,pData(pd))
design
fit <- lmFit(eset, design)
fit
fit2 <- eBayes(fit)
fit2
dim(fit2)
options(digits = 2)
fit2$p.value
deg <- topTable(fit2, coef = 2, p.value = 0.001, adjust.method = 'fdr', number = nrow(eset))
dim(deg)[1]

and I'm getting total 5201 deferentially expressed genes, but i am still confused whether its right or wrong.

and i have to plot a heatmap of these DEGs, 

for wich i followed one biostar question answer help, the code is as follows for heatmap

##heatmap of DEG
idx = rownames(deg)
head(idx)
heatmap.2(exprs(fit2)[idx,],trace='none',scale='row')

but it is showing following error,

Error in (function (classes, fdef, mtable)  : 
  unable to find an inherited method for function ‘exprs’ for signature ‘"MArrayLM"’

can you tell me how to make heatmap after limma and is my number of DEGs is correct or not?

 

ADD REPLY
0
Entering edit mode

Having code is really one of the best ways to get a clear answer. To get your heatmap, you might want to try using `exprs(my_expressionset)` rather than `exprs(fit2)`

ADD REPLY
0
Entering edit mode

Hello Sir,

Thank you for the reply and the answer,

below i have made changes to the R code for the heatmap generation i have used for DEG analysis, along with the result of each step.

Expression data from tumor-associated endothelial cells of hepatocellular carcinoma GEO accession, (GSE51401) and there are total 24 Tumar and 24 Non tumor samples (total 48 samples) 

1. my first question is, is this code is right ?

2. i am getting total 5201 DEGs, so is this number of DEGs is right or am i doing any mistake in the code ?

3. do i have to make contrast matrix for the DEG analysis?

4. i want to take only unregulated genes from the DEGs.

library(affy)
library(limma)
library(gplots)

pd <- read.AnnotatedDataFrame("tnt.txt",header = T,sep = "\t")
pData(pd)

                                      Substance
GSM1244729_BH2010140_5_NT105.CEL            TEC
GSM1244731_BH2010140_5_NP105.CEL            NEC
my_expressionset <- ReadAffy(filenames =rownames(pData(pd)))
my_expressionseteset <- rma(my_expressionset)

pData(eset)

                                     sample
GSM1244729_BH2010140_5_NT105.CEL           1
GSM1244733_BH2010140_5_OT105.CEL           2

# design matrix

design <- model.matrix(~Substance,pData(pd))

design

                           (Intercept) SubstanceTEC

GSM1244729_BH2010140_5_NT105.CEL                1            1
GSM1244731_BH2010140_5_NP105.CEL                1            0


fit <- lmFit(eset, design)
fitAn object of class "MArrayLM"

$coefficients
          (Intercept) SubstanceTEC
1007_s_at         7.1       -0.025
1053_at           6.4        0.546
117_at            8.5        0.138
121_at            7.0        0.280
1255_g_at         2.1        0.063
54670 more rows ...

fit2 <- eBayes(fit)
fit2

dim(fit2)
options(digits = 2)
fit2$p.value
deg <- topTable(fit2, coef = 2, p.value = 0.001, adjust.method = 'fdr', number = nrow(eset))

head(deg)

            logFC AveExpr   t P.Value adj.P.Val  B
223597_at    -4.4     6.9 -21 4.7e-26   2.5e-21 48
205819_at    -4.8     9.2 -20 7.4e-25   2.0e-20 45
dim(deg)[1]

[1] 5201

##heatmap of DEG
idx = rownames(deg)
eset[idx]
eset2 <- exprs(eset[idx])
head(eset2)
hr <- hclust(as.dist(1-cor(t(eset2), method="pearson")),method="complete")
heatmap(eset2, Rowv = as.dendrogram(hr),Colv = NA)

 

Thank you sir in advance,

 

ADD REPLY
1
Entering edit mode
@gordon-smyth
Last seen 1 hour ago
WEHI, Melbourne, Australia

CEL files do not generate NA values. If you have NA values, then you're doing some wrong.

Even if you did have NA values, then limma would handle them automatically. You wouldn't have to make any intervention yourself.

ADD COMMENT

Login before adding your answer.

Traffic: 220 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6