关于ID转换
2019-08-17 本文已影响0人
BioJournal_Link
背景
不同的数据库，标识同一个分子所用的的”代码“可能是不一样的，因为不同的数据库有自己标识分子的规则。当我们对一种ID不熟悉时，就经常需要转换成我们熟悉的ID。比如，在芯片数据里面，我们最常用的就是PROBE_ID与GENE_SYMBOL的转换，即探针基因名转换。通过转换，我们就可以知道差异表达的基因是哪些基因，而不是那些我们完全不熟悉的探针名。
芯片数据下载并提取表达矩阵后，我们得到的是探针表达矩阵，即变量（列）是样本GSM，观测（行）是探针PROBE。要想把PROBE转换成其他ID，我们就需要有一个PROBE_ID与其他ID有对应关系的文件。
前人把有这个对应关系的文件写成了R包（根据功能暂且叫他注释包），这样我们只要加载这个注释包，我们就可以获取ID之间的“关系”文件。因为不同的PLATFORM所对应的PROBE_ID是不一样的，所以不同的PLATFORM所对应的R包不一样。所以，又有人总结了常见的PLATFORM与注释包对应关系的文件。只要我们确定了PLATFORM，就可以从文件中找出相应的R包。
PLATFORM与注释包文件——生信技能树——传送门
PLATFORM与注释包文件——我的简书——修改后的传送门
PLATFORM与注释包platformMap——果子学生信——传送门
当有专门的注释包时，有以下3种方法可以得到【ID转换文件】
安装注释R包
#因为是根据注释R包得到ID转换文件，所以必须先安装注释R包
#常用的注释R包在【PLATFORM与注释包文件——生信技能树】和【PLATFORM与注释包platformMap——果子学生信】有很好的总结
#安装【生信技能树】里面的R包
#【PLATFORM与注释包文件——生信技能树】这个链接里的文本保存并读取后，按里面的代码读取后并不是目标格式，所以自己改进了文本与代码，见【PLATFORM与注释包文件——我的简书】
#打开【PLATFORM与注释包文件——我的简书】链接，复制文本于目标文件夹，粘贴并保存名为GPL_info.txt的文本文件
GPL_info <- read.table("GPL_info.txt",header = T)#读取该文件
#设置镜像
options("repos"="https://mirrors.ustc.edu.cn/CRAN/")
options(BioC_mirror="https://mirrors.ustc.edu.cn/bioc/")
#安装没有安装过的注释包
for (i in 1:nrow(GPL_info)){
print(i)
platformDB=GPL_info$R_package[i]
if( !(platformDB %in% rownames(installed.packages()))) {
BiocManager::install(platformDB,update = F,ask = F)
}
}
#安装【果子学生信】里面的R包
#打开【PLATFORM与注释包platformMap——果子学生信】链接
#于目标文件夹中创建platformMap.txt文本文档
#复制下面文本于platformMap.txt文本文档中，保存为UTF-8编码格式
plateformMap <- data.table::fread("platformMap.txt")
platformMap <- platformMap[,-1]
for (i in 1:nrow(platformMap)){
print(i)
platformDB=platformMap$R_package[i]
if( !(platformDB %in% rownames(installed.packages()))) {
BiocManager::install(platformDB,update = F,ask = F)
}
}
根据R包制作ID转换文件
#R包方法一
index = lareset@annotation#获取ExpressionSet的平台信息，即GPL
platformDB <- GPL_info$R_package[grep(index,GPL_info$gpl)]#根据GPL找出该注释包名称
#platformDB <- platformMap$R_package[grep(index,platformMap$gpl)]
platformDB <- as.character(platformDB)#因子转换成字符
#加载或安装该注释包
if(!require(platformDB,character.only = T)) BiocManager::install(platformDB,update = F,ask = F)
PROBE_ID <- rownames(express)#获取表达矩阵的行名，就是探针名称
#再根据PROBE_ID找出该注释包中对应的SYMBOL_ID
SYMBOL_ID <- annotate::lookUp(PROBE_ID,platformDB, "SYMBOL")
#SYMBOL_ID是个list，其包含了PROBE_ID数目的子list，每个子list以每个PROBE_ID命名，每个子list的内容就是对应的SYMBOL_ID。
SYMBOL_ID <- unlist(SYMBOL_ID)#提取SYMBOL_ID，转换成字符向量
probe2symbol_1 <- data.frame(PROBE_ID,SYMBOL_ID,stringsAsFactors = F)#制作probe2symbol转换文件
dim(probe2symbol_1)#33297行，2列
#R包方法二：
index = lareset@annotation#获取ExpressionSet的平台信息，即GPL
platform <- GPL_info$bioc_package[grep(index,GPL_info$gpl)]
#platform <- platformMap$bioc_package[grep(index,platformMap$gpl)]
probe2symbol_2 <- toTable(get(paste0(platform,"SYMBOL")))#一步到位
colnames(probe2symbol_2) <- c("PROBE_ID","SYMBOL_ID")
dim(probe2symbol_2)#19827行，2列
方法一与方法二的区别
#那么上述两种方法得出的【ID转换文件】有什么区别
#首先，我们发现行数不一样，probe2symbol_1有33297行，probe2symbol_2有19827行
#因为都是从同一个注释包里提取出来的，理论上应该是一样的
length(unique(probe2symbol_1$PROBE_ID))
#33297,说明probe2symbol_1里的PROBE_ID都是唯一值
length(unique(probe2symbol_1$SYMBOL_ID))
#18835<33297，说明SYMBOL_ID不全是唯一值，里面可能有写symbol是重复的或者是NA
table(is.na(probe2symbol_1$SYMBOL_ID))#看下有多少个缺失值
#FALSE TRUE
#19827 13470
#有13470个缺失值，有19827个SYMBOL_ID，刚好与probe2symbol_2的行数对应
#也就是说，probe2symbol_2很可能是toTable函数直接删除了SYMBOL_ID为缺失值后的结果
probe2symbol_1_na <- probe2symbol_1[!is.na(probe2symbol_1$SYMBOL_ID),]
probe2symbol_2_na <- probe2symbol_2[!is.na(probe2symbol_2$SYMBOL_ID),]
#去除掉含有缺失值symbol的行
table(probe2symbol_1_na$SYMBOL_ID %in% probe2symbol_2_na$SYMBOL_ID)# TRUE 19827
table(probe2symbol_2_na$SYMBOL_ID %in% probe2symbol_1_na$SYMBOL_ID)# TRUE 19827
length(unique(probe2symbol_1_na$SYMBOL_ID))#18834
length(unique(probe2symbol_2_na$SYMBOL_ID))#18834
#结论：这两个ID转换文件本质上是一模一样的，而显然方法二的代码比较简洁
【问题】为什么probe2symboll_1里面会有probe对应无symbol的情况，那该probe的存在有何意义
【猜想】平台套餐不一样，用同样的平台，贵一点的套餐可以给你测更多的probe，这样对结果可能更准确
#R包方法三：
index = lareset@annotation#获取ExpressionSet的平台信息，即GPL
platformDB <- GPL_info$R_package[grep(index,GPL_info$gpl)]#根据GPL找出该注释包名称
#platformDB <- platformMap$R_package[grep(index,platformMap$gpl)]
platformDB <- as.character(platformDB)#因子转换成字符
#加载或安装该注释包
if(!require(platformDB,character.only = T)) BiocManager::install(platformDB,update = F,ask = F)
ls(paste0("package:",platformDB))
probe2symbol_3 <- eval(parse(text=paste0('select(',platformDB,',keys = keys(',platformDB,'),columns=','"SYMBOL")')))
#为了代码通用，这里把字符串当做命令运行
#probe2symbol_ENTREID <- select(hugene10sttranscriptcluster.db,keys = keys(hugene10sttranscriptcluster.db),columns = c('SYMBOL','ENTREZID'))
#select(注释R包，keys（源数据），columns（目标数据）)，即根据【源数据】提取出【目标数据】
#keys = keys(注释R包)，即在该注释R包提取PROBE_ID，默认PROBE_ID
#keys = keys(注释R包，keytype=c("SYMBOL","ENTREZID"))，你也可以提取出SYMBOL_ID等其他ID
colnames(probe2symbol_3) <- c("PROBE_ID","SYMBOL_ID")
dim(probe2symbol_3)#40284行
方法三与前两种方法的区别
length(unique(probe2symbol_3$PROBE_ID))#33297 PROBE_ID竟然不是唯一值
table(is.na(probe2symbol_3$PROBE_ID))#FALSE 40284 说明没有缺失值
#结论：probe2symbol的PROBE_ID多了几千个重复项
#【问题】但为什么probe2symbol的PROBE_ID会有重复项呢？
#【猜想】会不会是一个PROBE_ID对应多个SYMBOL_ID？
table(is.na(probe2symbol$SYMBOL_ID))#FALSE 29122 TRUE 11162 ,有11162个NA
length(unique(probe2symbol$SYMBOL_ID))#21312，包括NA
probe2symbol_na <- probe2symbol_3[!is.na(probe2symbol_3$SYMBOL_ID),]#去掉含有缺失值的行
dim(probe2symbol)#29122 2
length(unique(probe2symbol$PROBE_ID))#22135个唯一值，其他6987（29122-22135）个是重复项
length(unique(probe2symbol$SYMBOL_ID))#21311，不包括NA，其他7811（29122-21311）个是重复项
#假如一个PROBE_ID对应多个SYMBOL_ID
#那么对应多个SYMBOL_ID的PROBE_ID在probe2symbol$SYMBOL_ID里必定是重复项
#先挑选出有重复项的PROBE_ID,因为只有有重复项的PROBE_ID才有可能一对多
probe_rep <- names(table(probe2symbol$PROBE_ID)>1)[table(probe2symbol$PROBE_ID)>1]
length(probe_rep)#2210个探针拥有重复项
probe_one <- names(table(probe2symbol$PROBE_ID)==1)[table(probe2symbol$PROBE_ID)==1]
length(probe_one)#19925,19925+2210=22135
#(29122-19925)/2210=4.16,2210个重复项平均每个重复4个左右
#写一个循环，挑选出probe2symbol里PROBE_ID一对多PROBE_ID
#然后看SYMBOL_ID的种类，如果种类大于1，那么说明该PROBE_ID属于一对多的情况
plus <- c()
for (i in 1:length(probe_rep)) {
rep <- probe2symbol$PROBE_ID %in% probe_rep[i]
rep_probe <- probe2symbol[rep,]
plus[i] <- length(table(rep_probe[,2]))>1
}
length(plus)
table(plus)#说明这2210个重复项全部都是一对多
#验证一下
i <- 1
probe_1 <- probe_rep[i]
probe_2 <- probe2symbol[probe2symbol$PROBE_ID %in% probe_1,]
table(probe_2[,2])
# i=1时（基因名差不多）
# OR4F17 OR4F4 OR4F5
# 1 1 1
# i=35时（基因名差不多）
# ZBTB8A ZBTB8B
# 1 1
#我们可以发现，这些一对多的PROBE_ID所对应的SYMBOL_ID，很多名字是“相似”的
# i=66时（基因名差不多）
# NBPF1 NBPF10 NBPF11 NBPF12 NBPF14 NBPF15
# 1 1 1 1 1 1
# NBPF19 NBPF20 NBPF25P NBPF26 NBPF3 NBPF8 NBPF9
# 1 1 1 1 1 1 1
# i=190时（基因名有区别）
# LOC100996724 LOC653513 PDE4DIP PFN1P2
# 1 1 1 1
#曾老师的解释：因为多个基因位置会重叠，且一个基因会有多个名字
#所以会出现一个PROBE_ID会对应多个SYMBOL_ID的情况
当没有专门的注释R包时
#方法一
library(GEOquery)
lareset <- gset[[1]]
lareset@annotation
GPLXXX <- getGEO(lareset@annotation,destdir = getwd())
#目标文件夹将会得到GPLXXX.soft的文件
#小插曲，因为是从外国网站下载，所以下载速度会特别慢，如果网速不好，建议去网吧。。。
#另外，下载时会出现下载不完全的情况，但他依旧会下载，但会给你warning
#如果下载完全，它会出现 |======================| 100%
GPL_anno <- eval(parse(text = paste0("Table(",lareset@annotation,")")))
#GPL_anno <- Table(GPLXXX)
#GPL_anno <- GPLXXX@dataTable@table
#GPLXXX@dataTable@columns#可以了解GPL_anno各个变量的含义
#Columns(GPL6244)[,]#可以了解GPL_anno各个变量的含义
#GPL6244@header#可以了解GPL信息
GPL_anno <- GPL_anno[,c("ID","gene_assignment")]
GPL_anno <- GPL_anno[(GPL_anno$gene_assignment != "---"),]
View(GPL_anno)#发现我们要的PROBE_ID和SYMBOL_ID在ID列和gene_assigment列
GPL_anno <- GPL_anno[,c("ID","gene_assignment")]#提取这两列
GPL_anno <- GPL_anno[(GPL_anno$gene_assignment != "---"),]#删除含有“---”的行
#仔细观察gene_assignment，发现每个名词都以“//”分隔，我们要的SYMBOL_ID就在第二个
SYMBOL <- c()
SYMBOL_ID <- c()
for (i in 1:nrow(GPL_anno)) {
SYMBOL[i] <- strsplit(GPL_anno$gene_assignment[i],split = "//")[[1]][2]
SYMBOL_ID[i] <- strsplit(SYMBOL[i],split = " ")[[1]][2]
}
length(SYMBOL)
length(SYMBOL_ID)
PROBE_ID <- GPL_anno$ID
probe2symbol_4 <- data.frame(PROBE_ID,SYMBOL_ID)
#以下是更简洁的代码，取自【果子学生信】
library(dplyr)
library(tidyr)
probe2symbol_4 <- GPL_anno %>%
select(ID,gene_assignment) %>%
filter(gene_assignment != "---") %>%
separate(gene_assignment,c("drop","symbol"),sep="//") %>%
select(-drop)
#方法二
#网站下载soft格式文件，读取
GPL_anno <-data.table::fread("./GSEXXX_family.soft/GSEXXX_family.soft",skip ="ID")
#其余同前，本质上是一个方法
library(dplyr)
library(tidyr)
probe2symbol_4 <- GPL_anno %>%
select(ID,gene_assignment) %>%
filter(gene_assignment != "---") %>%
separate(gene_assignment,c("drop","symbol"),sep="//") %>%
select(-drop)
