前言
看了这么多期的 circos
,都有点乏了,换换口味也是好的。
下面,咱就开始吧。
做生信分析,总是免不了要给基因 ID
和 Symbol
转换来转换去。
方法
一般要进行 ID
和 Symbol
的转换呢,主要有两种方式:
- 网站提供的工具,比如
biodbnet
- 编写代码
1. 用网站转换
如果不会编写代码的话,可以使用这个网站 biodbnet
这种方式比较简单,比如上面的例子,我们输入的是人类(9606
)基因 symbol
,需要对应的基因 id
,提交之后
可以下载转换的结果。
但是以我的经验来说,这个网站如果输入的基因很多,速度非常慢,而且很多基因 symbol
无法转换到 id
的,所以对于有编程基础的朋友,并不推荐这种方式
2. 编程实现
编程的话,很多语言都可以实现,看自己比较喜欢,比较擅长用什么语言
下面主要介绍一下 R
以及 Python
两种语言实现方式
2.1 R
用 R
实现的话,一般都是使用 org.Hs.eg.db
这个模块提供的数据来进行转换
安装和导入
# install
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("org.Hs.eg.db")
# library
library('org.Hs.eg.db')
安装这个包之后,对于的路径下面会有一个 org.Hs.eg.sqlite
文件,存储了人类基因数据,后面的各种转换其实都是对这个文件进行操作。
查看基本信息
# 获取所有可用的表
columns(org.Hs.eg.db)
# [1] "ACCNUM" "ALIAS" "ENSEMBL" "ENSEMBLPROT" "ENSEMBLTRANS" "ENTREZID"
# [7] "ENZYME" "EVIDENCE" "EVIDENCEALL" "GENENAME" "GO" "GOALL"
# [13] "IPI" "MAP" "OMIM" "ONTOLOGY" "ONTOLOGYALL" "PATH"
# [19] "PFAM" "PMID" "PROSITE" "REFSEQ" "SYMBOL" "UCSCKG"
# [25] "UNIGENE" "UNIPROT"
从上面的输出信息可以看出,包含了很多数据表,如 ENSEMBL、ENTREZID、SYMBOL
等
# keytype 配合 keys 使用,在 select 函数中匹配 keys 参数指定的 id
keytypes(org.Hs.eg.db)
# [1] "ACCNUM" "ALIAS" "ENSEMBL" "ENSEMBLPROT" "ENSEMBLTRANS" "ENTREZID"
# [7] "ENZYME" "EVIDENCE" "EVIDENCEALL" "GENENAME" "GO" "GOALL"
# [13] "IPI" "MAP" "OMIM" "ONTOLOGY" "ONTOLOGYALL" "PATH"
# [19] "PFAM" "PMID" "PROSITE" "REFSEQ" "SYMBOL" "UCSCKG"
# [25] "UNIGENE" "UNIPROT"
查看数据库或数据表的键
# keys 返回数据库或表的键
head(keys(org.Hs.eg.db))
# [1] "1" "2" "3" "9" "10" "11"
head(keys(org.Hs.eg.db, keytype = 'SYMBOL'))
# [1] "A1BG" "A2M" "A2MP1" "NAT1" "NAT2" "NATP"
好了,看完了这些信息,我们就可以开工啦!
先读取想要转换的基因的 symbol
# read gene symbol
symbol <- read.table(file = '~/Downloads/symbol.txt', sep = '\t', header = FALSE)
symbol <- as.character(unique(symbol$V1))
读取完成,将 symbol
转换为 entrezid
# 将 symbol 对应到 entrezid
entrezid <- select(org.Hs.eg.db, keys=symbol, columns = 'ENTREZID', keytype = 'SYMBOL')
# 'select()' returned 1:1 mapping between keys and columns
可以看到最后的输出信息,表示是一对一匹配的
那到这是不是就结束了呢,我们来看看结果
SYMBOL ENTREZID
1 COL10A1 1300
2 CTHRC1 115908
3 POSTN 10631
4 COL11A1 1301
... ...
120 MURC <NA>
121 H2AFX <NA>
122 HIST1H1T <NA>
123 C14orf80 <NA>
咦,怎么没匹配到 ID
呢,这可咋办呢。
在这里,我们就要引出一个基因 “别名(alias
)”:
通常,基因 symbol
是由 HUGO(Human Genome Organisation)
基因命名法给出的权威性的命名,但是在这之前,许多研究中对基因的命名并没有那么规范,不同研究中可能会对同一个基因有不同的称呼,其中一些名称已经被广泛使用,
因此会存在一个基因或其对应的蛋白质会有不同的别名,不同的别名可能会对应于同一个基因,这种一对多或多对一的关系。
详情请自行维基百科:Gene nomenclature
好了,既然 symbol
找不对,那就试试 alias
吧
# 是否存在未匹配的 SYMBOL
no_map <- sort(as.character(entrezid[is.na(entrezid$ENTREZID),'SYMBOL']))
先把未匹配上的基因挑出来
# 进一步查看是否是基因别名 alias
alias <- select(org.Hs.eg.db, keys=no_map, columns = c('SYMBOL', 'ENTREZID'), keytype = 'ALIAS')
# 'select()' returned 1:many mapping between keys and columns
我们把 keytype
换成了 ALIAS
,与 keys
参数,也就是我们认为是别名的基因。
然后要对应到的是 SYMBOL
和 ENTREZID
。
看看输出信息,many mapping
?出现多对一了?
看看 alias
长啥样
# >alias
#
# ALIAS SYMBOL ENTREZID
# 1 FAM63A MINDY1 55793
# 2 FAM129B NIBAN2 64855
# 3 MB21D1 CGAS 115004
# 4 AIM1 CRYBG1 202
# 5 AIM1 AURKB 9212
# 6 AIM1 SLC45A2 51151
# 7 TMEM57 MACO1 55219
# 8 WISP1 CCN4 8840
# 9 PYCRL PYCR3 65263
# 10 C16orf59 TEDC2 80178
# 11 SDCCAG3 ENTR1 10807
# 12 GATSL3 CASTOR1 652968
# 13 C11orf84 SPINDOC 144097
# 14 DOPEY2 DOP1B 9980
# 15 AIM1L CRYBG2 55057
# 16 FAM109A PHETA1 144717
# 17 TMEM2 CEMIP2 23670
# 18 KIAA1524 CIP2A 57650
# 19 FAM64A PIMREG 54478
# 20 GSG2 HASPIN 83903
# 21 KIAA1468 RELCH 57614
# 22 MURC CAVIN4 347273
# 23 H2AFX H2AX 3014
# 24 HIST1H1T H1-6 3010
# 25 C14orf80 TEDC1 283643
可以看到 4-6
行输出结果,别名 AIM1
对应到了 3
个基因 symbol
确实出现了我们上面说到的情况。那这种情况要怎么处理呢?
一般对我来说,我会选择删掉,毕竟这种无法确定这个基因别名到底对应的是哪个 symbol
# 删除多重配对的结果
uni_alias <- mapIds(org.Hs.eg.db, keys = no_map, column = 'SYMBOL', keytype = 'ALIAS', multiVals = 'filter')
我们使用 mapIds
,用法和 select
差不多,并设置 multiVals='filter'
,意思是删除这些重复匹配,你也可以设置其他值,如 first
保留第一个值等等。
最后返回的 uni_alias
为删除多匹配结果的 symbol
# 重新匹配到 id
alias_symbol_id <- select(org.Hs.eg.db, keys = uni_alias, columns = 'ENTREZID', keytype = 'SYMBOL')
# 'select()' returned 1:1 mapping between keys and columns
从输出信息可以看出,已经变成一对一了
最后,将两个结果合并,并输出
# 合并结果
res <- rbind(entrezid[!is.na(entrezid$ENTREZID),], alias_symbol_id)
# 输出结果
write.table(res, file = '~/Downloads/symbol_id.txt', sep = '\t', row.names = FALSE)
2.2 Python
Python
版本的话,作为一个进阶。下面我就简单介绍一下我之前用过的方法。
我之前是直接去 NCBI ftp
,下载对应的基因信息文件,然后利用正则表达式提取自己想要的信息,重新存为一个 Excel
。如 id
和 symbol
或其他像 ensemble
等基因或蛋白质的信息。
需要的时候,直接从存储的文件中进行匹配。这些操作比较复杂,感兴趣的可以私聊。
下面我就直接把前面安装 R
包的时候下载的文件拿来用了,加入一些数据库查询语句,简单匹配一下,大家作为例子了解一下
import pandas as pd
import sqlite3
# org.Hs.eg.db 包中的 sqlite 数据文件
db = "org.Hs.eg.db/extdata/org.Hs.eg.sqlite"
# 建立连接
conn = sqlite3.connect(db)
导入模块,并对数据文件建立连接
查询文件中所包含的所有表
pd.read_sql('select * from sqlite_master where type="table"', con=conn)
查询文件中所包含的所有视图
pd.read_sql('select * from sqlite_master where type="view"', con=conn)
查询文件中所包含的所有索引
pd.read_sql('select * from sqlite_master where type="index"', con=conn)
可以看到,类似 R
,存在许多表,例如
pd.read_sql('select * from gene_info', con=conn)
获取基因 symbol
及其 id
df = pd.read_sql('select gene_id,symbol from gene_info inner join genes on gene_info._id = genes._id', con=conn)
type(df)
# pandas.core.frame.DataFrame
最后,这就变成一个 pandas DataFrame
格式数据了
symbol = pd.read_csv('~/Downloads/symbol.txt', header=None, names=['symbol'])
df.loc[df.symbol.isin(symbol.symbol)]
可以看到匹配到了 100
个基因
后续代码
# 获取基因 symbol、别名列表
alias = pd.read_sql('select symbol, alias_symbol from alias inner join gene_info on alias._id = gene_info._id', con=conn)
# 获取为匹配的别名
no_map = symbol.loc[~symbol.symbol.isin(entrezid.symbol)]
# 未匹配的别名再匹配到 symbol
tmp = alias.loc[alias.alias_symbol.isin(no_map.symbol)]
left_symbol = tmp.loc[tmp.alias_symbol.isin(tmp.alias_symbol.drop_duplicates(keep=False))]
# 再用 symbol 匹配 id
left_id = df.loc[df.symbol.isin(left_symbol.symbol)]
# 合并输出并输出
res = pd.concat([entrezid, left_id])
res.to_csv('~/Downloads/symbol_id.p.txt', index=
大功告成!