芯片注释：探针与基因的对应关系
注释来源：
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富集分析，将基因按照不同功能归类，对基因进行注释和分类
GeneRatio ：我所输入基因数据对应该通路的个数/我所输入的基因数据匹配到数据库中的基因个数
BgRatio：在数据库中对应通路的基因个数/数据库中所有通路的基因个数
富集分析目的：这条通路的gene是否在差异基因中足够多，BgRatio随机抽样抽中GeneRatio的比例
下面的图用enrichplot画
代码分析流程
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") #看包里面的全部内容（注释包）
找到对应关系:hgu133plus2SYMBOL,并用totable处理获取对应关系里面的内容
ids <- toTable(hgu133plus2SYMBOL)
head(ids)
方法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 的注释矩阵
