Takeru Elysia news!¶
大規模GWASやPRS研究を可能にするHail(2)¶
Takeru ElysiaはApache Sparkが利用可能な大規模計算プラットフォームです。従来のジョブスケジューラーを用いたものと異なり、複数台の計算機を使う同じcluster構成でも、次のような特徴を持っています。
- 自分でデータ分割をする必要がない
- 並列処理を自動実行してくれる
- 計算にnfsを使用しないため、ファイルI/Oがボトルネックになりにくい
- 高い拡張性を持っており、数ノードからスタートして数百ノード以上にスケールすることができる
動作するソフトウェアやライブラリの中にはGATKといったメジャーなものもありますが今回は大規模なGWASやPRS(Polygenic Risk Score)を用いた研究ができる、Hailについてお話しします。近年、各種バイオバンクが充実してきており、数年まではできなかったことが可能になってきているようです。Apache Spark上で動作するHailなら、いくつものバイオバンクのデータを読み込み、自在にフィルターや変換、アノテーションをしていくことで新しい知見をより早く得ることができます。
今日の内容: Hail上にplink形式、任意のテキスト形式、bgen形式のファイルを読み込み、データ構造や内容を閲覧してみます。¶
使用するファイルについて:
plink形式と任意のテキストファイルは、PRSのTutorial、
「Basic Tutorial for Polygenic Risk Score Analysis」にあるものを使いました。
https://choishingwan.github.io/PRS-Tutorial/
https://github.com/choishingwan/PRS-Tutorial/raw/master/resources/EUR.zip
をダウンロードし、展開、
次の3つのファイルを読み込みます。
EUR.bed,EUR.bim,EUR.fam
EUR.bedはバイナリで、pedファイルに変換して中を見てみると次のようになっています。
とても横に長いファイルです。
0 HG00096 0 0 1 -9 GA AG AT GA TA ....
EUR.bimは次のような内容が入っています。
$ head EUR.bim
1 rs3094315 0 752566 G A
1 rs2905035 0 775659 A G
1 rs2980319 0 777122 A T
1 rs4040617 0 779322 G A
1 rs2977612 0 780785 T A
1 rs2905062 0 785050 G A
1 rs2980300 0 785989 T C
1 rs3121561 0 990380 T C
1 rs3813193 0 998501 C G
1 rs4075116 0 1003629 C T
EUR.famは次のような内容が入っています。
$ head EUR.fam
0 HG00096 0 0 1 -9
0 HG00097 0 0 2 -9
0 HG00099 0 0 2 -9
0 HG00100 0 0 2 -9
0 HG00101 0 0 1 -9
0 HG00102 0 0 2 -9
0 HG00103 0 0 1 -9
0 HG00105 0 0 1 -9
0 HG00106 0 0 2 -9
0 HG00107 0 0 1 -9
また、任意のテキストファイルとして同じTutorialで使用している、EUR.heightを使います。
$ head EUR.height
FID IID Height
0 HG00096 167.18
0 HG00097 169.85
0 HG00099 170.6
0 HG00100 166.78
0 HG00101 165.03
0 HG00102 170.2
0 HG00103 163.2
0 HG00105 164.92
0 HG00106 174.18
bgenファイルはUK Biobankで配布される形式ですが、次の場所に小さいサンプルデータがありますのでそれを使っています。
https://bitbucket.org/gavinband/bgen/src/release/example/
HailはPython上で動作します。まず、notebook上でHailを動作させます。
この試験ではkubernetes上にspark clusterをデプロイしています。
plinkファイルの読み込み¶
plink_data = hl.import_plink(
bed="data/EUR.bed", bim="data/EUR.bim", fam="data/EUR.fam",
quant_pheno=True, missing='-9')
データの確認¶
hl.import_plinkはplink形式のbed,bim,famファイルをmatrixtableとして読み込みます。
plink_data
データの内容確認¶
まず、データの中身を少し見てみましょう。
hailのmatrixtableは次の4つのfieldsからできています。
- column field
- row field
- entry field
- global fieldです。
まず、データの構造を確認します。
describe¶
plink_data.describe()
show¶
データの中身を見てみます
plink_data.show()
plink_data.entries().show()
plink_data.cols().show()
plink_data.rows().show()
plink_data.globals_table().show()
テキストファイルの読み込み¶
hailではタブ区切りやカンマ区切り、スペース区切りといった任意のテキストファイルを読み込むことができます。
任意の大規模な表現型のデータを簡単に扱うことができます。
delimiter=’\s+’ を指定することでスペース区切りのファイルを読むことができます。
height = hl.import_table('data/EUR.height',
impute=True,
delimiter='\s+',
key='IID')
height.show()
hail table¶
ここで読み込んだデータは、matrixtableではなく、tableと呼ばれる形式で扱われます。
hailには、matrixtableとtableがあります。
height
tableでもfilterなどの操作が可能で、例えば次のようにすると身長が180cm以上のデータに絞り込むことができます。
height.filter(height.Height > 180 ).show()
bgen形式のファイルを読み込み¶
hl.index_bgen("data/example.8bits.bgen",
contig_recoding={"01": "1"},
reference_genome='GRCh37')
bgen = hl.import_bgen("data/example.8bits.bgen",
entry_fields=['GT', 'GP'],
sample_file="data/example.sample")
bgen形式はhailではmatrixtableとして扱われます。
vcf形式やplink形式と同様に、describeやshowを使用して中身を見てみます。
bgen
bgen.describe()
bgen.show()
bgen.col.show()
bgen.row.show()
bgen.entry.show()
今回は次のような操作を行いました。
- plink形式(bed,bim,fam)を読み込み
- 任意のテキストファイルを読み込み
- bgen形式を読み込み
- 各ファイルの閲覧
次回の記事では今回読み込んだplink形式のデータに、任意のテキストファイルとして読み込んだ身長のデータを結合してみようと思います。