ggetの使い方。
論文は以下。
ggetは Ensembl, UniProt, NCBI などの遺伝子配列データベースを対象にした遺伝子名や遺伝子IDによる配列検索、BLAST検索などを簡単に実行するインターフェースを提供してくれるPythonパッケージ。コマンドラインバージョンもあるがここでは解析プログラムにも組み込みやすいPython APIを紹介する。
R言語の場合は使いやすいツールとしてbiomaRtがあるけど、biomaRtが検索時にDBをひとつ固定するのに対して、ggetの場合はデータベース横断的に検索してくれるのが特徴。
ここでは遺伝子検索に関連する一連の操作として、以下の基本的な検索、
- リファレンス生物種の情報取得（ gget ref ）
- フリーテキストによる遺伝子検索（ gget search ）
- 遺伝子メタ情報の取得（ gget info ）
- 遺伝子配列の取得（ gget seq ）
- BLAST検索（ gget blast ）
さらに、発展的な検索？というより外部ツールの簡単な利用インターフェースとして提供されている以下、
- マッピングによるゲノム座標取得（ gget blat ）
- KEGGやGOによるエンリッチメント解析（ gget enrichr ）
- ARCHS4を利用した発現相関遺伝子の探索（ gget archs4 ）
- タンパク質立体構造の取得（ gget pdb ）
などの機能も紹介する。
インストール
pip install gget
でオーケー。パッケージ依存性は最小限で、依存してるのもbeautifulsoup4, MySQL-connector, requestsくらいなのでインストールも簡単。
基本的な検索
import gget
gget ref : リファレンス生物種の情報取得
リファレンス生物種のメタ情報の取得、ゲノム配列や遺伝子アノテーションファイルのダウンロードリンクなどの表示。
利用できるリファレンスゲノムのリストアップ。とくに指定しない場合はEnsemblの最新版から取得。
gget.ref(species=None, list_species=True)[:10]
['acanthochromis_polyacanthus',
'accipiter_nisus',
'ailuropoda_melanoleuca',
'amazona_collaria',
'amphilophus_citrinellus',
'amphiprion_ocellaris',
'amphiprion_percula',
'anabas_testudineus',
'anas_platyrhynchos',
'anas_platyrhynchos_platyrhynchos']
生物種名を指定すると、利用できるファイルのFTP LINKとファイルサイズなどを表示してくれる。
gget.ref(species='mus_musculus')
{'mus_musculus': {'transcriptome_cdna': {'ftp': 'http://ftp.ensembl.org/pub/release-108/fasta/mus_musculus/cdna/Mus_musculus.GRCm39.cdna.all.fa.gz',
'ensembl_release': 108,
'release_date': '2022-10-04',
'release_time': '19:32',
'bytes': '49M'},
'genome_dna': {'ftp': 'http://ftp.ensembl.org/pub/release-108/fasta/mus_musculus/dna/Mus_musculus.GRCm39.dna.primary_assembly.fa.gz',
'ensembl_release': 108,
'release_date': '2022-10-04',
'release_time': '18:37',
'bytes': '769M'},
'annotation_gtf': {'ftp': 'http://ftp.ensembl.org/pub/release-108/gtf/mus_musculus/Mus_musculus.GRCm39.108.gtf.gz',
'ensembl_release': 108,
'release_date': '2022-10-04',
'release_time': '19:16',
'bytes': '31M'},
'coding_seq_cds': {'ftp': 'http://ftp.ensembl.org/pub/release-108/fasta/mus_musculus/cds/Mus_musculus.GRCm39.cds.all.fa.gz',
'ensembl_release': 108,
'release_date': '2022-10-04',
'release_time': '19:32',
'bytes': '16M'},
'non-coding_seq_ncRNA': {'ftp': 'http://ftp.ensembl.org/pub/release-108/fasta/mus_musculus/ncrna/Mus_musculus.GRCm39.ncrna.fa.gz',
'ensembl_release': 108,
'release_date': '2022-10-04',
'release_time': '19:45',
'bytes': '7.6M'},
'protein_translation_pep': {'ftp': 'http://ftp.ensembl.org/pub/release-108/fasta/mus_musculus/pep/Mus_musculus.GRCm39.pep.all.fa.gz',
'ensembl_release': 108,
'release_date': '2022-10-04',
'release_time': '19:32',
'bytes': '11M'}}}
ダウンロードする場合は以下のように指定すると簡単。whichで必要なファイルだけ指定。ftpフラグ立てるとリンクURLだけをリストで返す。
dl_links = gget.ref(species='mus_musculus',
which=['ncrna'],
ftp=True)
import urllib.request
urllib.request.urlretrieve(dl_links[0],
'./GRCm39_rna.fa.gz')
('./GRCm39_rna.fa.gz', <http.client.HTTPMessage at 0x13d3b7430>)
gget search : フリーテキストによる遺伝子検索
適当なフリーテキストでEnsemblを検索する。結果はpandasのデータフレームで返す。
gget.search(searchwords='ACE2', species='homo_sapiens')
| ensembl_id | gene_name | ensembl_description | ext_ref_description | biotype | url | |
|---|---|---|---|---|---|---|
| 0 | ENSG00000130234 | ACE2 | angiotensin converting enzyme 2 [Source:HGNC S... | angiotensin converting enzyme 2 | protein_coding | https://uswest.ensembl.org/homo_sapiens/Gene/S... |
| 1 | ENSG00000182240 | BACE2 | beta-secretase 2 [Source:HGNC Symbol;Acc:HGNC:... | beta-secretase 2 | protein_coding | https://uswest.ensembl.org/homo_sapiens/Gene/S... |
| 2 | ENSG00000224388 | BACE2-IT1 | BACE2 intronic transcript 1 [Source:HGNC Symbo... | BACE2 intronic transcript 1 | lncRNA | https://uswest.ensembl.org/homo_sapiens/Gene/S... |
もう少し長い英文でも可能。その場合は全単語を含むか、いずれかを含むか、AND/ORをコントロールできる。
df = gget.search(searchwords='angiotensin converting enzyme 2',
species='homo_sapiens',
andor='and',
limit=1)
display(df)
| ensembl_id | gene_name | ensembl_description | ext_ref_description | biotype | url | |
|---|---|---|---|---|---|---|
| 0 | ENSG00000130234 | ACE2 | angiotensin converting enzyme 2 [Source:HGNC S... | angiotensin converting enzyme 2 | protein_coding | https://uswest.ensembl.org/homo_sapiens/Gene/S... |
ens_id = df['ensembl_id'][0]
ens_id
'ENSG00000130234'
gget info : 遺伝子メタ情報の取得
Ensembl IDを利用して、Ensembl, UniProt, NCBIから横断的に遺伝子のメタ情報を取得する。結果はpandasのデータフレームで返ってくる。
df = gget.info(ens_ids=[ens_id])
以下のような情報が含まれる。この情報使うとensembl, uniprot, pdbなどとのID変換も容易。
df.columns
Index(['ensembl_id', 'uniprot_id', 'pdb_id', 'ncbi_gene_id', 'species',
'assembly_name', 'primary_gene_name', 'ensembl_gene_name', 'synonyms',
'parent_gene', 'protein_names', 'ensembl_description',
'uniprot_description', 'ncbi_description', 'subcellular_localisation',
'object_type', 'biotype', 'canonical_transcript', 'seq_region_name',
'strand', 'start', 'end', 'all_transcripts', 'transcript_biotypes',
'transcript_names', 'transcript_strands', 'transcript_starts',
'transcript_ends', 'all_exons', 'exon_starts', 'exon_ends',
'all_translations', 'translation_starts', 'translation_ends'],
dtype='object')
たとえば以下のような感じ。
df[['ensembl_gene_name', 'uniprot_id', 'uniprot_description']]
| ensembl_gene_name | uniprot_id | uniprot_description | |
|---|---|---|---|
| ENSG00000130234 | ACE2 | Q9BYF1 | (Microbial infection) Acts as a receptor for h... |
gget seq : 遺伝子配列の取得
Ensembl IDを利用して、遺伝子の配列をFASTAで返す。
result = gget.seq(ens_ids=[ens_id])
print(result[0])
print(result[1][:60])
>ENSG00000130234 chromosome:GRCh38:X:15494566:15607236:-1
GTGTTTGTCTTGACTGGCTAGCAACTTAGAACTTTTTAAAAGAGGCAAAGGCAGAGGAGA
塩基配列、タンパク質配列どちらも取得できる。translateフラグを立てるとタンパク質配列。塩基配列の場合はEnsemblから、タンパク質配列の場合はUniProtから取得。
アイソフォームを全部取得する場合はisoformsフラグを立てる。saveフラグ立てるとカレントディレクトリにFASTAファイルを保存する。
result = gget.seq(ens_ids=[ens_id],
translate=True,
isoforms=True,
save=True)
gget blast : BLAST検索
たとえば以下の配列をBLAST検索してみる。
test_seq = result[1]
test_seq[:60]
'MSSSSWLLLSLVAVTAAQSTIEEQAKTFLDKFNHEAEDLFYQSSLASWNYNTNITEENVQ'
sequenceに遺伝子の配列、あるいは検索するFASTAファイルのパスを指定する。
program は、'blastn', 'blastp', 'blastx', 'tblastn', 'tblastx' から選び、
database は、'nt', 'nr', 'refseq_rna', 'refseq_protein', 'swissprot', 'pdbaa', 'pdbnt' から選ぶ。
デフォルトでは、塩基配列の場合、 'nt' に対する 'blastn' 、タンパク質配列の場合、 'nr' に対する 'blastp' となる。
トップNヒットに限定、E-valueによるフィルタリングも可能。詳細は help(gget.blast) で。
結果はBLAST結果を含むpandasのデータフレームで返ってくる。
result = gget.blast(test_seq)
result[:5]
| Description | Scientific Name | Common Name | Taxid | Max Score | Total Score | Query Cover | E value | Per. Ident | Acc. Len | Accession | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | angiotensin-converting enzyme 2 isoform 1 prec... | Homo sapiens | human | 9606 | 1614 | 1614 | 99% | 0.0 | 98.60% | 805 | NP_001358344.1 |
| 1 | angiotensin I converting enzyme (peptidyl-dipe... | synthetic construct | NaN | 32630 | 1612 | 1612 | 99% | 0.0 | 98.47% | 805 | BAJ21180.1 |
| 2 | ACE2 [Homo sapiens] | Homo sapiens | human | 9606 | 1612 | 1612 | 99% | 0.0 | 98.47% | 805 | BAB40370.1 |
| 3 | human ACE2/FLAGA(R) fusion protein [Vector pLe... | Vector pLenti-hACE2 | NaN | 2913627 | 1612 | 1612 | 99% | 0.0 | 98.60% | 813 | UJZ92616.1 |
| 4 | unnamed protein product [Homo sapiens] | Homo sapiens | human | 9606 | 1612 | 1612 | 99% | 0.0 | 98.47% | 805 | BAG37592.1 |
発展的な検索
gget blat : マッピングによるゲノム座標取得
BLATによるマッピングを利用して、なんらかの配列のゲノム上の座標を推定する。
内部的には、UCSCのウェブBLATを利用してるだけなので、リファレンスゲノムに制限がある。
たとえば以下のACE2遺伝子について。
seq = gget.seq(ens_ids=[ens_id])
seq = seq[1][:60]
print(seq)
GTGTTTGTCTTGACTGGCTAGCAACTTAGAACTTTTTAAAAGAGGCAAAGGCAGAGGAGA
result = gget.blat(sequence=seq,
assembly='human')
結果はpandasのデータフレームで返る。
result[:3]
| genome | query_size | aligned_start | aligned_end | matches | mismatches | %_aligned | %_matched | chromosome | strand | start | end | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | hg38 | 60 | 1 | 60 | 60 | 0 | 100.0 | 100.00 | chrX | - | 15607177 | 15607236 |
| 1 | hg38 | 60 | 1 | 60 | 59 | 1 | 100.0 | 98.33 | chr4 | - | 79206646 | 79206705 |
| 2 | hg38 | 60 | 1 | 60 | 58 | 2 | 100.0 | 96.67 | chr4_GL000257v2_alt | - | 531132 | 531191 |
gget enrichr : KEGGやGOによるエンリッチメント解析
遺伝子シンボルのリストで、EnrichrウェブAPIへのジョブサブミットを実行する。機能カテゴリのエンリッチメントをざっくりと確認したいときに便利。
たとえば以下のような適当な遺伝子セットを用意する。
df = gget.search(searchwords='MAPK',
species='homo_sapiens')
df[:3]
| ensembl_id | gene_name | ensembl_description | ext_ref_description | biotype | url | |
|---|---|---|---|---|---|---|
| 0 | ENSG00000008735 | MAPK8IP2 | mitogen-activated protein kinase 8 interacting... | mitogen-activated protein kinase 8 interacting... | protein_coding | https://uswest.ensembl.org/homo_sapiens/Gene/S... |
| 1 | ENSG00000050130 | JKAMP | JNK1/MAPK8 associated membrane protein [Source... | JNK1/MAPK8 associated membrane protein | protein_coding | https://uswest.ensembl.org/homo_sapiens/Gene/S... |
| 2 | ENSG00000050748 | MAPK9 | mitogen-activated protein kinase 9 [Source:HGN... | mitogen-activated protein kinase 9 | protein_coding | https://uswest.ensembl.org/homo_sapiens/Gene/S... |
遺伝子のシンボルを渡して、データベースを設定して、Enrichrを実行する。
バックグラウンドのデータベースはEnrichrが提供している一覧から選べる。
plot=True とすると、簡単なプロットを出力してくれる。
gget.enrichr(genes=df['gene_name'].values,
database='pathway',
plot=True)
| rank | path_name | p_val | z_score | combined_score | overlapping_genes | adj_p_val | database | |
|---|---|---|---|---|---|---|---|---|
| 0 | 1 | MAPK signaling pathway | 6.389385e-22 | 39.975187 | 1950.878438 | [MAPK8IP2, MAPK14, MAPK8IP3, MAPK12, MAPK8IP1,... | 8.625670e-20 | KEGG_2021_Human |
| 1 | 2 | IL-17 signaling pathway | 1.661678e-19 | 79.709259 | 3446.731116 | [MAPK14, MAPK15, MAPK12, MAPK13, MAPK10, MAPK1... | 1.121632e-17 | KEGG_2021_Human |
| 2 | 3 | Neurotrophin signaling pathway | 1.274299e-14 | 48.110450 | 1539.235873 | [MAPK10, MAPK11, MAPK9, MAPK8, MAPK7, MAPKAPK2... | 5.734345e-13 | KEGG_2021_Human |
| 3 | 4 | GnRH signaling pathway | 4.753301e-14 | 55.657047 | 1707.410809 | [MAPK10, MAPK11, MAPK9, MAPK8, MAPK7, MAPK1, M... | 1.604239e-12 | KEGG_2021_Human |
| 4 | 5 | Fc epsilon RI signaling pathway | 1.378735e-13 | 68.949153 | 2041.752629 | [MAPK10, MAPK11, MAPK9, MAPK8, MAPK1, MAPK14, ... | 3.381377e-12 | KEGG_2021_Human |
| ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 130 | 131 | Chemokine signaling pathway | 9.194986e-02 | 4.077812 | 9.731747 | [MAPK1, MAPK3] | 9.475749e-02 | KEGG_2021_Human |
| 131 | 132 | Regulation of actin cytoskeleton | 1.136693e-01 | 3.582244 | 7.789452 | [MAPK1, MAPK3] | 1.162527e-01 | KEGG_2021_Human |
| 132 | 133 | Chemical carcinogenesis | 1.320810e-01 | 3.261355 | 6.602091 | [MAPK1, MAPK3] | 1.340672e-01 | KEGG_2021_Human |
| 133 | 134 | Human papillomavirus infection | 2.187139e-01 | 2.338399 | 3.554345 | [MAPK1, MAPK3] | 2.203461e-01 | KEGG_2021_Human |
| 134 | 135 | PI3K-Akt signaling pathway | 2.412301e-01 | 2.183044 | 3.104297 | [MAPK1, MAPK3] | 2.412301e-01 | KEGG_2021_Human |
135 rows × 8 columns
gget archs4 : ARCHS4を利用した発現相関遺伝子の探索
ARCHS4検索の実行。
統一パイプラインで処理された大量の公共RNA-seqデータ（ヒトかマウス）を対象に、全サンプルで発現量相関の高い遺伝子を検索する。
gget.archs4('ACE2')[:10]
| gene_symbol | pearson_correlation | |
|---|---|---|
| 1 | SLC5A1 | 0.579634 |
| 2 | CYP2C18 | 0.576577 |
| 3 | ADH1C | 0.565597 |
| 4 | SLC13A2 | 0.563400 |
| 5 | TSPAN8 | 0.562570 |
| 6 | GSTA1 | 0.561725 |
| 7 | MOGAT2 | 0.557799 |
| 8 | SLC44A3 | 0.554614 |
| 9 | S100A14 | 0.553783 |
| 10 | HMGCS2 | 0.553401 |
gget pdb : タンパク質立体構造の取得
PDBからタンパク質の構造情報を取得。PDB IDが必要なので、遺伝子名から検索する場合はたとえば以下のようにPDB IDまで辿ってダウンロードする。
ens_id = gget.search(searchwords='ACE2', species='homo_sapiens').loc[0, 'ensembl_id']
pdb_id = gget.info(ens_ids=ens_id)['pdb_id'].values[0][0]
result = gget.pdb(pdb_id=pdb_id,
save=True)
!head ./7fdi.pdb
HEADER VIRAL PROTEIN 16-JUL-21 7FDI
TITLE SARS-COV-2 SPIKE RBDMACSP36 BINDING TO HACE2
COMPND MOL_ID: 1;
COMPND 2 MOLECULE: ANGIOTENSIN-CONVERTING ENZYME 2;
COMPND 3 CHAIN: A;
COMPND 4 SYNONYM: ANGIOTENSIN-CONVERTING ENZYME HOMOLOG,ACEH,ANGIOTENSIN-
COMPND 5 CONVERTING ENZYME-RELATED CARBOXYPEPTIDASE,ACE-RELATED
COMPND 6 CARBOXYPEPTIDASE,METALLOPROTEASE MPROT15;
COMPND 7 EC: 3.4.17.23,3.4.17.-;
COMPND 8 ENGINEERED: YES;
その他
わざわざggetから使う必要もない気もするけど、MUSCLEによるマルチプルアラインメントや、Colab版と同等のalphafoldも実行できる。
それぞれ詳細はドキュメントを参照。
結局
gget.search, gget.infoは日常的に使えるかも。
あとはさくっとnrにBLASTかけたいときはgget.blastするとか。でもEnsemblベースだからせっかくBLASTしても当たった先の配列とってきたりできない。バクテリア関連の解析も対象外。
今後の拡張に期待。