GEO数据挖掘4
sunqi
2020/7/12
概述
对GEO数据进行差异分析
简单比较
代码语言：javascript复制rm(list = ls())
options(stringsAsFactors = F)
options(digits = 4) #设置全局的数字有效位数为4
load(file = 'step1-output.Rdata')
#导入的数据中，dat为表达矩阵，group_list为分组信息
#按照group_list分组画箱线图
boxplot(dat[1,]~group_list)
# 定义函数，用于和绘制箱式图
bp=function(g){
# 高级绘图包，用于绘制自定义图
# 比如添加P值之类的操作
library(ggpubr)
# 准备需要的数据
df=data.frame(gene=g,stage=group_list)
p <- ggboxplot(df, x = "stage", y = "gene",
color = "stage", palette = "jco",
add = "jitter")
# 添加p值
p stat_compare_means()# 添加比较函数，默认使用wilcox
}
# 对4个样本的箱式图绘制
bp(dat[1,])
bp(dat[2,])
bp(dat[3,])
bp(dat[4,])
基因差异分析
这里需要使用差异比较用到的limma包，在使用这个包进行分析之前，需要准备三个矩阵 * 表达矩阵 * 分组矩阵 * 差异比较矩阵
代码语言：javascript复制# 如果没有这个包，需要进行install，应该是biocmanger下载
library(limma)
# 数据预处理
# 分组矩阵准备
design <- model.matrix(~0 factor(group_list))
# 给分组矩阵添加分组名
colnames(design)=levels(factor(group_list))
head(design)## Control Vemurafenib
## 1 1 0
## 2 1 0
## 3 1 0
## 4 0 1
## 5 0 1
## 6 0 1# 表达矩阵
exprSet=dat
# 给分组矩阵添加样本名
rownames(design)=colnames(exprSet)
design## Control Vemurafenib
## GSM1052615 1 0
## GSM1052616 1 0
## GSM1052617 1 0
## GSM1052618 0 1
## GSM1052619 0 1
## GSM1052620 0 1
## attr(,"assign")
## [1] 1 1
## attr(,"contrasts")
## attr(,"contrasts")$`factor(group_list)`
## [1] "contr.treatment"# 差异比较矩阵，也就是哪个和那个比的矩阵
contrast.matrix<-makeContrasts("Vemurafenib-Control",
levels = design)
# 进行差异比较的主要函数
# 这里需要输入三个文件，表达矩阵，分组矩阵，比较矩阵
deg = function(exprSet,design,contrast.matrix){
##线性模型拟合
fit <- lmFit(exprSet,design)
##添加比较矩阵，进行差值计算
fit2 <- contrasts.fit(fit, contrast.matrix)
## 贝叶斯检验
fit2 <- eBayes(fit2)
## 生成检验结果报告
tempOutput = topTable(fit2, coef=1, n=Inf)
## 去除缺失的数据
nrDEG = na.omit(tempOutput)
return(nrDEG)
}
# 调用函数
deg = deg(exprSet,design,contrast.matrix)
# 保存比较结果后续使用
save(deg,file = 'deg.Rdata')logFC为输入的表达矩阵中case一组的平均表达量减去control一组的平均表达量的值
火山图绘制
进行差异比较滞后，绘制火山图，对差异基因进行可视化
代码语言：javascript复制# 这里的if函数用于折叠代码，方便阅读，没什么重要意义
nrDEG=deg# 重新赋值
attach(nrDEG)# 方便写代码
plot(logFC,-log10(P.Value))# 对p取对数
# 上面的图不太美观，进行美化
# 有点饿了，想吃烤串，一会儿学习
library(ggpubr)
df=nrDEG
df$v= -log10(P.Value) #对p进行取负对数
# 散点图
ggscatter(df, x = "logFC", y = "v",size=0.5)
# 火山图倒置之后，继续对基因标注
# 对上调和下调的基因进行标注
df$g=ifelse(df$P.Value>0.01,'stable', #判断是否为稳定基因
ifelse( df$logFC >2,'up', #判断是否上调
ifelse( df$logFC < -2,'down','stable') ))#判断是否下调
# 查看有多少上调和下调基因
table(df$g)##
## down stable up
## 51 18686 84# 添加name列
df$name=rownames(df)
#散点图
ggscatter(df, x = "logFC", y = "v",size=0.5,color = 'g')
# 添加分组信息
ggscatter(df, x = "logFC", y = "v", color = "g",size = 0.5,
label = "name", repel = T,
label.select = head(rownames(deg)), #挑选一些基因在图中显示出来
palette = c("#00AFBB", "#E7B800", "#FC4E07") )
# 保存图片
# ggsave('volcano.png')
# MA图绘制
# MA图用于展示两组变化关系，(x y)/2作为横轴，(y-x)作为纵轴
# 方便观察两组数值偏离程度（差异程度）
ggscatter(df, x = "AveExpr", y = "logFC",size = 0.2)
df$p_c = ifelse(df$P.Value<0.001,'p<0.001',
ifelse(df$P.Value<0.01,'0.001<p<0.01','p>0.01'))
table(df$p_c )##
## 0.001<p<0.01 p<0.001 p>0.01
## 2213 5291 11317ggscatter(df,x = "AveExpr", y = "logFC", color = "p_c",size=0.2,
palette = c("green", "red", "black") )
# 保存结果
# ggsave('MA.png')差异基因热图
代码语言：javascript复制load(file = 'step1-output.Rdata')
x=deg$logFC # deg为差异比较的结果
names(x)=rownames(deg) # 将基因名添加
# 对差异排序，并取前后100的基因名，赋值到cg
cg=c(names(head(sort(x),100)),
names(tail(sort(x),100)))
# 绘制热图
library(pheatmap)
# 对差异基因绘制热图
pheatmap(dat[cg,],show_colnames =F,show_rownames = F)
# 对数据归一化后绘制热图
n=t(scale(t(dat[cg,])))
n[n>2]=2
n[n< -2]= -2
pheatmap(n,show_colnames =F,show_rownames = F)
# 添加分组信息
ac=data.frame(g=group_list)
rownames(ac)=colnames(n)
pheatmap(n,show_colnames =F,
show_rownames = F,
cluster_cols = F,
annotation_col=ac)
结束语
这里对GEO数据的差异分析已经结束，后续为kegg和go分析
love&peace
