大規模GWASやPRS研究を可能にするHail(3)

2020-06-12 00:14未分類

Takeru Elysia news!¶

大規模GWASやPRS研究を可能にするHail(3)¶

Takeru ElysiaはApache Sparkが利用可能な大規模計算プラットフォームです。従来のジョブスケジューラーを用いたものと異なり、複数台の計算機を使う同じcluster構成でも、次のような特徴を持っています。

自分でデータ分割をする必要がない
並列処理を自動実行してくれる
計算にnfsを使用しないため、ファイルI/Oがボトルネックになりにくい
高い拡張性を持っており、数ノードからスタートして数百ノード以上にスケールすることができる

動作するソフトウェアやライブラリの中にはGATKといったメジャーなものもありますが今回は大規模なGWASやPRS(Polygenic Risk Score)を用いた研究ができる、Hailについてお話しします。近年、各種バイオバンクが充実してきており、数年まではできなかったことが可能になってきているようです。Apache Spark上で動作するHailなら、いくつものバイオバンクのデータを読み込み、自在にフィルターや変換、アノテーションをしていくことで新しい知見をより早く得ることができます。

今日の内容: Hail上にplink形式および任意のテキスト形式のファイルを読み込み、データを結合します。¶

使用するファイルについて: 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

HailはPython上で動作します。まず、notebook上でHailを動作させます。この試験ではkubernetes上にspark clusterをデプロイしています。

plinkファイルの読み込み¶

In [2]:

plink_data = hl.import_plink(
    bed="data/EUR.bed", bim="data/EUR.bim", fam="data/EUR.fam",
    quant_pheno=True, missing='-9')

2020-06-12 08:56:09 Hail: INFO: Found 503 samples in fam file.
2020-06-12 08:56:09 Hail: INFO: Found 2395907 variants in bim file.

データの確認¶

hl.import_plinkはplink形式のbed,bim,famファイルをmatrixtableとして読み込みます。

In [3]:

plink_data

Out[3]:

<hail.matrixtable.MatrixTable at 0x7fed76b4b650>

データの内容確認¶

まず、データの中身を少し見てみましょう。 hailのmatrixtableは次の4つのfieldsからできています。

column field
row field
entry field
global fieldです。

まず、データの構造を確認します。

describe¶

In [4]:

plink_data.describe()

----------------------------------------
Global fields:
    None
----------------------------------------
Column fields:
    's': str
    'fam_id': str
    'pat_id': str
    'mat_id': str
    'is_female': bool
    'quant_pheno': float64
----------------------------------------
Row fields:
    'locus': locus<GRCh37>
    'alleles': array<str>
    'rsid': str
    'cm_position': float64
----------------------------------------
Entry fields:
    'GT': call
----------------------------------------
Column key: ['s']
Row key: ['locus', 'alleles']
----------------------------------------

summarize¶

summarizeはデータの概要を示してくれます。

In [5]:

plink_data.summarize()

2020-06-12 08:56:11 Hail: INFO: Coerced sorted dataset

Columns

503 records.

s (str):

Non-missing 503 (100.00%)

Missing 0

Min Size 7

Max Size 7

Mean Size 7.00

Sample Values ['HG00096', 'HG00097', 'HG00099', 'HG00100', 'HG00101']
fam_id (str):

Non-missing 0

Missing 503 (100.00%)
pat_id (str):

Non-missing 0

Missing 503 (100.00%)
mat_id (str):

Non-missing 0

Missing 503 (100.00%)
is_female (bool):

Non-missing 503 (100.00%)

Missing 0

Counts {False: 240, True: 263}
quant_pheno (float64):

Non-missing 0

Missing 503 (100.00%)

2020-06-12 08:56:22 Hail: INFO: Ordering unsorted dataset with network shuffle

Rows

2395907 records.

locus (locus<GRCh37>):

Non-missing	2395907 (100.00%)
Missing	0
Contig Counts	{'X': 22390, '12': 116320, '8': 139501, '19': 33391, '4': 153096, '15': 67018, '11': 121831, '9': 114354, '22': 30599, '13': 96863, '16': 65957, '5': 154779, '10': 129397, '21': 31875, '6': 170124, '1': 179598, '17': 53725, '14': 78549, '20': 58784, '2': 207900, '18': 71918, '7': 133966, '3': 163972}

alleles (array<str>):

Non-missing 2395907 (100.00%)

Missing 0

Min Size 2

Max Size 2

Mean Size 2.00
- alleles[<elements>] (str):
  
  Non-missing 4791814 (100.00%)
  
  Missing 0
  
  Min Size 1
  
  Max Size 41
  
  Mean Size 1.00
  
  Sample Values ['A', 'G', 'G', 'A', 'T']
rsid (str):

Non-missing 2395907 (100.00%)

Missing 0

Min Size 3

