基因 ID 和 Symbol 转换

前言

看了这么多期的 circos,都有点乏了,换换口味也是好的。

下面,咱就开始吧。

做生信分析,总是免不了要给基因 IDSymbol 转换来转换去。

方法

一般要进行 IDSymbol 的转换呢,主要有两种方式:

  1. 网站提供的工具,比如 biodbnet
  2. 编写代码

1. 用网站转换

如果不会编写代码的话,可以使用这个网站 biodbnet

image.png

这种方式比较简单,比如上面的例子,我们输入的是人类(9606)基因 symbol,需要对应的基因 id,提交之后

image.png

可以下载转换的结果。

但是以我的经验来说,这个网站如果输入的基因很多,速度非常慢,而且很多基因 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 参数,也就是我们认为是别名的基因。

然后要对应到的是 SYMBOLENTREZID

看看输出信息,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。如 idsymbol 或其他像 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)
image.png

查询文件中所包含的所有视图

pd.read_sql('select * from sqlite_master where type="view"', con=conn)
image.png

查询文件中所包含的所有索引

pd.read_sql('select * from sqlite_master where type="index"', con=conn)
image.png

可以看到,类似 R,存在许多表,例如

pd.read_sql('select * from gene_info', con=conn)
image.png

获取基因 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 格式数据了

image.png

symbol = pd.read_csv('~/Downloads/symbol.txt', header=None, names=['symbol'])
df.loc[df.symbol.isin(symbol.symbol)]
image.png

可以看到匹配到了 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=
image.png

大功告成!

©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 213,014评论 6 492
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 90,796评论 3 386
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 158,484评论 0 348
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 56,830评论 1 285
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 65,946评论 6 386
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 50,114评论 1 292
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 39,182评论 3 412
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 37,927评论 0 268
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 44,369评论 1 303
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 36,678评论 2 327
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 38,832评论 1 341
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 34,533评论 4 335
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 40,166评论 3 317
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 30,885评论 0 21
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 32,128评论 1 267
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 46,659评论 2 362
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 43,738评论 2 351

推荐阅读更多精彩内容