Entering edit mode
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.
Have you read the limma user guide? If so, post the code you have tried. If not, that is the place to start.
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?
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)`
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.
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,