基因 ID 比對利器

一、背景

對于每個生物資訊分析的人來說，ID 比對（映射）是一項非常常見，但又很繁瑣的任務。假設，我們有一個來自上遊分析的 gene symbol 或報告的 ID 清單，然後我們的下一個分析卻需要使用基因 ID（例如 Entrez gene id 或 Ensembl gene id）。這時候，我們就希望将基因符号或報告的 ID 的清單轉換為相應的基因 ID。

在開始介紹今天的主角 mygene 前，我們先來認識一下 MyGene.info。

MyGene.info 是一個由 NIH(美國國立衛生研究院)/NIGMS 資助，用于提供簡單易用的 REST Web 服務來查詢/檢索基因注釋資料的 API。 MyGene.info 目前包含了 NCBI Entrez、Ensembl、Uniprot、UCSC 在内的 20 多個資料庫，MyGene.info 會每周從這些資料庫中進行資料更新。雖然 MyGene.info 中包含的各個資料源可能有資料使用限制，但 MyGene.info 本身的服務是免費的，其源碼托管在：https://github.com/biothings/mygene.info。

MyGene.info 提供兩種簡單的 Web 服務：一種用于基因查詢，另一種用于基因注釋檢索。兩者都以 JSON 格式傳回結果。截至 2019 年 3 月 6 日，MyGene.info 最新的 API 為 v3，相比 v1、v2，v3 新增了以下幾點内容：

• Refseq accession number now contains version
• "ensembl", "refseq" and "accession" contains associations between RNA and protein
• Better mapping between Ensembl and Entrez gene IDs
• JSON structure slightly changed
• and more bugfixes

雖然 MyGene.info 是一個線上的 web 服務，但它同時也提供了基于 R 和 Python 的第三方子產品，源碼均在 GitHub 上開源：

• MyGene R Client：https://github.com/biothings/mygene.R
• MyGene Python Client：https://github.com/biothings/mygene.py

在這裡，我們簡單展示如何在 Python 中使用 mygene 子產品來快速輕松地進行類似的 ID 比對（映射）。mygene  本質上是一個友善的 Python 子產品，通過這個子產品我們可以通路 MyGene.info 的基因查詢 Web 服務。

二、mygene 安裝與使用

1. 安裝 mygene

在 Python 中 mygene 的安裝非常簡單，直接使用 pip 就可以安裝。

mygene 安裝完成後，現在我們隻需要導入它并執行個體化 mygeneinfo 類：

import mygene

mg = mygene.MyGeneInfo()
           

2. 把 gene symbols 轉換成 Entrez gene ids

假設 xli 是我們要轉換為 entrez gene id 的 gene symbol 清單：

然後我們調用 

querymany

 方法，并告訴它我們輸入的是 "符号(symbol)"，我們想要傳回的是 "entrezgene"(Entrez gene ids)：

上面程式的比對（映射）結果作為字典清單傳回。每個字典都包含我們要求傳回的字段，在本例中為 "entrezgene" 字段。 每個字典還傳回比對的查詢詞 "query" 和内部 id("_id")，大部分時間與 "entrezgene" 相同（如果基因僅來自 Ensembl，則為 ensembl 基因 id）。

3. 把 gene symbols 轉換成 Ensembl gene ids

如果我們隻需要傳回 Ensembl gene ids 時，隻需要把 fields 參數值改成 'ensembl.gene' 即可：

>>> out = mg.querymany(xli, scopes='symbol', fields='ensembl.gene', species='human')
querying 1-9...done.
Finished.

>>> for i in out:
...     print i
...
{u'ensembl': {u'gene': u'ENSG00000150676'}, u'query': u'CCDC83', u'_id': u'220047', u'_score': 87.86632}
{u'ensembl': {u'gene': u'ENSG00000099308'}, u'query': u'MAST3', u'_id': u'23031', u'_score': 88.42681}
{u'ensembl': [{u'gene': u'ENSG00000230143'}, {u'gene': u'ENSG00000236271'}, {u'gene': u'ENSG00000137312'}, {u'gene': u'ENSG00000206379'}, {u'gene': u'ENSG00000232280'}, {u'gene': u'ENSG00000206480'}, {u'gene': u'ENSG00000224740'}, {u'gene': u'ENSG00000223654'}], u'query': u'FLOT1', u'_id': u'10211', u'_score': 90.23538}
{u'ensembl': {u'gene': u'ENSG00000142676'}, u'query': u'RPL11', u'_id': u'6135', u'_score': 82.40764}
{u'ensembl': {u'gene': u'ENSG00000180776'}, u'query': u'ZDHHC20', u'_id': u'253832', u'_score': 87.6894}
{u'ensembl': {u'gene': u'ENSG00000108848'}, u'query': u'LUC7L3', u'_id': u'51747', u'_score': 86.635506}
{u'ensembl': {u'gene': u'ENSG00000277370'}, u'query': u'SNORD49A', u'_id': u'26800', u'_score': 107.55141}
{u'ensembl': {u'gene': u'ENSG00000103811'}, u'query': u'CTSH', u'_id': u'1512', u'_score': 85.88113}
{u'ensembl': {u'gene': u'ENSG00000101473'}, u'query': u'ACOT8', u'_id': u'10005', u'_score': 83.99602}
           

4. ID 與基因不比對

如果輸入 id 沒有比對的基因，mygene 将在螢幕輸出中列印相關的通知。此查詢條目傳回的字典中将包含 "notfound" 值為 True 的關鍵字。

>>> xli = ['CCDC83', 'MAST3', 'FLOT1', 'RPL11', 'Gm10494']
>>> out = mg.querymany(xli, scopes='symbol', fields='entrezgene', species='human')
querying 1-5...done.
Finished.
1 input query terms found no hit:
        [u'Gm10494']
