title: “GEO data analysis”
author: “yuluyang”
date: “2024-03-22”
生信技能树数据挖掘课程笔记~小洁老师授课
主要内容:
- 数据分组的内容
- 关键词的分组和多分组比较
- id map报错的原因及解决方法
- 基因组的注释流程
- 数据的行列互换及方差数值
- 画图示例代码
示例数据
library(GEOquery)
eSet = getGEO("GSE7305", destdir = '.', getGPL = F)
library(stringr)
# 紧接数据下载后
# 第一步:提取表达矩阵exp
eSet = eSet[[1]] # 要先打破列表的框架方便后续提取信息
exp <- exprs(eSet)
range(exp) # 看数据范围决定是否需要log,是否有负值,异常值
exp = log2(exp+1)
# 第二步:提取临床信息
pd <- pData(eSet)
## 要记得自己命名了哪些内容,'exp'是表达矩阵、'pd'是临床/分组信息
可以自己上GEO筛选一个芯片数据集,然后按以下流程复现一遍
只要命名方式不改变,后续作图代码可以直接copy然后run
提取后如何处理分组信息
if(F){ # 括号里是F,即不运行,这里只是示例有这种方法可供选择
# 第一种方法,有现成的可以用来分组的列
}else if(F){
# 第二种方法,眼睛数,自己生成
Group = rep(c("Disease","Normal"),each = 10) # 这个数据是很整齐排好的,前10个为`disease`,后10个为`normal`
# 如果两个组个数不一样的话咋整?
}else if(T){
# 第三种方法,使用字符串处理的函数获取分组
k = str_detect(pd$title,"Normal");table(k) # 选中信息这一列,找一个能区分组的关键词"Normal",然后`str_detect`分成`true`和`false`两组
## 记得检查`table`的结果
Group = ifelse(k,"Normal","Disease") # 对`true`、`false`结果进行分组分配
}
## 想起什么名字都可以
## 分组关系一定要对应清楚
# 需要把Group转换成因子,并设置参考水平,指定levels,对照组在前,处理组在后
Group = factor(Group,levels = c("Normal","Disease")) ## `factor()`转化为因子型向量,`levels=`指定因子的水平向量
Group
## `levels=`如果不指定的话,默认是按照字母顺序排列的
推荐使用第三种,优先考虑代码梳理
factor
因子型向量是R中一类特殊的向量类别
作用是给其他类别元素进行分组、计数,相当于统计学中分类变量,分可枚举【有序】和类别型【无序】两类
补充:如果两个组样本数量不一致
# 例如 对照组有10个,样本组15个
rep(c("control","disease"),times=c(10,15))
## 要记得检查生成的向量是不是先10个`control`再15个`disease`
rep(c("control","disease"),c(10,15))
## `times=`可以省略
## 不可以用`each=`
补充:如何检查分组设置与临床信息是否对应
# 新建一个数据框
data.frame(a=pd$title,
b=Group)
## 两列对照着看
补充:多分组里面取两个分组
# 举🌰,现编一个数据
pd$test=paste0(rep(c("a","b","c","d"),each=5),rep(c(1:5),times=4)) # 新增一列'test',分成a、b、c、d 4类,每一类5个
k1=str_detect(pd$test,"a");table(k1) # "a"为分组关键词
k2=str_detect(pd$test,"b");table(k2) # "b"为分组关键词
pd=pd[k1|k2,];pd # 再把这两列取子集出来
# 再上述取一遍
## k = str_detect(pd$,"");table(k)
## Group = ifelse(k,"","")
探针注释来源
- Bioconductor的注释包(通常可以用代码直接下载)
- GPL页面的表格文件解析(手动下载到工作目录然后按列读取)
- 官网下载对应产品的注释表(部分芯片是定制的,需要到对应平台官网查询,也有可能找不到)
- 自主注释(哒咩❌费劲)
用代码获取探针注释
# 推荐:用代码下载
library(tinyarray)
find_anno(gpl_number)
library(hgu133plus2.db);ids <- toTable(hgu133plus2SYMBOL)
## 可能会报错,提示你"use ids <- AnnoProbe::idmap('GPL570')...please try different 'type'"
## '?AnnoProbe::idmap'查看帮助文档,有哪些不同type
## 当你成功获取注释,`environment`会出现`ids`
## 检查`ids`列名是否为'probe_id'和'symbol',如不是,需手动修改
# 有些探针不对应具体的基因,需要去掉
## 去除空缺值
ids=na.omit(ids)
save(exp,Group,ids,file = "step2output.Rdata")
自己处理数据的时候,按分析的内容分段式以Rdata
存放好
另一种可能的报错
# 返回 no annotation available in Bioconductor and AnnoProbe
# 那么就需要在GEO网页上找到对应的GPL平台,打开该平台的表格找到ID、symbol等信息
# 下载文件再导入
遇到这种有难度的芯片就参考一下别人的代码
https://www.yuque.com/xiaojiewanglezenmofenshen/kzgwzl/sv262capcgg9o8s5
补充:其他获取探针注释的方法
# 方法1 BioconductorR包(最常用,已全部收入find_anno里面,不用看啦)
if(F){
gpl_number #看看编号是多少
#http://www.bio-info-trainee.com/1399.html #在这里搜索,找到对应的R包
library(hgu133plus2.db)
ls("package:hgu133plus2.db") #列出R包里都有啥
ids <- toTable(hgu133plus2SYMBOL) #把R包里的注释表格变成数据框
}
# 方法2 读取GPL网页的表格文件,按列取子集
##https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GPL570
# 方法3 官网下载注释文件并读取
# 方法4 自主注释,了解一下
#https://mp.weixin.qq.com/s/mrtjpN8yDKUdCSvSUuUwcA
PCA画图
dat=as.data.frame(t(exp))
library(FactoMineR)
library(factoextra)
dat.pca <- PCA(dat, graph = FALSE)
fviz_pca_ind(dat.pca,
geom.ind = "point", # show points only (nbut not "text")
col.ind = Group, # color by groups
palette = c("#00AFBB", "#E7B800"),
addEllipses = TRUE, # Concentration ellipses
legend.title = "Groups"
)
以上仅为作图代码示例
详细教程可参考STHDA网站,我后续也会进一步整理
http://www.sthda.com/english/articles/31-principal-component-methods-in-r-practical-guide/112-pca-principal-component-analysis-essentials/
heatmap 画图
g = names(tail(sort(apply(exp,1,sd)),1000)) # 按方差筛选差异最大的1000个基因
n = exp[g,]
library(pheatmap)
annotation_col = data.frame(row.names = colnames(n),
Group = Group)
pheatmap(n,
show_colnames =F,
show_rownames = F,
annotation_col=annotation_col,
scale = "row", #按行标准化,只保留行内差别,不保留行间差别,会把数据范围缩放到大概-5~5之间
breaks = seq(-3,3,length.out = 100) #设置色带分布范围为-3~3之间,超出此范围的数字显示极限颜色
)
差异分析——画火山图
library(limma)
design = model.matrix(~Group)
fit = lmFit(exp,design)
fit = eBayes(fit)
deg = topTable(fit,coef = 2,number = Inf)
#为deg数据框添加几列
#1.加probe_id列,把行名变成一列
library(dplyr)
deg = mutate(deg,probe_id = rownames(deg))
#2.加上探针注释
ids = distinct(ids,symbol,.keep_all = T)
#其他去重方式在zz.去重方式.R
deg = inner_join(deg,ids,by="probe_id")
nrow(deg) #如果行数为0就是你找的探针注释是错的
#3.加change列,标记上下调基因
logFC_t = 2
p_t = 0.005
#思考,如何使用padj而非p值
k1 = (deg$P.Value < p_t)&(deg$logFC < -logFC_t)
k2 = (deg$P.Value < p_t)&(deg$logFC > logFC_t)
deg = mutate(deg,change = ifelse(k1,"down",ifelse(k2,"up","stable")))
table(deg$change)
## 画火山图
library(ggplot2)
ggplot(data = deg, aes(x = logFC, y = -log10(P.Value))) +
geom_point(alpha=0.4, size=3.5, aes(color=change)) +
scale_color_manual(values=c("blue", "grey","red"))+
geom_vline(xintercept=c(-logFC_t,logFC_t),lty=4,col="black",linewidth=0.8) +
geom_hline(yintercept = -log10(p_t),lty=4,col="black",linewidth=0.8) +
theme_bw()
差异分析——画差异基因热图
# 表达矩阵行名替换为基因名
exp = exp[deg$probe_id,]
rownames(exp) = deg$symbol
diff_gene = deg$symbol[deg$change !="stable"]
n = exp[diff_gene,]
library(pheatmap)
annotation_col = data.frame(group = Group)
rownames(annotation_col) = colnames(n)
pheatmap(n,show_colnames =F,
show_rownames = F,
scale = "row",
#cluster_cols = F,
annotation_col=annotation_col,
breaks = seq(-3,3,length.out = 100)
)
富集分析
# 加ENTREZID列,用于富集分析(symbol转entrezid,然后inner_join)
library(clusterProfiler)
library(org.Hs.eg.db)
s2e = bitr(deg$symbol,
fromType = "SYMBOL",
toType = "ENTREZID",
OrgDb = org.Hs.eg.db)#人类,注意物种
# 一部分基因没匹配上是正常的。<30%的失败都没事。
# 其他物种http://bioconductor.org/packages/release/BiocViews.html#___OrgDb
nrow(deg)
deg = inner_join(deg,s2e,by=c("symbol"="SYMBOL"))
# 多了几行少了几行都正常,SYMBOL与ENTREZID不是一对一的。
nrow(deg)
save(exp,Group,deg,logFC_t,p_t,file = "step4output.Rdata")
library(ggthemes)
library(org.Hs.eg.db)
library(dplyr)
library(enrichplot)
# 第一步:输入数据
gene_diff = deg$ENTREZID[deg$change != "stable"]
# 第二步:富集
ekk <- enrichKEGG(gene = gene_diff,organism = 'hsa')
ekk <- setReadable(ekk,OrgDb = org.Hs.eg.db,keyType = "ENTREZID")
ego <- enrichGO(gene = gene_diff,OrgDb= org.Hs.eg.db,
ont = "ALL",readable = TRUE)
# setReadable和readable = TRUE都是把富集结果表格里的基因名称转为symbol
class(ekk)
# 第三步:可视化
barplot(ego, split = "ONTOLOGY") +
facet_grid(ONTOLOGY ~ ., space = "free_y",scales = "free_y")
barplot(ekk)
#或者是dotplot
— 更多资料
GSEA:https://www.yuque.com/docs/share/a67a180f-dd2b-4f6f-96c2-68a4b86fe862?#
Y叔的书:http://yulab-smu.top/clusterProfiler-book/index.html
GOplot:https://mp.weixin.qq.com/s/LonwdDhDn8iFUfxqSJ2Wew
网上的资料和宝藏无穷无尽,学好R语言慢慢发掘~