前言
看了這么多期的 circos，都有點(diǎn)乏了，換換口味也是好的。
下面，咱就開始吧。
做生信分析，總是免不了要給基因 ID 和 Symbol 轉(zhuǎn)換來轉(zhuǎn)換去。
方法
一般要進(jìn)行 ID 和 Symbol 的轉(zhuǎn)換呢，主要有兩種方式：
- 網(wǎng)站提供的工具，比如
biodbnet - 編寫代碼
1. 用網(wǎng)站轉(zhuǎn)換
如果不會編寫代碼的話，可以使用這個(gè)網(wǎng)站 biodbnet
這種方式比較簡單，比如上面的例子，我們輸入的是人類（9606）基因 symbol，需要對應(yīng)的基因 id，提交之后
可以下載轉(zhuǎn)換的結(jié)果。
但是以我的經(jīng)驗(yàn)來說，這個(gè)網(wǎng)站如果輸入的基因很多，速度非常慢，而且很多基因 symbol 無法轉(zhuǎn)換到 id 的，所以對于有編程基礎(chǔ)的朋友，并不推薦這種方式
2. 編程實(shí)現(xiàn)
編程的話，很多語言都可以實(shí)現(xiàn)，看自己比較喜歡，比較擅長用什么語言
下面主要介紹一下 R 以及 Python 兩種語言實(shí)現(xiàn)方式
2.1 R
用 R 實(shí)現(xiàn)的話，一般都是使用 org.Hs.eg.db 這個(gè)模塊提供的數(shù)據(jù)來進(jìn)行轉(zhuǎn)換
安裝和導(dǎo)入
# install
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("org.Hs.eg.db")
# library
library('org.Hs.eg.db')
安裝這個(gè)包之后，對于的路徑下面會有一個(gè) org.Hs.eg.sqlite 文件，存儲了人類基因數(shù)據(jù)，后面的各種轉(zhuǎn)換其實(shí)都是對這個(gè)文件進(jìn)行操作。
查看基本信息
# 獲取所有可用的表
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"
從上面的輸出信息可以看出，包含了很多數(shù)據(jù)表，如 ENSEMBL、ENTREZID、SYMBOL 等
# keytype 配合 keys 使用，在 select 函數(shù)中匹配 keys 參數(shù)指定的 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"
查看數(shù)據(jù)庫或數(shù)據(jù)表的鍵
# keys 返回?cái)?shù)據(jù)庫或表的鍵
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"
好了，看完了這些信息，我們就可以開工啦！
先讀取想要轉(zhuǎn)換的基因的 symbol
# read gene symbol
symbol <- read.table(file = '~/Downloads/symbol.txt', sep = '\t', header = FALSE)
symbol <- as.character(unique(symbol$V1))
讀取完成，將 symbol 轉(zhuǎn)換為 entrezid
# 將 symbol 對應(yīng)到 entrezid
entrezid <- select(org.Hs.eg.db, keys=symbol, columns = 'ENTREZID', keytype = 'SYMBOL')
# 'select()' returned 1:1 mapping between keys and columns
可以看到最后的輸出信息，表示是一對一匹配的
那到這是不是就結(jié)束了呢，我們來看看結(jié)果
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 呢，這可咋辦呢。
在這里，我們就要引出一個(gè)基因 “別名（alias）”：
通常，基因 symbol 是由 HUGO(Human Genome Organisation) 基因命名法給出的權(quán)威性的命名，但是在這之前，許多研究中對基因的命名并沒有那么規(guī)范，不同研究中可能會對同一個(gè)基因有不同的稱呼，其中一些名稱已經(jīng)被廣泛使用，
因此會存在一個(gè)基因或其對應(yīng)的蛋白質(zhì)會有不同的別名，不同的別名可能會對應(yīng)于同一個(gè)基因，這種一對多或多對一的關(guān)系。
詳情請自行維基百科：Gene nomenclature
好了，既然 symbol 找不對，那就試試 alias 吧
# 是否存在未匹配的 SYMBOL
no_map <- sort(as.character(entrezid[is.na(entrezid$ENTREZID),'SYMBOL']))
先把未匹配上的基因挑出來
# 進(jìn)一步查看是否是基因別名 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 參數(shù)，也就是我們認(rèn)為是別名的基因。
然后要對應(yīng)到的是 SYMBOL 和 ENTREZID。
看看輸出信息，many mapping？出現(xiàn)多對一了？
看看 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 行輸出結(jié)果，別名 AIM1 對應(yīng)到了 3 個(gè)基因 symbol
確實(shí)出現(xiàn)了我們上面說到的情況。那這種情況要怎么處理呢？
一般對我來說，我會選擇刪掉，畢竟這種無法確定這個(gè)基因別名到底對應(yīng)的是哪個(gè) symbol
# 刪除多重配對的結(jié)果
uni_alias <- mapIds(org.Hs.eg.db, keys = no_map, column = 'SYMBOL', keytype = 'ALIAS', multiVals = 'filter')
我們使用 mapIds，用法和 select 差不多，并設(shè)置 multiVals='filter'，意思是刪除這些重復(fù)匹配，你也可以設(shè)置其他值，如 first 保留第一個(gè)值等等。
最后返回的 uni_alias 為刪除多匹配結(jié)果的 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
從輸出信息可以看出，已經(jīng)變成一對一了
最后，將兩個(gè)結(jié)果合并，并輸出
# 合并結(jié)果
res <- rbind(entrezid[!is.na(entrezid$ENTREZID),], alias_symbol_id)
# 輸出結(jié)果
write.table(res, file = '~/Downloads/symbol_id.txt', sep = '\t', row.names = FALSE)
2.2 Python
Python 版本的話，作為一個(gè)進(jìn)階。下面我就簡單介紹一下我之前用過的方法。
我之前是直接去 NCBI ftp ，下載對應(yīng)的基因信息文件，然后利用正則表達(dá)式提取自己想要的信息，重新存為一個(gè) Excel。如 id 和 symbol 或其他像 ensemble 等基因或蛋白質(zhì)的信息。
需要的時(shí)候，直接從存儲的文件中進(jìn)行匹配。這些操作比較復(fù)雜，感興趣的可以私聊。
下面我就直接把前面安裝 R 包的時(shí)候下載的文件拿來用了，加入一些數(shù)據(jù)庫查詢語句，簡單匹配一下，大家作為例子了解一下
import pandas as pd
import sqlite3
# org.Hs.eg.db 包中的 sqlite 數(shù)據(jù)文件
db = "org.Hs.eg.db/extdata/org.Hs.eg.sqlite"
# 建立連接
conn = sqlite3.connect(db)
導(dǎo)入模塊，并對數(shù)據(jù)文件建立連接
查詢文件中所包含的所有表
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
最后，這就變成一個(gè) pandas DataFrame 格式數(shù)據(jù)了
symbol = pd.read_csv('~/Downloads/symbol.txt', header=None, names=['symbol'])
df.loc[df.symbol.isin(symbol.symbol)]
可以看到匹配到了 100 個(gè)基因
后續(xù)代碼
# 獲取基因 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=
大功告成！
