富集分析15-1/2
芯片注释：探针与基因的对应关系
注释来源：
-
GEO数据库中GPL页面的表格soft文件解析
-
Bioconductor 的注释包
用R获取芯片探针与基因的对应关系三部曲 -
官网下载对应产品的注视表格
-
自主注释
芯片探针序列的基因组注释
image
富集分析
输入数据：差异基因的entrezid; 所有基因的entrezid
id转换：bitr()
symbol：常用基因名（PD-L1/CD36等
ENTREZID:富集分析指定id
之间的相互转换
富集分析：
1.输入数据：差异基因及所有基因的ENTREZID
2.基因SYMBOL和ENTREZID之间的转换：
所需转换基因SYMBOL组成向量
告诉函数给出的是什么类型的id
输出数据是什么类型的id
转换数据依据的数据库
*完成数据转换后基因数量发生改变：本来两者并非一一对应的关系，属于正常现象
KEGG:系统分析基因功能，基因组信息的数据库，“理解生物系统的高级功能和实用程序资源库“
GO:
细胞组分：细胞每个部分和细胞外环境
分子功能：分子水平的活性，ex：催化和结合活性
生物过程：一个或多个分子功能有序结合而产生的系列事件
通过差异基因做GO富集分析，将基因按照不同功能归类，对基因进行注释和分类
image
GeneRatio ：我所输入基因数据对应该通路的个数/我所输入的基因数据匹配到数据库中的基因个数
BgRatio：在数据库中对应通路的基因个数/数据库中所有通路的基因个数
富集分析目的：这条通路的gene是否在差异基因中足够多，BgRatio随机抽样抽中GeneRatio的比例
image
image
下面的图用enrichplot画
image
image
image
image
image
代码分析流程
image
1.数据下载,提取矩阵和临床信息：
rm(list = ls()) #清除环境里的所有东西
options(stringsAsFactors = F) #避免因子的影响
library(GEOquery)
gse_number = "GSE56649" #指定后面编号的变量名，后续方便改动
下载数据：
eSet <- getGEO(gse_number,
destdir = '.', #下载文件从哪里获取，“·”：当前目录
getGPL = F) #需不需要下载gpl注释文件
class(eSet) #数据类型
length(eSet) #列表里元素数量
eSet = eSet[[1]]
#(1)提取表达矩阵exp
exp <- exprs(eSet) #获得表达矩阵
exp[1:4,1:4] #大致观察数据是否取过log（一般不超过20），或者range（exp）得到函数取值范围
boxplot(exp) #或者用箱线图也能看数据大致范围
exp = log2(exp+1) #取log，括号里为啥要➕1，避免矩阵出现0而影响数据，因为其他数据+1取log的影响不大
boxplot(exp)
#(2)提取临床信息
pd <- pData(eSet) #行名是不同患者，列名是项目
#(3)调整pd的行名顺序与exp列名完全一致
p = identical(rownames(pd),colnames(exp));p #确定是否一致
if(!p) exp = exp[,match(rownames(pd),colnames(exp))] #不一致的话按照exp的列名排序
#(4)提取芯片平台编号
gpl_number <- eSet@annotation
gpl_number
save(gse_number,pd,exp,gpl_number,file = "step1output.Rdata")
贴一个Jimmy老师的github GEO数据挖掘代码
https://github.com/jmzeng1314/GEO
2.分组信息
# Group(实验分组)和ids(探针注释)
rm(list = ls())
load(file = "step1output.Rdata")
library(stringr)
# 1.Group----
# 第一类，有现成的可以用来分组的列
if(F) Group = pd$`disease state:ch1`
#第二类，自己生成
if(F){
Group=c(rep("RA",times=13),
rep("control",times=9))
}
rep(c("RA","control"),times = c(13,9)) #两种方式结果一样
>[1] "RA" "RA" "RA" "RA" "RA" "RA" "RA" "RA" "RA"
[10] "RA" "RA" "RA" "RA" "control" "control" "control" "control" "control"
[19] "control" "control" "control" "control"
#第三类，匹配关键词，自行分类
Group=ifelse(str_detect(pd$source_name_ch1,"control"),"control","RA") #是否含有“control”，是输出“control”，否“ra“
#设置参考水平，指定levels，对照组在前，处理组在后
Group = factor(Group, #生成/转换因子的函数
levels = c("control","RA")) #指定因子前后顺序，确保对照组在前
> Group
[1] RA RA RA RA RA RA RA RA RA RA RA
[12] RA RA control control control control control control control control control
Levels: control RA #因子里面有那几个取值，确保生成的因子，control组在第一个位置，treat组在level第二个位置，如果多组的话依次排开
因子的第一个位置相当于整组的一个管家（参考），所以需要放到level的第一个位置。所以需要自己指定
注意levels与因子内容必须对应一致
Group = pd$`disease state:ch1`
Group = factor(Group,
levels = c("healthy control","rheumatoid arthritis"))
#如果数据列表中，没有很好的分组信息，需要加以处理：
pd$characteristics_ch1
[1] "disease state: rheumatoid arthritis" "disease state: rheumatoid arthritis"
[3] "disease state: rheumatoid arthritis" "disease state: rheumatoid arthritis"...
#方法1
str_remove(pd$characteristics_ch1,"disease state: ")
[1] "rheumatoid arthritis" "rheumatoid arthritis" "rheumatoid arthritis"
[4] "rheumatoid arthritis" "rheumatoid arthritis" "rheumatoid arthritis"
#方法2
str_split(pd$characteristics_ch1,": ")
[[1]]
[1] "disease state" "rheumatoid arthritis"
[[2]]
[1] "disease state" "rheumatoid arthritis"
此处生成的是列表，需交代继续生成矩阵
str_split(pd$characteristics_ch1,": ",simplify=T)
[,1] [,2]
[1,] "disease state" "rheumatoid arthritis"
[2,] "disease state" "rheumatoid arthritis"
[3,] "disease state" "rheumatoid arthritis"
[4,] "disease state" "rheumatoid arthritis"
再截取矩阵的第二列
str_split(pd$characteristics_ch1,": ",simplify=T)[,2]
[1] "rheumatoid arthritis" "rheumatoid arthritis" "rheumatoid arthritis"
[4] "rheumatoid arthritis" "rheumatoid arthritis" "rheumatoid arthritis"
3. 探针注释
ids-----------------
方法1 BioconductorR包(最常用)
gpl_number
[1] "GPL570"
去网站找到芯片所需要的R包 ，hgu133plus2（包里含有探针和基因的关系） ，加上“.db”就是所需的r包
http://www.bio-info-trainee.com/1399.html
image.png
安装包并加载
if(!require(hgu133plus2.db))BiocManager::install("hgu133plus2.db")
library(hgu133plus2.db)
ls("package:hgu133plus2.db") #看包里面的全部内容（注释包）
image.png
找到对应关系:hgu133plus2SYMBOL,并用totable处理获取对应关系里面的内容
ids <- toTable(hgu133plus2SYMBOL)
head(ids)
image.png
方法2 读取GPL平台的soft文件，按列取子集https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GPL570
if(F){
注：soft文件列名不统一，活学活用，有的GPL平台没有提供注释，如GPL16956
a = getGEO(gpl_number,destdir = ".")
b = a@dataTable@table
colnames(b)
ids2 = b[,c("ID","Gene Symbol")]
colnames(ids2) = c("probe_id","symbol")
ids2 = ids2[ids2$symbol!="" & !str_detect(ids2$symbol,"///"),]
}
#将矩阵中对应多个基因的探针，和无基因的探针去掉
方法3 官网下载,文件读取
http://www.affymetrix.com/support/technical/byproduct.affx?product=hg-u133-plus
方法4 自主注释
https://mp.weixin.qq.com/s/mrtjpN8yDKUdCSvSUuUwcA
save(exp,Group,ids,gse_number,file = "step2output.Rdata")
翻一下这篇文章前面提到的四种方法
用R获取芯片探针与基因的对应关系三部曲
果子老师的经验分享
一些 tips
- ls() 查看环境中存在的变量
- ls("package:tidyr") 查看R包中存在的变量
- toTable() 得到探针 ID 的注释矩阵