Max Size 11

Mean Size 9.21

Sample Values ['rs3094315', 'rs2905035', 'rs2980319', 'rs4040617', 'rs2977612']
cm_position (float64):

Non-missing 2395907 (100.00%)

Missing 0

Minimum 0.00

Maximum 0.00

Mean 0.00

Std Dev 0.00

2020-06-12 08:56:39 Hail: INFO: Ordering unsorted dataset with network shuffle

Entries

1205141221 records.

GT (call):

Non-missing	1202745190 (99.80%)
Missing	2396031 (0.20%)
Homozygous Reference	732575805
Heterozygous	379388567
Homozygous Variant	90780818
Ploidy	{2: 1202745190}
Phased	{False: 1202745190}

show¶

データの中身を見てみます

In [6]:

plink_data.show()

2020-06-12 08:57:59 Hail: INFO: Ordering unsorted dataset with network shuffle
2020-06-12 08:58:07 Hail: INFO: Ordering unsorted dataset with network shuffle

locus	alleles	HG00096.GT	HG00097.GT	HG00099.GT	HG00100.GT
locus<GRCh37>	array<str>	call	call	call	call
1:752566	["A","G"]	0/1	0/0	0/1	0/0
1:775659	["G","A"]	0/1	0/0	0/1	0/0
1:777122	["T","A"]	0/1	0/0	0/1	0/0
1:779322	["A","G"]	0/1	0/0	0/1	0/0
1:780785	["A","T"]	0/1	0/0	0/1	0/0
1:785050	["A","G"]	0/1	0/0	0/1	0/0
1:785989	["C","T"]	0/1	0/0	0/1	0/0
1:990380	["C","T"]	0/1	0/1	0/0	0/1
1:998501	["G","C"]	0/0	0/1	0/0	0/1
1:1003629	["T","C"]	1/1	0/0	0/0	0/0

showing top 10 rows

showing the first 4 of 503 columns

In [7]:

plink_data.entries().show()

2020-06-12 08:58:55 Hail: WARN: entries(): Resulting entries table is sorted by '(row_key, col_key)'.
    To preserve row-major matrix table order, first unkey columns with 'key_cols_by()'
2020-06-12 08:59:03 Hail: INFO: Ordering unsorted dataset with network shuffle

locus	alleles	rsid	cm_position	s	fam_id	pat_id	mat_id	is_female	quant_pheno	GT
locus<GRCh37>	array<str>	str	float64	str	str	str	str	bool	float64	call
1:752566	["A","G"]	"rs3094315"	0.00e+00	"HG00096"	NA	NA	NA	false	NA	0/1
1:752566	["A","G"]	"rs3094315"	0.00e+00	"HG00097"	NA	NA	NA	true	NA	0/0
1:752566	["A","G"]	"rs3094315"	0.00e+00	"HG00099"	NA	NA	NA	true	NA	0/1
1:752566	["A","G"]	"rs3094315"	0.00e+00	"HG00100"	NA	NA	NA	true	NA	0/0
1:752566	["A","G"]	"rs3094315"	0.00e+00	"HG00101"	NA	NA	NA	false	NA	0/0
1:752566	["A","G"]	"rs3094315"	0.00e+00	"HG00102"	NA	NA	NA	true	NA	0/0
1:752566	["A","G"]	"rs3094315"	0.00e+00	"HG00103"	NA	NA	NA	false	NA	NA
1:752566	["A","G"]	"rs3094315"	0.00e+00	"HG00105"	NA	NA	NA	false	NA	0/0
1:752566	["A","G"]	"rs3094315"	0.00e+00	"HG00106"	NA	NA	NA	true	NA	0/0
1:752566	["A","G"]	"rs3094315"	0.00e+00	"HG00107"	NA	NA	NA	false	NA	0/1

showing top 10 rows

In [8]:

plink_data.cols().show()

2020-06-12 08:59:49 Hail: WARN: cols(): Resulting column table is sorted by 'col_key'.
    To preserve matrix table column order, first unkey columns with 'key_cols_by()'

s	fam_id	pat_id	mat_id	is_female	quant_pheno
str	str	str	str	bool	float64
"HG00096"	NA	NA	NA	false	NA
"HG00097"	NA	NA	NA	true	NA
"HG00099"	NA	NA	NA	true	NA
"HG00100"	NA	NA	NA	true	NA
"HG00101"	NA	NA	NA	false	NA
"HG00102"	NA	NA	NA	true	NA
"HG00103"	NA	NA	NA	false	NA
"HG00105"	NA	NA	NA	false	NA
"HG00106"	NA	NA	NA	true	NA
"HG00107"	NA	NA	NA	false	NA

showing top 10 rows

In [9]:

plink_data.rows().show()

2020-06-12 09:00:02 Hail: INFO: Ordering unsorted dataset with network shuffle

