数据来源：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数据挖掘（二）基因差异分析
网友评论