GEO实战演练（一）

《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