《Epidermal Growth Factor Receptor Activation: An Upstream Signal for Transition of Quiescent Astrocytes into Reactive Astrocytes after Neural Injury》
首先ncbi中检索该文章,得到数据GSE5282
第一步:Start:获取数据并得到表达矩阵
rm(list = ls())
options(stringsAsFactors = F)
library(GEOquery)
eSet <- getGEO("GSE5282",
destdir = '.',
getGPL = F)
class(eSet)
length(eSet)
class(eSet[[1]])
exp <- exprs(eSet[[1]])
exp[1:4,1:4]
exp1 <- exp[,c(4,5,7,1,2)]
exp1[1:4,1:3]
pd <- pData(eSet[[1]])
pd
save(pd,exp,exp1,file = "step1output.Rdata")
eSet[[1]]@annotation
第二步:设置分组信息
##control:EGF4h
m(list = ls())
options(stringsAsFactors = F)
load(file = "step1output.Rdata")
class(pd)
dim(pd)
colnames(pd)
group_list1=rep(c('control','EGF_4h'),c(3,2))
group_list1
save(group_list1,file = "step2output.Rdata")
第三步:检验组间数据的显著性(PCA)
rm(list = ls()) ## 魔幻操作,一键清空~
options(stringsAsFactors = F)
load(file = "step1output.Rdata")
load(file = "step2output.Rdata")
exp1[1:4,1:4]
dat1=as.data.frame(t(exp1))
dim(exp1)
dim(dat1)
dat1=cbind(dat1,group_list1)
library(FactoMineR)
library(factoextra)
dat.pca <- PCA(dat1[,-ncol(dat1)], graph = FALSE)
fviz_pca_ind(dat.pca,
geom.ind = "point",
col.ind = dat1$group_list1,
palette = c("#00AFBB", "#E7B800"),
addEllipses = TRUE,
legend.title = "Groups")
cg=names(tail(sort(apply(exp,1,sd)),1000))
n=t(scale(t(exp[cg,])))
thr=2
n[n>thr]=thr
n[n< -thr]= -thr
n[1:4,1:4]
pheatmap(n,show_colnames =F,show_rownames = F)
Result:
PCA
通过图示,很明显发现两者的组间差异明显,数据合理。
第四步:筛选差异基因
rm(list = ls())
options(stringsAsFactors = F)
options(digits = 4)
load(file = "step1output.Rdata")
load(file = "step2output.Rdata")
exp1[1:4,1:4]
table(group_list1)
boxplot(exp1[1,]~group_list1)
bp=function(g){
library(ggpubr)
df=data.frame(gene=g,group=group_list1)
p <- ggboxplot(df, x = "group", y = "gene",
color = "group", palette = "jco",
add = "jitter")
p + stat_compare_means(label.y = 8)
}
bp(exp1[1,])
library(limma)
design=model.matrix(~factor(group_list1))
fit=lmFit(exp1,design)
fit=eBayes(fit)
deg1=topTable(fit,coef=2,number = Inf)
head(deg1)
deg1$logFC <- (deg1$logFC)
bp(exp1[rownames(deg1)[1],])
bp(exp1[rownames(deg1)[2],])
library(dplyr)
deg1 <- mutate(deg1,probe_id=rownames(deg1))
head(deg1)
library(rat2302.db)
ls("package:rat2302.db")
ids <- toTable(rat2302SYMBOL)
head(ids)
deg1 <- inner_join(deg1,ids,by="probe_id")
head(deg1)
logFC_t=1
change=ifelse(deg1$P.Value>0.01,'stable',
ifelse( deg1$logFC >logFC_t,'up',
ifelse( deg1$logFC < -logFC_t,'down','stable') )
)
deg1 <- mutate(deg1,change)
head(deg1)
table(deg1$change)
library(ggplot2)
library(clusterProfiler)
library(org.Hs.eg.db)
#s2e <- bitr(unique(deg1$symbol), fromType = "SYMBOL",
toType = c( "ENTREZID"),
OrgDb = org.Hs.eg.db)
#######################################################
head(s2e)
head(deg1)
deg1 <- inner_join(deg1,s2e,by=c("symbol"="SYMBOL"))
head(deg1)
save(group_list1,deg1,file = "step4output.Rdata")
Q:
s2e <- bitr(unique(deg1$symbol), fromType = "SYMBOL", toType = c( "ENTREZID"), OrgDb = org.Hs.eg.db)
和deg1 <- inner_join(deg1,s2e,by=c("symbol"="SYMBOL"))
得到的很尴尬,全部成了stable基因,如果用这个deg1
去画热图 就天下大乱了。
第五步:火山图
rm(list = ls())
options(stringsAsFactors = F)
load(file = "step4output.Rdata")
head(deg1)
plot(deg1$logFC,-log10(deg1$P.Value))
library(ggpubr)
table(dat1$change)
ggscatter(dat1, x = "logFC", y = "v",size=0.5,color = "change")
ggscatter(dat1, x = "logFC", y = "v", color = "change",size = 0.5,
label = "symbol", repel = T,
label.select = dat1$symbol[1:30] ,
#label.select = c('CD36','DUSP6'), #挑选一些基因在图中显示出来
palette = c("#00AFBB", "#999999", "#FC4E07") )
Result:
volcano
第五步:heatmap
rm(list = ls()) ## 魔幻操作,一键清空~
options(stringsAsFactors = F)
load(file = 'step1output.Rdata')
load(file = 'step2output.Rdata')
load(file = 'step4output.Rdata')
x=deg1$logFC
names(x)=deg1$probe_id
cg=c(names(head(sort(x),100)),
names(tail(sort(x),100)))
library(pheatmap)
pheatmap(exp1[cg,],show_colnames =F,show_rownames = F)
n=t(scale(t(exp1[cg,])))
thr=2
n[n>thr]=thr
n[n< -thr]= -thr
n[1:4,1:4]
pheatmap(n,show_colnames =F,show_rownames = F)
#显示分组
ac=data.frame(group=group_list1)
rownames(ac)=colnames(n)
pheatmap(n,show_colnames =T,
show_rownames = F,
annotation_col = ac,color = colorRampPalette(c("navy", "white", "firebrick3"))(50))
#保存
pdf(file = "heatmap.pdf")
pheatmap(n,show_colnames =F,
show_rownames = F,
#cluster_cols = F,
annotation_col=ac
)
dev.off()
Result:
heatmap
网友评论