Pass "returnall=True" to return complete lists of duplicate or missing query terms.

>>> for i in out:
...     print(i)
...
{u'entrezgene': u'220047', u'query': u'CCDC83', u'_id': u'220047', u'_score': 87.6894}
{u'entrezgene': u'23031', u'query': u'MAST3', u'_id': u'23031', u'_score': 88.89522}
{u'entrezgene': u'10211', u'query': u'FLOT1', u'_id': u'10211', u'_score': 89.862946}
{u'entrezgene': u'6135', u'query': u'RPL11', u'_id': u'6135', u'_score': 82.584694}
{u'query': u'Gm10494', u'notfound': True}
           

5. 輸入 ID 不僅僅是符号

上面的 xli 清單包含了 symbols, reporters 和 accession numbers，現在我們想要找回 Entrez gene id 和 uniprot id。 幸運的是，mygene 的參數範圍，字段，物種都足夠靈活，可以支援多個值，清單或逗号分隔的字元串：

>>> mg.querymany(xli, scopes='symbol,reporter,accession', fields='entrezgene,uniprot', species='human')
querying 1-8...done.
Finished.
1 input query terms found dup hits:
        [('1007_s_at', 2)]
2 input query terms found no hit:
        ['DDX26B', 'Gm10494']
Pass "returnall=True" to return complete lists of duplicate or missing query terms.
[{'query': 'DDX26B', 'notfound': True}, 
 {'query': 'CCDC83', '_id': '220047', '_score': 88.13916, 'entrezgene': '220047', 'uniprot': {'Swiss-Prot': 'Q8IWF9', 'TrEMBL': 'H0YDV3'}}, 
 {'query': 'MAST3', '_id': '23031', '_score': 88.42681, 'entrezgene': '23031', 'uniprot': {'Swiss-Prot': 'O60307', 'TrEMBL': 'V9GYV0'}}, 
 {'query': 'FLOT1', '_id': '10211', '_score': 90.16039, 'entrezgene': '10211', 'uniprot': {'Swiss-Prot': 'O75955', 'TrEMBL': ['A2AB09', 'Q5ST80', 'A2ABJ5', 'A2AB10', 'A2AB12', 'A2AB13', 'A2AB11']}}, 
 {'query': 'RPL11', '_id': '6135', '_score': 82.44751, 'entrezgene': '6135', 'uniprot': {'Swiss-Prot': 'P62913', 'TrEMBL': ['Q5VVC8', 'Q5VVD0', 'A0A2R8Y447']}}, 
 {'query': 'Gm10494', 'notfound': True}, 
 {'query': '1007_s_at', '_id': '100616237', '_score': 12.442588, 'entrezgene': '100616237'}, 
 {'query': '1007_s_at', '_id': '780', '_score': 11.8290205, 'entrezgene': '780', 'uniprot': {'Swiss-Prot': 'Q08345', 'TrEMBL': ['A0A0A0MSX3', 'A0A024RCL1', 'A0A024RCQ1', 'A0A024RCJ0', 'Q96T61', 'Q96T62', 'A2ABM8', 'A2ABL2', 'A2ABL0', 'E7ERN0', 'A0A0G2JIA2', 'D6RB35', 'A0A0G2JI85', 'D6RAJ3', 'A0A0G2JHK4', 'A0A0G2JJA0', 'H0YAH6', 'A0A140T972', 'E7EUD5', 'E7EXB0', 'E7EPN2', 'E7ETI3', 'E7EVT1', 'E7EVW6', 'A0A0G2JNZ7', 'H0Y717', 'E7ESR9', 'D6R9C4', 'E7EQ23', 'E7EUP7', 'E7EQ30', 'E7EPH4', 'H0Y9F4', 'E7EN94', 'D6RBU7', 'D6RGW5', 'D6RB82', 'E7ETX3', 'E7EX99', 'E7ERI6', 'E7ES06', 'E7ENJ2']}}, 
 {'query': 'AK125780', '_id': '2978', '_score': 10.087812, 'entrezgene': '2978', 'uniprot': {'Swiss-Prot': 'P43080', 'TrEMBL': ['A6PVH5', 'B2R9P6', 'A0A0A0MTF5']}}
]
           

6. ID 比對多個基因

從上一個結果中，你可能已經注意到查詢詞 "1007_s_at" 與兩個基因比對。在這種情況下，我們将從輸出中收到通知，傳回的結果将包括兩個比對的基因。

通過傳遞 "returnall = True"，我們将從傳回的結果中獲得重複或缺少的查詢條目以及比對的輸出：

上面代碼傳回的結果包含了用于比對輸出的清單 "out"，用于缺少查詢項（清單）的 "missing"，以及用于具有多個比對（包括比對數）查詢項的 "dup"。

根據 mygene 的說法，mygene 可以支援大批量的清單 ID 轉換。如傳遞一個大于 1000 ids 的 id 清單（即上面的xli），mygene 将一次批量映射 1000 個 ID，然後将所有結果連接配接起來。是以，從使用者端來看，它與傳遞較短清單完全相同。同時我們不必擔心 mygene 後端伺服器的飽和。

mygene 是一款強大的基因 ID 比對轉換子產品，以 MyGene.info 為背景也可以有更多的玩法，感興趣的可以參考MyGene.info 的官方文檔：http://docs.mygene.info/en/latest/index.html，也歡迎留言交流。

本文分享自微信公衆号 - 生信科技愛好者（bioitee）。

如有侵權，請聯系 [email protected] 删除。

本文參與“OSC源創計劃”，歡迎正在閱讀的你也加入，一起分享。