ggetの使い方。

論文は以下。

Laura Luebbert, Lior Pachter, Efficient querying of genomic reference databases with gget, Bioinformatics, 2023;, btac836, https://doi.org/10.1093/bioinformatics/btac836

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	Per. Ident	Acc. Len	Accession
0	angiotensin-converting enzyme 2 isoform 1 prec...	Homo sapiens	human	9606	1614	1614	99%	98.60%	805	NP_001358344.1
1	angiotensin I converting enzyme (peptidyl-dipe...	synthetic construct	NaN	32630	1612	1612	99%	98.47%	805	BAJ21180.1
2	ACE2 [Homo sapiens]	Homo sapiens	human	9606	1612	1612	99%	98.47%	805	BAB40370.1
3	human ACE2/FLAGA(R) fusion protein [Vector pLe...	Vector pLenti-hACE2	NaN	2913627	1612	1612	99%	98.60%	813	UJZ92616.1
4	unnamed protein product [Homo sapiens]	Homo sapiens	human	9606	1612	1612	99%	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しても当たった先の配列とってきたりできない。バクテリア関連の解析も対象外。

今後の拡張に期待。