GEO数据挖掘（二）基因差异分析

数据来源：https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE54236

基因差异分析的目的是为了找到实验组与对照组表达有差异的基因，通过对这些基因进行一些文献调研或利用一些在线数据库的分析，可能会找到一个新的课题方向。

rm(list = ls())
options(stringsAsFactors = F)

load('GSE54236_download.Rdata')
#install.packages("limma")
library(limma)
design=model.matrix(~factor( group_list ))
fit=lmFit(dat,design)
fit=eBayes(fit)
options(digits = 4) #设置全局的数字有效位数为4
deg = topTable(fit,coef=2,adjust='BH', n=Inf) 
head(deg)
#         logFC AveExpr      t   P.Value adj.P.Val     B
#CAP2      1.515   7.905 11.123 9.928e-22 1.945e-17 38.65
#CDH13     1.599   6.741 10.154 4.540e-19 4.448e-15 32.70
#CENPF     2.128   8.419 10.086 6.955e-19 4.543e-15 32.28
#NEK2      2.231   4.454  9.969 1.446e-18 7.084e-15 31.57
#CENPA     2.438   3.895  9.922 2.142e-18 8.393e-15 31.19
#LOC344887 3.542   6.200  9.806 4.003e-18 1.307e-14 30.58
#设定上下调基因,以｜logFC｜=1为阈值，也就是基因表达量翻倍的情况
deg$g=ifelse(deg$P.Value>0.05,'stable', 
            ifelse( deg$logFC >1,'up', 
                    ifelse( deg$logFC < -1,'down','stable') ))
#统计上下调基因数量
table(deg$g)
#  down stable     up 
#   475  18843    278
save(deg,file = "GSE54236_deg.Rdata")
#保存为excle表
write.table(as.data.frame(deg), file="GSE54236_deg.xls", sep="\t", row.names=T)

火山图展示差异基因

#install.packages("ggplot2")
library(ggplot2)
p <- ggplot(
  # 数据、映射、颜色
  deg, aes(x = logFC, y = -log10(P.Value), colour=g)) +
  geom_point(alpha=0.4, size=3.5) +
  scale_color_manual(values=c("#c8e09f","#bdbdbf", "#ec5141"))+
  # 辅助线
  geom_vline(xintercept=c(-1,1),lty=4,col="black",lwd=0.8) +
  geom_hline(yintercept = -log10(0.05),lty=4,col="black",lwd=0.8) +
  # 坐标轴
  labs(x="log2(fold change)",
       y="-log10 (p-value)")+
  ggtitle("GSE54236")+
  theme_bw()+
  # 图例
  theme(plot.title = element_text(hjust = 0.5), 
        legend.position="right", 
        legend.title = element_blank())
p

热图

#热图（由于样本数量和基因都比较多，只取部分数据进行绘图）
#install.packages("pheatmap")
library(pheatmap)
deg = deg[order(deg$logFC),]
deg$name = rownames(deg)
#取|logFc|取大的10个上调基因和10个下调基因画图
up_gene <- tail(deg$name[which(deg$g == 'up')],10)
down_gene <- head(deg$name[which(deg$g == 'down')],10)
top10genes <- c(as.character(up_gene),as.character(down_gene))
diff=dat[top10genes,]
colnames(diff)
#取8个样本，4个TUMOR,4个NORMAL来画图示范
diff = diff[,c("GSM1310570","GSM1310571","GSM1310572","GSM1310573",
               "GSM1310726","GSM1310727","GSM1310728","GSM1310729")]
annotation_col=data.frame(rep(c("Tumor","Nomal"),each = 4))
rownames(annotation_col)=colnames(diff)
pheatmap(diff,
         annotation_col = annotation_col,
         scale = "row",
         show_rownames = F,
         show_colnames = F,
         fontsize = 10,
         fontsize_row = 3,
         fontsize_col = 3) 

画图这部分大家可以看看画图所使用包的说明书，或者去网上找一些画的好看的图的代码，再修改里面的一些参数就可以了，这部分内容主要是要拿到差异分析的结果。

3.png

参考

https://github.com/jmzeng1314/GEO

GEO数据挖掘

GEO数据挖掘（一）数据下载及基因ID转换

GEO数据挖掘（二）基因差异分析

GEO数据挖掘（三）使用DAVID数据库进行GO、KEGG富集分析

GEO数据挖掘（四）使用STRING数据库进行PPI分析