locus	alleles	rsid	cm_position
locus<GRCh37>	array<str>	str	float64
1:752566	["A","G"]	"rs3094315"	0.00e+00
1:775659	["G","A"]	"rs2905035"	0.00e+00
1:777122	["T","A"]	"rs2980319"	0.00e+00
1:779322	["A","G"]	"rs4040617"	0.00e+00
1:780785	["A","T"]	"rs2977612"	0.00e+00
1:785050	["A","G"]	"rs2905062"	0.00e+00
1:785989	["C","T"]	"rs2980300"	0.00e+00
1:990380	["C","T"]	"rs3121561"	0.00e+00
1:998501	["G","C"]	"rs3813193"	0.00e+00
1:1003629	["T","C"]	"rs4075116"	0.00e+00

showing top 10 rows

テキストファイルの読み込み¶

hailではタブ区切りやカンマ区切り、スペース区切りといった任意のテキストファイルを読み込むことができます。任意の大規模な表現型のデータを簡単に扱うことができます。

delimiter='\s+' を指定することでスペース区切りのファイルを読むことができます。

In [10]:

height = hl.import_table('data/EUR.height', 
                      impute=True,
                      delimiter='\s+',
                      key='IID')

2020-06-12 09:00:09 Hail: INFO: Reading table to impute column types
2020-06-12 09:00:09 Hail: INFO: Finished type imputation
  Loading column 'FID' as type 'int32' (imputed)
  Loading column 'IID' as type 'str' (imputed)
  Loading column 'Height' as type 'float64' (imputed)

hail table¶

ここで読み込んだデータは、matrixtableではなく、tableと呼ばれる形式で扱われます。 hailには、matrixtableとtableがあります。

In [11]:

height

Out[11]:

<hail.table.Table at 0x7fed76a40610>

summaraize¶

In [12]:

height.summarize()

2020-06-12 09:00:09 Hail: INFO: Coerced sorted dataset

503 records.

FID (int32):

Non-missing 503 (100.00%)

Missing 0

Minimum 0

Maximum 0

Mean 0.00

Std Dev 0.00
IID (str):

Non-missing 503 (100.00%)

Missing 0

Min Size 7

Max Size 7

Mean Size 7.00

Sample Values ['HG00096', 'HG00097', 'HG00099', 'HG00100', 'HG00101']
Height (float64):

Non-missing 503 (100.00%)

Missing 0

Minimum 154.30

Maximum 181.79

Mean 164.92

Std Dev 4.47

filter¶

tableでもfilterなどの操作が可能で、例えば次のようにすると身長が180cm以上のデータに絞り込むことができます。

In [13]:

height.filter(height.Height > 180 ).show()

FID	IID	Height
int32	str	float64
0	"HG00154"	1.82e+02
0	"NA12890"	1.80e+02

annotate_cols¶

annotate_colsを使うことで先に読み込んだplink_dataにアノテーションをつけることができます。以下の例ではplink_data.sを使いheight.Heightの値を呼び出して、そのデータをHeightという名前で追加しています。

In [14]:

plink_data2 = plink_data.annotate_cols(Height = height[plink_data.s].Height)

In [15]:

plink_data2.cols().show()

2020-06-12 09:00:21 Hail: INFO: Ordering unsorted dataset with network shuffle
2020-06-12 09:00:22 Hail: INFO: Coerced sorted dataset

s	fam_id	pat_id	mat_id	is_female	quant_pheno	Height
str	str	str	str	bool	float64	float64
"HG00096"	NA	NA	NA	false	NA	1.67e+02
"HG00097"	NA	NA	NA	true	NA	1.70e+02
"HG00099"	NA	NA	NA	true	NA	1.71e+02
"HG00100"	NA	NA	NA	true	NA	1.67e+02
"HG00101"	NA	NA	NA	false	NA	1.65e+02
"HG00102"	NA	NA	NA	true	NA	1.70e+02
"HG00103"	NA	NA	NA	false	NA	1.63e+02
"HG00105"	NA	NA	NA	false	NA	1.65e+02
"HG00106"	NA	NA	NA	true	NA	1.74e+02
"HG00107"	NA	NA	NA	false	NA	1.59e+02

showing top 10 rows

今回は次のような操作を行いました。

plink形式(bed,bim,fam)を読み込み
任意のテキストファイルを読み込み
読み込んだplink形式だったデータに任意のテキストファイルにあるアノテーションデータを追加

今回は小さなデータでしたが、Hailの特徴はとても大きなデータを扱えることにあります。

次回の記事ではHailのtutorialを参考にして、大きなデータを扱うと実際どのように動くのかをお示しできればと思います。

Nabe International | Takeru Boost >> Takeruを回せ！ >>

大規模GWASやPRS研究を可能にするHail(3)

Takeru Elysia news!¶