このノートブックでは、 染色体ごとに保存されているゲノムデータ (VCF形式)をマージし、matrix table 形式で保存します。
その後、imputation quality の低いバリアント (DR2 < 0.3) を除外します。
またその後のゲノムデータを用いて、PRS の計算手順を説明します。
2種類の脳梗塞 PRS モデルを事例として用います。
export PYSPARK_SUBMIT_ARGS='--driver-memory 48g --executor-memory 48g pyspark-shell'
のようにこの計算用のメモリを大きく確保してください。(この場合48ギガのメモリを割り当てています。) これを行っていただかないと、いずれかのセルで OutOfMemory エラーが発生し、それ以降のセル実行は機能しなくなります。下記のコードを実行してください。 ページ上側のメニューバーにある 実行 ボタンを押下することで、実行することができます。
import hail as hl
hl.init()
from hail.plot import show
from pprint import pprint
hl.plot.output_notebook()
2022-11-23 12:14:48.075 WARN NativeCodeLoader:60 - Unable to load native-hadoop library for your platform... using builtin-java classes where applicable
Setting default log level to "WARN". To adjust logging level use sc.setLogLevel(newLevel). For SparkR, use setLogLevel(newLevel). Running on Apache Spark version 3.1.3 SparkUI available at http://guaca021:4040 Welcome to __ __ <>__ / /_/ /__ __/ / / __ / _ `/ / / /_/ /_/\_,_/_/_/ version 0.2.105-acd89e80c345 LOGGING: writing to /lustre8/home/kozonishida-pg/prs-on-hail/hail-20221123-1214-0.2.105-acd89e80c345.log
!ls outputs/
chr1.beagle.log chr20.conform-gt.vcf.gz chr1.beagle.vcf.gz chr20.mt chr1.beagle.vcf.gz.tbi chr21.beagle.log chr1.conform-gt.log chr21.beagle.vcf.gz chr1.conform-gt.vcf.gz chr21.beagle.vcf.gz.tbi chr1.mt chr21.conform-gt.log chr10.beagle.log chr21.conform-gt.vcf.gz chr10.beagle.vcf.gz chr21.mt chr10.beagle.vcf.gz.tbi chr22.beagle.log chr10.conform-gt.log chr22.beagle.vcf.gz chr10.conform-gt.vcf.gz chr22.beagle.vcf.gz.tbi chr10.mt chr22.conform-gt.log chr11.beagle.log chr22.conform-gt.vcf.gz chr11.beagle.vcf.gz chr22.mt chr11.beagle.vcf.gz.tbi chr3.beagle.log chr11.conform-gt.log chr3.beagle.vcf.gz chr11.conform-gt.vcf.gz chr3.beagle.vcf.gz.tbi chr11.mt chr3.conform-gt.log chr12.beagle.log chr3.conform-gt.vcf.gz chr12.beagle.vcf.gz chr3.mt chr12.beagle.vcf.gz.tbi chr4.beagle.log chr12.conform-gt.log chr4.beagle.vcf.gz chr12.conform-gt.vcf.gz chr4.beagle.vcf.gz.tbi chr12.mt chr4.conform-gt.log chr13.beagle.log chr4.conform-gt.vcf.gz chr13.beagle.vcf.gz chr4.mt chr13.beagle.vcf.gz.tbi chr5.beagle.log chr13.conform-gt.log chr5.beagle.vcf.gz chr13.conform-gt.vcf.gz chr5.beagle.vcf.gz.tbi chr13.mt chr5.conform-gt.log chr14.beagle.log chr5.conform-gt.vcf.gz chr14.beagle.vcf.gz chr5.mt chr14.beagle.vcf.gz.tbi chr6.beagle.log chr14.conform-gt.log chr6.beagle.vcf.gz chr14.conform-gt.vcf.gz chr6.beagle.vcf.gz.tbi chr14.mt chr6.conform-gt.log chr15.beagle.log chr6.conform-gt.vcf.gz chr15.beagle.vcf.gz chr6.mt chr15.beagle.vcf.gz.tbi chr7.beagle.log chr15.conform-gt.log chr7.beagle.vcf.gz chr15.conform-gt.vcf.gz chr7.beagle.vcf.gz.tbi chr15.mt chr7.conform-gt.log chr16.beagle.log chr7.conform-gt.vcf.gz chr16.beagle.vcf.gz chr7.mt chr16.beagle.vcf.gz.tbi chr8.beagle.log chr16.conform-gt.log chr8.beagle.vcf.gz chr16.conform-gt.vcf.gz chr8.beagle.vcf.gz.tbi chr16.mt chr8.conform-gt.log chr17.beagle.log chr8.conform-gt.vcf.gz chr17.beagle.vcf.gz chr8.mt chr17.beagle.vcf.gz.tbi chr9.beagle.log chr17.conform-gt.log chr9.beagle.vcf.gz chr17.conform-gt.vcf.gz chr9.beagle.vcf.gz.tbi chr17.mt chr9.conform-gt.log chr18.beagle.log chr9.conform-gt.vcf.gz chr18.beagle.vcf.gz chr9.mt chr18.beagle.vcf.gz.tbi chrAll.filtered.matched_PGS002724.mt chr18.conform-gt.log chrAll.filtered.matched_PGS002725.mt chr18.conform-gt.vcf.gz chrAll.filtered.mt chr18.mt chrAll.mt chr19.beagle.log chrX_PAR1.beagle.log chr19.beagle.vcf.gz chrX_PAR1.beagle.vcf.gz chr19.beagle.vcf.gz.tbi chrX_PAR1.beagle.vcf.gz.tbi chr19.conform-gt.log chrX_PAR1.conform-gt.log chr19.conform-gt.vcf.gz chrX_PAR1.conform-gt.vcf.gz chr19.mt chrX_PAR2.beagle.log chr2.beagle.log chrX_PAR2.beagle.vcf.gz chr2.beagle.vcf.gz chrX_PAR2.beagle.vcf.gz.tbi chr2.beagle.vcf.gz.tbi chrX_PAR2.conform-gt.log chr2.conform-gt.log chrX_PAR2.conform-gt.vcf.gz chr2.conform-gt.vcf.gz chrX_nonPAR.beagle.log chr2.mt chrX_nonPAR.beagle.vcf.gz chr20.beagle.log chrX_nonPAR.beagle.vcf.gz.tbi chr20.beagle.vcf.gz chrX_nonPAR.conform-gt.log chr20.beagle.vcf.gz.tbi chrX_nonPAR.conform-gt.vcf.gz chr20.conform-gt.log
ゲノムデータは、1番染色体から22番染色体まで、染色体ごとに異なるファイルに保存されています。
例えば、1番染色体のゲノムデータは、 outputs/chr1.beagle.vcf.gz
に保存されています。
また、2番染色体のゲノムデータは、 outputs/chr2.beagle.vcf.gz
に 22番染色体のゲノムデータは、 outputs/chr22.beagle.vcf.gz
に保存されています。
ゲノムデータは、よく利用されるファイル形式である VCF フォーマットで保存されています。 これを、matrix table 形式に変換して、保存します。
下記コマンドを実行してください。 このセルの実行にはかなりの待ち時間が生じますが「このノートブックを実行する前に下記をご注意ください」のように途中でセル実行を停止することはおすすめしません。
for chr in range(1,23):
infile = 'outputs/chr' + str(chr) + '.beagle.vcf.gz'
outfile = 'outputs/chr' + str(chr) + '.mt'
print(chr, infile, outfile)
hl.import_vcf(infile, force_bgz=True).write(outfile, overwrite=True)
1 outputs/chr1.beagle.vcf.gz outputs/chr1.mt
2022-11-24 00:12:03.759 Hail: INFO: scanning VCF for sortedness... 2022-11-24 00:12:36.902 Hail: INFO: Coerced sorted VCF - no additional import work to do 2022-11-24 00:26:13.044 Hail: INFO: wrote matrix table with 2428653 rows and 2318 columns in 7 partitions to outputs/chr1.mt
2 outputs/chr2.beagle.vcf.gz outputs/chr2.mt
2022-11-24 00:26:13.449 Hail: INFO: scanning VCF for sortedness... 2022-11-24 00:26:45.883 Hail: INFO: Coerced sorted VCF - no additional import work to do 2022-11-24 00:41:24.089 Hail: INFO: wrote matrix table with 2627240 rows and 2318 columns in 7 partitions to outputs/chr2.mt
3 outputs/chr3.beagle.vcf.gz outputs/chr3.mt
2022-11-24 00:41:24.562 Hail: INFO: scanning VCF for sortedness... 2022-11-24 00:41:50.082 Hail: INFO: Coerced sorted VCF - no additional import work to do 2022-11-24 00:55:55.758 Hail: INFO: wrote matrix table with 2186425 rows and 2318 columns in 6 partitions to outputs/chr3.mt
4 outputs/chr4.beagle.vcf.gz outputs/chr4.mt
2022-11-24 00:55:56.162 Hail: INFO: scanning VCF for sortedness... 2022-11-24 00:56:25.917 Hail: INFO: Coerced sorted VCF - no additional import work to do 2022-11-24 01:08:23.150 Hail: INFO: wrote matrix table with 2212857 rows and 2318 columns in 7 partitions to outputs/chr4.mt
5 outputs/chr5.beagle.vcf.gz outputs/chr5.mt
2022-11-24 01:08:23.583 Hail: INFO: scanning VCF for sortedness... 2022-11-24 01:08:50.201 Hail: INFO: Coerced sorted VCF - no additional import work to do 2022-11-24 01:21:52.499 Hail: INFO: wrote matrix table with 1986979 rows and 2318 columns in 6 partitions to outputs/chr5.mt
6 outputs/chr6.beagle.vcf.gz outputs/chr6.mt
2022-11-24 01:21:52.956 Hail: INFO: scanning VCF for sortedness... 2022-11-24 01:22:18.976 Hail: INFO: Coerced sorted VCF - no additional import work to do 2022-11-24 01:35:31.727 Hail: INFO: wrote matrix table with 1964598 rows and 2318 columns in 6 partitions to outputs/chr6.mt
7 outputs/chr7.beagle.vcf.gz outputs/chr7.mt
2022-11-24 01:35:32.141 Hail: INFO: scanning VCF for sortedness... 2022-11-24 01:35:58.592 Hail: INFO: Coerced sorted VCF - no additional import work to do 2022-11-24 01:50:01.298 Hail: INFO: wrote matrix table with 1801231 rows and 2318 columns in 5 partitions to outputs/chr7.mt
8 outputs/chr8.beagle.vcf.gz outputs/chr8.mt
2022-11-24 01:50:01.775 Hail: INFO: scanning VCF for sortedness... 2022-11-24 01:50:31.947 Hail: INFO: Coerced sorted VCF - no additional import work to do 2022-11-24 02:04:55.201 Hail: INFO: wrote matrix table with 1722793 rows and 2318 columns in 5 partitions to outputs/chr8.mt
9 outputs/chr9.beagle.vcf.gz outputs/chr9.mt
2022-11-24 02:04:55.609 Hail: INFO: scanning VCF for sortedness... 2022-11-24 02:05:19.902 Hail: INFO: Coerced sorted VCF - no additional import work to do 2022-11-24 02:18:18.739 Hail: INFO: wrote matrix table with 1342561 rows and 2318 columns in 4 partitions to outputs/chr9.mt
10 outputs/chr10.beagle.vcf.gz outputs/chr10.mt
2022-11-24 02:18:19.246 Hail: INFO: scanning VCF for sortedness... 2022-11-24 02:18:45.990 Hail: INFO: Coerced sorted VCF - no additional import work to do 2022-11-24 02:34:04.301 Hail: INFO: wrote matrix table with 1532460 rows and 2318 columns in 4 partitions to outputs/chr10.mt
11 outputs/chr11.beagle.vcf.gz outputs/chr11.mt
2022-11-24 02:34:04.802 Hail: INFO: scanning VCF for sortedness... 2022-11-24 02:34:30.324 Hail: INFO: Coerced sorted VCF - no additional import work to do 2022-11-24 02:48:51.586 Hail: INFO: wrote matrix table with 1520309 rows and 2318 columns in 4 partitions to outputs/chr11.mt
12 outputs/chr12.beagle.vcf.gz outputs/chr12.mt
2022-11-24 02:48:51.976 Hail: INFO: scanning VCF for sortedness... 2022-11-24 02:49:26.205 Hail: INFO: Coerced sorted VCF - no additional import work to do 2022-11-24 03:03:38.788 Hail: INFO: wrote matrix table with 1467858 rows and 2318 columns in 4 partitions to outputs/chr12.mt
13 outputs/chr13.beagle.vcf.gz outputs/chr13.mt
2022-11-24 03:03:39.549 Hail: INFO: scanning VCF for sortedness... 2022-11-24 03:04:03.862 Hail: INFO: Coerced sorted VCF - no additional import work to do 2022-11-24 03:18:20.910 Hail: INFO: wrote matrix table with 1099285 rows and 2318 columns in 3 partitions to outputs/chr13.mt
14 outputs/chr14.beagle.vcf.gz outputs/chr14.mt
2022-11-24 03:18:21.432 Hail: INFO: scanning VCF for sortedness... 2022-11-24 03:18:56.995 Hail: INFO: Coerced sorted VCF - no additional import work to do 2022-11-24 03:32:45.903 Hail: INFO: wrote matrix table with 1002655 rows and 2318 columns in 3 partitions to outputs/chr14.mt
15 outputs/chr15.beagle.vcf.gz outputs/chr15.mt
2022-11-24 03:32:46.434 Hail: INFO: scanning VCF for sortedness... 2022-11-24 03:33:15.573 Hail: INFO: Coerced sorted VCF - no additional import work to do 2022-11-24 03:46:02.317 Hail: INFO: wrote matrix table with 911380 rows and 2318 columns in 3 partitions to outputs/chr15.mt
16 outputs/chr16.beagle.vcf.gz outputs/chr16.mt
2022-11-24 03:46:02.789 Hail: INFO: scanning VCF for sortedness... 2022-11-24 03:46:30.289 Hail: INFO: Coerced sorted VCF - no additional import work to do 2022-11-24 03:59:20.950 Hail: INFO: wrote matrix table with 988654 rows and 2318 columns in 3 partitions to outputs/chr16.mt
17 outputs/chr17.beagle.vcf.gz outputs/chr17.mt
2022-11-24 03:59:21.566 Hail: INFO: scanning VCF for sortedness... 2022-11-24 03:59:40.999 Hail: INFO: Coerced sorted VCF - no additional import work to do 2022-11-24 04:11:07.839 Hail: INFO: wrote matrix table with 864311 rows and 2318 columns in 3 partitions to outputs/chr17.mt
18 outputs/chr18.beagle.vcf.gz outputs/chr18.mt
2022-11-24 04:11:08.461 Hail: INFO: scanning VCF for sortedness... 2022-11-24 04:11:29.928 Hail: INFO: Coerced sorted VCF - no additional import work to do 2022-11-24 04:22:58.193 Hail: INFO: wrote matrix table with 864327 rows and 2318 columns in 3 partitions to outputs/chr18.mt
19 outputs/chr19.beagle.vcf.gz outputs/chr19.mt
2022-11-24 04:22:58.737 Hail: INFO: scanning VCF for sortedness... 2022-11-24 04:23:31.187 Hail: INFO: Coerced sorted VCF - no additional import work to do 2022-11-24 04:37:33.986 Hail: INFO: wrote matrix table with 706126 rows and 2318 columns in 2 partitions to outputs/chr19.mt
20 outputs/chr20.beagle.vcf.gz outputs/chr20.mt
2022-11-24 04:37:34.601 Hail: INFO: scanning VCF for sortedness... 2022-11-24 04:37:55.931 Hail: INFO: Coerced sorted VCF - no additional import work to do 2022-11-24 04:51:29.869 Hail: INFO: wrote matrix table with 679241 rows and 2318 columns in 2 partitions to outputs/chr20.mt
21 outputs/chr21.beagle.vcf.gz outputs/chr21.mt
2022-11-24 04:51:30.363 Hail: INFO: scanning VCF for sortedness... 2022-11-24 04:52:11.583 Hail: INFO: Coerced sorted VCF - no additional import work to do 2022-11-24 05:09:18.461 Hail: INFO: wrote matrix table with 427409 rows and 2318 columns in 1 partition to outputs/chr21.mt
22 outputs/chr22.beagle.vcf.gz outputs/chr22.mt
2022-11-24 05:09:18.993 Hail: INFO: scanning VCF for sortedness... 2022-11-24 05:09:54.515 Hail: INFO: Coerced sorted VCF - no additional import work to do 2022-11-24 05:26:20.144 Hail: INFO: wrote matrix table with 424147 rows and 2318 columns in 1 partition to outputs/chr22.mt
続いて、1番染色体から22番染色体の matrix table 形式のファイルを、ひとつの matrix table 形式のファイルにマージします。
下記コマンドを実行してください。
chr = 1
file = 'outputs/chr' + str(chr) + '.mt'
mt = hl.read_matrix_table(file)
print(chr, file, mt.count(), mt.count())
for chr in range(2,23):
file = 'outputs/chr' + str(chr) + '.mt'
tmpmt = hl.read_matrix_table(file)
mt = mt.union_rows(tmpmt)
print(chr, file, tmpmt.count(), mt.count())
1 outputs/chr1.mt (2428653, 2318) (2428653, 2318) 2 outputs/chr2.mt (2627240, 2318) (5055893, 2318) 3 outputs/chr3.mt (2186425, 2318) (7242318, 2318) 4 outputs/chr4.mt (2212857, 2318) (9455175, 2318) 5 outputs/chr5.mt (1986979, 2318) (11442154, 2318) 6 outputs/chr6.mt (1964598, 2318) (13406752, 2318) 7 outputs/chr7.mt (1801231, 2318) (15207983, 2318) 8 outputs/chr8.mt (1722793, 2318) (16930776, 2318) 9 outputs/chr9.mt (1342561, 2318) (18273337, 2318) 10 outputs/chr10.mt (1532460, 2318) (19805797, 2318) 11 outputs/chr11.mt (1520309, 2318) (21326106, 2318) 12 outputs/chr12.mt (1467858, 2318) (22793964, 2318) 13 outputs/chr13.mt (1099285, 2318) (23893249, 2318) 14 outputs/chr14.mt (1002655, 2318) (24895904, 2318) 15 outputs/chr15.mt (911380, 2318) (25807284, 2318) 16 outputs/chr16.mt (988654, 2318) (26795938, 2318) 17 outputs/chr17.mt (864311, 2318) (27660249, 2318) 18 outputs/chr18.mt (864327, 2318) (28524576, 2318) 19 outputs/chr19.mt (706126, 2318) (29230702, 2318) 20 outputs/chr20.mt (679241, 2318) (29909943, 2318) 21 outputs/chr21.mt (427409, 2318) (30337352, 2318) 22 outputs/chr22.mt (424147, 2318) (30761499, 2318)
ゲノムデータのマージが完了しました。 マージ後のゲノムデータのバリアント数を表示します。
mt.count()
(30761499, 2318)
これは、次のことを意味します。
今後のステップの実行時間を短縮するため、マージしたゲノムデータを保存し、再度読み込みます。
下記のコードは、マージしたゲノムデータを outputs/chrAll.mt に保存します。
mt.write('outputs/chrAll.mt', overwrite=True)
2022-11-15 20:03:09.905 Hail: INFO: wrote matrix table with 30761499 rows and 2318 columns in 89 partitions to outputs/chrAll.mt Total size: 18.99 GiB * Rows/entries: 18.99 GiB * Columns: 9.09 KiB * Globals: 11.00 B * Smallest partition: 282779 rows (181.07 MiB) * Largest partition: 424147 rows (311.72 MiB)
下記のコードは、 outputs/chrAll.mt
を読み込みます。
mt = hl.read_matrix_table('outputs/chrAll.mt')
下記のコードは、ゲノムデータのバリアント情報に variantID
を追加します。
mt = mt.annotate_rows(variantID = (hl.str(mt.locus.contig) + ":" + hl.str(mt.locus.position)) )
下記のコードは、 variantID
を追加した後のバリアント情報を表示します。
show(5)
は、先頭の 5 個のバリアントのみを表示する、ことを意味します。
mt.rows().show(5)
[Stage 0:> (0 + 4) / 4]
showing top 5 rows
variantID
のカラムが追加されていることが分かります。
mt には multi-allelic site が含まれます。次にその割合を確認してみます。
mt1 = mt.filter_rows(hl.len(mt.info.DR2) == 1)
mt1.count()
[Stage 2:========================================================>(88 + 1) / 89]
(30361461, 2318)
mtnot1 = mt.filter_rows(hl.len(mt.info.DR2) > 1)
mtnot1.count()
[Stage 3:========================================================>(88 + 1) / 89]
(400038, 2318)
1% ほどが multi-allelic site であることがわかります。 1% は無視するにはやや多いですが、
このチュートリアルでは内容をわかりやすくするために除外したもの(すなわちmt1
)を今後用います。
下記のコードでは、各バリアントの imputation quality
の分布を表示します。
p = hl.plot.histogram(mt1.info.DR2.first(), title='Imputation Quality Histogram', legend='Imputation Quality (DR2)')
show(p)
[Stage 5:=======================================================> (86 + 3) / 89]
imputation quality
が低いバリアントが少しあることが分かります。
下記のコードは、imputation quality
が低い(DR2 < 0.3)バリアントを除外します。
mt1_filt = mt1.filter_rows(mt1.info.DR2.first()>=0.3)
下記のコードは、imputation quality
が低いバリアントを除外した後の分布を表示します。
p = hl.plot.histogram(mt1_filt.info.DR2.first(), title='Imputation Quality Histogram', legend='Imputation Quality (DR2)')
show(p)
[Stage 15:======================================================> (87 + 2) / 89]
imputation quality
が低いバリアントが除外されたことが分かります。
下記のコードは、imputation quality
が低いバリアントを除外した後のバリアントの個数を表示します。
mt1_filt.count()
[Stage 10:=======================================================>(88 + 1) / 89]
(29799034, 2318)
(29799034, 2318) と表示されました。
これは、次のことを意味します。
研究対象者の人数が 2318 名
バリアントの個数が 29799034 個
下記のコードは、imputation quality
が低いバリアントを除外した後の allele frequency の分布を表示します。
p = hl.plot.histogram(mt1_filt.info.AF.first(), title='AF Histogram', legend='AF', bins=50)
show(p)
[Stage 12:======================================================> (87 + 2) / 89]
AF<1% のバリアントが多くあることが分かります。
下記のコードは、マージしたゲノムデータを outputs/chrAll.filtered.mt
に保存します。
ここまでのセルには実行完了までに時間がかかるものが存在します。
ここをチェックポイントとして同じことに時間をかけずに再開できるようにしておきましょう。
mt1_filt.write('outputs/chrAll.filtered.mt', overwrite=True)
2022-11-17 18:33:43.273 Hail: INFO: wrote matrix table with 29799034 rows and 2318 columns in 89 partitions to outputs/chrAll.filtered.mt
下記のコードでは、Step 7 で保存したゲノムデータを読み込みます。
mt1_filt = hl.read_matrix_table('outputs/chrAll.filtered.mt')
!ls -alh outputs
合計 14G drwxr-xr-x 26 kozonishida-pg oo-nig-pg 12K 11月 17 14:14 . drwxr-xr-x 5 kozonishida-pg oo-nig-pg 4.0K 11月 17 16:56 .. -rw-r--r-- 1 kozonishida-pg oo-nig-pg 6.7K 11月 14 14:44 chr1.beagle.log -rw-r--r-- 1 kozonishida-pg oo-nig-pg 914M 11月 14 14:44 chr1.beagle.vcf.gz -rw-r--r-- 1 kozonishida-pg oo-nig-pg 207K 11月 14 14:44 chr1.beagle.vcf.gz.tbi -rw-r--r-- 1 kozonishida-pg oo-nig-pg 16M 11月 14 14:46 chr1.conform-gt.log -rw-r--r-- 1 kozonishida-pg oo-nig-pg 88M 11月 14 14:46 chr1.conform-gt.vcf.gz drwxr-xr-x 8 kozonishida-pg oo-nig-pg 4.0K 11月 14 15:18 chr1.mt -rw-r--r-- 1 kozonishida-pg oo-nig-pg 4.5K 11月 14 14:44 chr10.beagle.log -rw-r--r-- 1 kozonishida-pg oo-nig-pg 574M 11月 14 14:45 chr10.beagle.vcf.gz -rw-r--r-- 1 kozonishida-pg oo-nig-pg 121K 11月 14 14:44 chr10.beagle.vcf.gz.tbi -rw-r--r-- 1 kozonishida-pg oo-nig-pg 10M 11月 14 14:46 chr10.conform-gt.log -rw-r--r-- 1 kozonishida-pg oo-nig-pg 59M 11月 14 14:46 chr10.conform-gt.vcf.gz drwxr-xr-x 8 kozonishida-pg oo-nig-pg 4.0K 11月 14 17:08 chr10.mt -rw-r--r-- 1 kozonishida-pg oo-nig-pg 4.5K 11月 14 14:44 chr11.beagle.log -rw-r--r-- 1 kozonishida-pg oo-nig-pg 554M 11月 14 14:45 chr11.beagle.vcf.gz -rw-r--r-- 1 kozonishida-pg oo-nig-pg 121K 11月 14 14:44 chr11.beagle.vcf.gz.tbi -rw-r--r-- 1 kozonishida-pg oo-nig-pg 9.7M 11月 14 14:46 chr11.conform-gt.log -rw-r--r-- 1 kozonishida-pg oo-nig-pg 57M 11月 14 14:46 chr11.conform-gt.vcf.gz drwxr-xr-x 8 kozonishida-pg oo-nig-pg 4.0K 11月 14 17:20 chr11.mt -rw-r--r-- 1 kozonishida-pg oo-nig-pg 4.6K 11月 14 14:44 chr12.beagle.log -rw-r--r-- 1 kozonishida-pg oo-nig-pg 546M 11月 14 14:45 chr12.beagle.vcf.gz -rw-r--r-- 1 kozonishida-pg oo-nig-pg 121K 11月 14 14:44 chr12.beagle.vcf.gz.tbi -rw-r--r-- 1 kozonishida-pg oo-nig-pg 9.4M 11月 14 14:46 chr12.conform-gt.log -rw-r--r-- 1 kozonishida-pg oo-nig-pg 55M 11月 14 14:46 chr12.conform-gt.vcf.gz drwxr-xr-x 8 kozonishida-pg oo-nig-pg 4.0K 11月 14 17:32 chr12.mt -rw-r--r-- 1 kozonishida-pg oo-nig-pg 3.8K 11月 14 14:44 chr13.beagle.log -rw-r--r-- 1 kozonishida-pg oo-nig-pg 405M 11月 14 14:45 chr13.beagle.vcf.gz -rw-r--r-- 1 kozonishida-pg oo-nig-pg 88K 11月 14 14:44 chr13.beagle.vcf.gz.tbi -rw-r--r-- 1 kozonishida-pg oo-nig-pg 6.9M 11月 14 14:46 chr13.conform-gt.log -rw-r--r-- 1 kozonishida-pg oo-nig-pg 41M 11月 14 14:46 chr13.conform-gt.vcf.gz drwxr-xr-x 8 kozonishida-pg oo-nig-pg 4.0K 11月 14 17:43 chr13.mt -rw-r--r-- 1 kozonishida-pg oo-nig-pg 3.9K 11月 14 14:44 chr14.beagle.log -rw-r--r-- 1 kozonishida-pg oo-nig-pg 387M 11月 14 14:45 chr14.beagle.vcf.gz -rw-r--r-- 1 kozonishida-pg oo-nig-pg 82K 11月 14 14:44 chr14.beagle.vcf.gz.tbi -rw-r--r-- 1 kozonishida-pg oo-nig-pg 6.4M 11月 14 14:46 chr14.conform-gt.log -rw-r--r-- 1 kozonishida-pg oo-nig-pg 38M 11月 14 14:46 chr14.conform-gt.vcf.gz drwxr-xr-x 8 kozonishida-pg oo-nig-pg 4.0K 11月 14 17:54 chr14.mt -rw-r--r-- 1 kozonishida-pg oo-nig-pg 3.9K 11月 14 14:44 chr15.beagle.log -rw-r--r-- 1 kozonishida-pg oo-nig-pg 347M 11月 14 14:45 chr15.beagle.vcf.gz -rw-r--r-- 1 kozonishida-pg oo-nig-pg 76K 11月 14 14:44 chr15.beagle.vcf.gz.tbi -rw-r--r-- 1 kozonishida-pg oo-nig-pg 6.0M 11月 14 14:46 chr15.conform-gt.log -rw-r--r-- 1 kozonishida-pg oo-nig-pg 36M 11月 14 14:46 chr15.conform-gt.vcf.gz drwxr-xr-x 8 kozonishida-pg oo-nig-pg 4.0K 11月 14 18:04 chr15.mt -rw-r--r-- 1 kozonishida-pg oo-nig-pg 3.8K 11月 14 14:44 chr16.beagle.log -rw-r--r-- 1 kozonishida-pg oo-nig-pg 374M 11月 14 14:45 chr16.beagle.vcf.gz -rw-r--r-- 1 kozonishida-pg oo-nig-pg 74K 11月 14 14:44 chr16.beagle.vcf.gz.tbi -rw-r--r-- 1 kozonishida-pg oo-nig-pg 6.5M 11月 14 14:46 chr16.conform-gt.log -rw-r--r-- 1 kozonishida-pg oo-nig-pg 40M 11月 14 14:46 chr16.conform-gt.vcf.gz drwxr-xr-x 8 kozonishida-pg oo-nig-pg 4.0K 11月 14 18:14 chr16.mt -rw-r--r-- 1 kozonishida-pg oo-nig-pg 3.9K 11月 14 14:44 chr17.beagle.log -rw-r--r-- 1 kozonishida-pg oo-nig-pg 327M 11月 14 14:45 chr17.beagle.vcf.gz -rw-r--r-- 1 kozonishida-pg oo-nig-pg 72K 11月 14 14:44 chr17.beagle.vcf.gz.tbi -rw-r--r-- 1 kozonishida-pg oo-nig-pg 5.7M 11月 14 14:46 chr17.conform-gt.log -rw-r--r-- 1 kozonishida-pg oo-nig-pg 34M 11月 14 14:46 chr17.conform-gt.vcf.gz drwxr-xr-x 8 kozonishida-pg oo-nig-pg 4.0K 11月 14 18:23 chr17.mt -rw-r--r-- 1 kozonishida-pg oo-nig-pg 3.9K 11月 14 14:44 chr18.beagle.log -rw-r--r-- 1 kozonishida-pg oo-nig-pg 322M 11月 14 14:45 chr18.beagle.vcf.gz -rw-r--r-- 1 kozonishida-pg oo-nig-pg 70K 11月 14 14:44 chr18.beagle.vcf.gz.tbi -rw-r--r-- 1 kozonishida-pg oo-nig-pg 5.7M 11月 14 14:46 chr18.conform-gt.log -rw-r--r-- 1 kozonishida-pg oo-nig-pg 35M 11月 14 14:46 chr18.conform-gt.vcf.gz drwxr-xr-x 8 kozonishida-pg oo-nig-pg 4.0K 11月 14 18:32 chr18.mt -rw-r--r-- 1 kozonishida-pg oo-nig-pg 3.2K 11月 14 14:44 chr19.beagle.log -rw-r--r-- 1 kozonishida-pg oo-nig-pg 284M 11月 14 14:45 chr19.beagle.vcf.gz -rw-r--r-- 1 kozonishida-pg oo-nig-pg 51K 11月 14 14:44 chr19.beagle.vcf.gz.tbi -rw-r--r-- 1 kozonishida-pg oo-nig-pg 4.2M 11月 14 14:46 chr19.conform-gt.log -rw-r--r-- 1 kozonishida-pg oo-nig-pg 26M 11月 14 14:46 chr19.conform-gt.vcf.gz drwxr-xr-x 8 kozonishida-pg oo-nig-pg 4.0K 11月 14 18:43 chr19.mt -rw-r--r-- 1 kozonishida-pg oo-nig-pg 6.6K 11月 14 14:44 chr2.beagle.log -rw-r--r-- 1 kozonishida-pg oo-nig-pg 957M 11月 14 14:44 chr2.beagle.vcf.gz -rw-r--r-- 1 kozonishida-pg oo-nig-pg 219K 11月 14 14:44 chr2.beagle.vcf.gz.tbi -rw-r--r-- 1 kozonishida-pg oo-nig-pg 16M 11月 14 14:46 chr2.conform-gt.log -rw-r--r-- 1 kozonishida-pg oo-nig-pg 93M 11月 14 14:46 chr2.conform-gt.vcf.gz drwxr-xr-x 8 kozonishida-pg oo-nig-pg 4.0K 11月 14 15:32 chr2.mt -rw-r--r-- 1 kozonishida-pg oo-nig-pg 3.1K 11月 14 14:44 chr20.beagle.log -rw-r--r-- 1 kozonishida-pg oo-nig-pg 249M 11月 14 14:45 chr20.beagle.vcf.gz -rw-r--r-- 1 kozonishida-pg oo-nig-pg 55K 11月 14 14:44 chr20.beagle.vcf.gz.tbi -rw-r--r-- 1 kozonishida-pg oo-nig-pg 4.7M 11月 14 14:46 chr20.conform-gt.log -rw-r--r-- 1 kozonishida-pg oo-nig-pg 29M 11月 14 14:46 chr20.conform-gt.vcf.gz drwxr-xr-x 8 kozonishida-pg oo-nig-pg 4.0K 11月 14 18:54 chr20.mt -rw-r--r-- 1 kozonishida-pg oo-nig-pg 2.4K 11月 14 14:44 chr21.beagle.log -rw-r--r-- 1 kozonishida-pg oo-nig-pg 174M 11月 14 14:45 chr21.beagle.vcf.gz -rw-r--r-- 1 kozonishida-pg oo-nig-pg 33K 11月 14 14:44 chr21.beagle.vcf.gz.tbi -rw-r--r-- 1 kozonishida-pg oo-nig-pg 2.7M 11月 14 14:46 chr21.conform-gt.log -rw-r--r-- 1 kozonishida-pg oo-nig-pg 17M 11月 14 14:46 chr21.conform-gt.vcf.gz drwxr-xr-x 8 kozonishida-pg oo-nig-pg 4.0K 11月 14 19:06 chr21.mt -rw-r--r-- 1 kozonishida-pg oo-nig-pg 2.4K 11月 14 14:44 chr22.beagle.log -rw-r--r-- 1 kozonishida-pg oo-nig-pg 185M 11月 14 14:45 chr22.beagle.vcf.gz -rw-r--r-- 1 kozonishida-pg oo-nig-pg 33K 11月 14 14:44 chr22.beagle.vcf.gz.tbi -rw-r--r-- 1 kozonishida-pg oo-nig-pg 2.9M 11月 14 14:46 chr22.conform-gt.log -rw-r--r-- 1 kozonishida-pg oo-nig-pg 19M 11月 14 14:46 chr22.conform-gt.vcf.gz drwxr-xr-x 8 kozonishida-pg oo-nig-pg 4.0K 11月 14 19:19 chr22.mt -rw-r--r-- 1 kozonishida-pg oo-nig-pg 5.2K 11月 14 14:44 chr3.beagle.log -rw-r--r-- 1 kozonishida-pg oo-nig-pg 793M 11月 14 14:44 chr3.beagle.vcf.gz -rw-r--r-- 1 kozonishida-pg oo-nig-pg 180K 11月 14 14:44 chr3.beagle.vcf.gz.tbi -rw-r--r-- 1 kozonishida-pg oo-nig-pg 14M 11月 14 14:46 chr3.conform-gt.log -rw-r--r-- 1 kozonishida-pg oo-nig-pg 79M 11月 14 14:46 chr3.conform-gt.vcf.gz drwxr-xr-x 8 kozonishida-pg oo-nig-pg 4.0K 11月 14 15:44 chr3.mt -rw-r--r-- 1 kozonishida-pg oo-nig-pg 5.2K 11月 14 14:44 chr4.beagle.log -rw-r--r-- 1 kozonishida-pg oo-nig-pg 837M 11月 14 14:44 chr4.beagle.vcf.gz -rw-r--r-- 1 kozonishida-pg oo-nig-pg 174K 11月 14 14:44 chr4.beagle.vcf.gz.tbi -rw-r--r-- 1 kozonishida-pg oo-nig-pg 13M 11月 14 14:46 chr4.conform-gt.log -rw-r--r-- 1 kozonishida-pg oo-nig-pg 74M 11月 14 14:46 chr4.conform-gt.vcf.gz drwxr-xr-x 8 kozonishida-pg oo-nig-pg 4.0K 11月 14 15:55 chr4.mt -rw-r--r-- 1 kozonishida-pg oo-nig-pg 5.3K 11月 14 14:44 chr5.beagle.log -rw-r--r-- 1 kozonishida-pg oo-nig-pg 708M 11月 14 14:44 chr5.beagle.vcf.gz -rw-r--r-- 1 kozonishida-pg oo-nig-pg 163K 11月 14 14:44 chr5.beagle.vcf.gz.tbi -rw-r--r-- 1 kozonishida-pg oo-nig-pg 12M 11月 14 14:46 chr5.conform-gt.log -rw-r--r-- 1 kozonishida-pg oo-nig-pg 68M 11月 14 14:46 chr5.conform-gt.vcf.gz drwxr-xr-x 8 kozonishida-pg oo-nig-pg 4.0K 11月 14 16:09 chr5.mt -rw-r--r-- 1 kozonishida-pg oo-nig-pg 5.2K 11月 14 14:44 chr6.beagle.log -rw-r--r-- 1 kozonishida-pg oo-nig-pg 752M 11月 14 14:45 chr6.beagle.vcf.gz -rw-r--r-- 1 kozonishida-pg oo-nig-pg 155K 11月 14 14:44 chr6.beagle.vcf.gz.tbi -rw-r--r-- 1 kozonishida-pg oo-nig-pg 13M 11月 14 14:46 chr6.conform-gt.log -rw-r--r-- 1 kozonishida-pg oo-nig-pg 73M 11月 14 14:46 chr6.conform-gt.vcf.gz drwxr-xr-x 8 kozonishida-pg oo-nig-pg 4.0K 11月 14 16:20 chr6.mt -rw-r--r-- 1 kozonishida-pg oo-nig-pg 4.6K 11月 14 14:44 chr7.beagle.log -rw-r--r-- 1 kozonishida-pg oo-nig-pg 684M 11月 14 14:45 chr7.beagle.vcf.gz -rw-r--r-- 1 kozonishida-pg oo-nig-pg 144K 11月 14 14:44 chr7.beagle.vcf.gz.tbi -rw-r--r-- 1 kozonishida-pg oo-nig-pg 11M 11月 14 14:46 chr7.conform-gt.log -rw-r--r-- 1 kozonishida-pg oo-nig-pg 63M 11月 14 14:46 chr7.conform-gt.vcf.gz drwxr-xr-x 8 kozonishida-pg oo-nig-pg 4.0K 11月 14 16:32 chr7.mt -rw-r--r-- 1 kozonishida-pg oo-nig-pg 4.5K 11月 14 14:44 chr8.beagle.log -rw-r--r-- 1 kozonishida-pg oo-nig-pg 619M 11月 14 14:45 chr8.beagle.vcf.gz -rw-r--r-- 1 kozonishida-pg oo-nig-pg 132K 11月 14 14:44 chr8.beagle.vcf.gz.tbi -rw-r--r-- 1 kozonishida-pg oo-nig-pg 11M 11月 14 14:46 chr8.conform-gt.log -rw-r--r-- 1 kozonishida-pg oo-nig-pg 61M 11月 14 14:46 chr8.conform-gt.vcf.gz drwxr-xr-x 8 kozonishida-pg oo-nig-pg 4.0K 11月 14 16:46 chr8.mt -rw-r--r-- 1 kozonishida-pg oo-nig-pg 4.5K 11月 14 14:44 chr9.beagle.log -rw-r--r-- 1 kozonishida-pg oo-nig-pg 496M 11月 14 14:45 chr9.beagle.vcf.gz -rw-r--r-- 1 kozonishida-pg oo-nig-pg 110K 11月 14 14:44 chr9.beagle.vcf.gz.tbi -rw-r--r-- 1 kozonishida-pg oo-nig-pg 8.5M 11月 14 14:46 chr9.conform-gt.log -rw-r--r-- 1 kozonishida-pg oo-nig-pg 52M 11月 14 14:46 chr9.conform-gt.vcf.gz drwxr-xr-x 8 kozonishida-pg oo-nig-pg 4.0K 11月 14 16:56 chr9.mt drwxr-xr-x 8 kozonishida-pg oo-nig-pg 4.0K 11月 17 15:29 chrAll.filtered.mt drwxr-xr-x 8 kozonishida-pg oo-nig-pg 4.0K 11月 15 20:03 chrAll.mt -rw-r--r-- 1 kozonishida-pg oo-nig-pg 1.8K 11月 14 14:44 chrX_PAR1.beagle.log -rw-r--r-- 1 kozonishida-pg oo-nig-pg 43M 11月 14 14:45 chrX_PAR1.beagle.vcf.gz -rw-r--r-- 1 kozonishida-pg oo-nig-pg 1.8K 11月 14 14:44 chrX_PAR1.beagle.vcf.gz.tbi -rw-r--r-- 1 kozonishida-pg oo-nig-pg 123K 11月 14 14:46 chrX_PAR1.conform-gt.log -rw-r--r-- 1 kozonishida-pg oo-nig-pg 971K 11月 14 14:46 chrX_PAR1.conform-gt.vcf.gz -rw-r--r-- 1 kozonishida-pg oo-nig-pg 1.8K 11月 14 14:44 chrX_PAR2.beagle.log -rw-r--r-- 1 kozonishida-pg oo-nig-pg 2.1M 11月 14 14:46 chrX_PAR2.beagle.vcf.gz -rw-r--r-- 1 kozonishida-pg oo-nig-pg 697 11月 14 14:44 chrX_PAR2.beagle.vcf.gz.tbi -rw-r--r-- 1 kozonishida-pg oo-nig-pg 11K 11月 14 14:46 chrX_PAR2.conform-gt.log -rw-r--r-- 1 kozonishida-pg oo-nig-pg 71K 11月 14 14:46 chrX_PAR2.conform-gt.vcf.gz -rw-r--r-- 1 kozonishida-pg oo-nig-pg 4.7K 11月 14 14:44 chrX_nonPAR.beagle.log -rw-r--r-- 1 kozonishida-pg oo-nig-pg 479M 11月 14 14:46 chrX_nonPAR.beagle.vcf.gz -rw-r--r-- 1 kozonishida-pg oo-nig-pg 134K 11月 14 14:44 chrX_nonPAR.beagle.vcf.gz.tbi -rw-r--r-- 1 kozonishida-pg oo-nig-pg 4.3M 11月 14 14:46 chrX_nonPAR.conform-gt.log -rw-r--r-- 1 kozonishida-pg oo-nig-pg 27M 11月 14 14:46 chrX_nonPAR.conform-gt.vcf.gz
下記のコードを実行すると、prs-models/PGS002724.txt
と prs-models/PGS002724.txt
のデータが読み込まれます。
!wget https://ftp.ebi.ac.uk/pub/databases/spot/pgs/scores/PGS002724/ScoringFiles/PGS002724.txt.gz
--2022-11-17 14:37:13-- https://ftp.ebi.ac.uk/pub/databases/spot/pgs/scores/PGS002724/ScoringFiles/PGS002724.txt.gz ftp.ebi.ac.uk (ftp.ebi.ac.uk) をDNSに問いあわせています... 193.62.193.138 ftp.ebi.ac.uk (ftp.ebi.ac.uk)|193.62.193.138|:443 に接続しています... 接続しました。 HTTP による接続要求を送信しました、応答を待っています... 200 OK 長さ: 16818717 (16M) [application/x-gzip] `PGS002724.txt.gz' に保存中 PGS002724.txt.gz 100%[===================>] 16.04M 7.16MB/s in 2.2s 2022-11-17 14:37:16 (7.16 MB/s) - `PGS002724.txt.gz' へ保存完了 [16818717/16818717]
[Stage 59:===============> (24 + 8) / 89]
!wget https://ftp.ebi.ac.uk/pub/databases/spot/pgs/scores/PGS002725/ScoringFiles/PGS002725.txt.gz
--2022-11-17 14:37:19-- https://ftp.ebi.ac.uk/pub/databases/spot/pgs/scores/PGS002725/ScoringFiles/PGS002725.txt.gz ftp.ebi.ac.uk (ftp.ebi.ac.uk) をDNSに問いあわせています... 193.62.193.138 ftp.ebi.ac.uk (ftp.ebi.ac.uk)|193.62.193.138|:443 に接続しています... 接続しました。 HTTP による接続要求を送信しました、応答を待っています... 200 OK 長さ: 81992802 (78M) [application/x-gzip] `PGS002725.txt.gz' に保存中 PGS002725.txt.gz 99%[==================> ] 77.85M 1009KB/s eta 0s
[Stage 59:===============> (24 + 8) / 89]
PGS002725.txt.gz 99%[==================> ] 78.12M 1.02MB/s eta 0s PGS002725.txt.gz 100%[===================>] 78.19M 1.04MB/s in 58s 2022-11-17 14:38:18 (1.34 MB/s) - `PGS002725.txt.gz' へ保存完了 [81992802/81992802]
!gunzip PGS002724.txt.gz
!gunzip PGS002725.txt.gz
!ls
PGS002724.txt PGS002725.txt Tutorial_2022_11_15.ipynb hail-20221111-1605-0.2.105-acd89e80c345.log hail-20221114-1432-0.2.105-acd89e80c345.log outputs prs-models
[Stage 59:===============> (25 + 8) / 89]
!mkdir prs-models
[Stage 59:===============> (24 + 8) / 89]
!mv PGS002724.txt prs-models/
[Stage 59:================> (26 + 8) / 89]
!mv PGS002725.txt prs-models/
!ls prs-models
PGS002724.txt PGS002725.txt
!head prs-models/PGS002724.txt -n 30
###PGS CATALOG SCORING FILE - see https://www.pgscatalog.org/downloads/#dl_ftp_scoring for additional information #format_version=2.0 ##POLYGENIC SCORE (PGS) INFORMATION #pgs_id=PGS002724 #pgs_name=GIGASTROKE_iPGS_EUR #trait_reported=Ischemic stroke #trait_mapped=stroke|Ischemic stroke #trait_efo=EFO_0000712|HP_0002140 #genome_build=GRCh37 #variants_number=1213574 #weight_type=NR ##SOURCE INFORMATION #pgp_id=PGP000333 #citation=Mishra A et al. Nature (2022). doi:10.1038/s41586-022-05165-3 chr_name chr_position effect_allele other_allele effect_weight 1 752721 G A 50.2009138795063 1 754182 G A 141.073654032741 1 760912 T C 180.556536852976 1 768448 A G -74.6438253333578 1 779322 G A -137.02495892717 1 838555 A C -128.635559661359 1 846808 T C -103.154370468722 1 853954 A C -128.44445344713 1 854250 G A -67.2278476132285 1 861808 G A 327.215410265433 1 863124 T G 403.185761038646 1 864938 A G -173.628855846674 1 870645 C T -80.5551733833478 1 873558 T G -261.898526036879 1 879317 T C 14.4910134312575
model_PGS002724 = hl.import_table('prs-models/PGS002724.txt', impute=True, force=True, comment='#')
2022-11-23 12:56:23.980 Hail: INFO: Reading table to impute column types 2022-11-23 12:56:31.052 Hail: INFO: Finished type imputation (0 + 1) / 1] Loading field 'chr_name' as type int32 (imputed) Loading field 'chr_position' as type int32 (imputed) Loading field 'effect_allele' as type str (imputed) Loading field 'other_allele' as type str (imputed) Loading field 'effect_weight' as type float64 (imputed)
model_PGS002725 = hl.import_table('prs-models/PGS002725.txt', impute=True, force=True, comment='#')
2022-11-23 12:57:00.648 Hail: INFO: Reading table to impute column types 2022-11-23 12:57:16.612 Hail: INFO: Finished type imputation (1 + 1) / 2] Loading field 'chr_name' as type int32 (imputed) Loading field 'chr_position' as type int32 (imputed) Loading field 'effect_allele' as type str (imputed) Loading field 'other_allele' as type str (imputed) Loading field 'effect_weight' as type float64 (imputed)
下記のコードは、PRSモデルに含まれるバリアントの個数を表示します。
model_PGS002724.count()
[Stage 23:> (0 + 1) / 1]
1213574
model_PGS002725.count()
[Stage 24:> (0 + 2) / 2]
6010730
1213574
と 6010730
と表示されました。
これは、PRSモデル 002724 と 002725 に含まれるバリアントの個数がそれぞれ 1213574, 6010730 個であることを意味します。
下記のコードは、読み込んだ PRS モデルの最初の 5 行を表示します。
model_PGS002724.show(5)
showing top 5 rows
model_PGS002725.show(5)
showing top 5 rows
下記のコードは、読み込んだ PRS モデルに variantID
を追加します。
model_PGS002724 = model_PGS002724.annotate(
variantID = hl.str(model_PGS002724.chr_name) + ":" + hl.str(model_PGS002724.chr_position)
)
model_PGS002725 = model_PGS002725.annotate(
variantID = hl.str(model_PGS002725.chr_name) + ":" + hl.str(model_PGS002725.chr_position)
)
下記のコードは variantID
を追加した後の最初の 5 行を表示します。
model_PGS002724.show(5)
showing top 5 rows
model_PGS002725.show(5)
showing top 5 rows
variantID
のカラムが追加されていることが分かります。
下記のコードは、PRSモデルのバリアント情報を variantID で検索できるようにします。
model_PGS002724 = model_PGS002724.key_by('variantID')
model_PGS002725 = model_PGS002725.key_by('variantID')
mt1_filt.rows().show(5)
[Stage 33:===========================================> (3 + 1) / 4]
showing top 5 rows
下記のコードは、ゲノムデータと PRS モデルに共通するバリアントを抽出します。
mt_match_PGS002724 = mt1_filt.annotate_rows(**model_PGS002724[mt1_filt.variantID])
mt_match_PGS002724 = mt_match_PGS002724.filter_rows(hl.is_defined(mt_match_PGS002724.effect_weight))
mt_match_PGS002725 = mt1_filt.annotate_rows(**model_PGS002725[mt1_filt.variantID])
mt_match_PGS002725 = mt_match_PGS002725.filter_rows(hl.is_defined(mt_match_PGS002725.effect_weight))
mt_match_PGS002724.rows().show(5)
2022-11-23 13:05:08.949 Hail: INFO: Ordering unsorted dataset with network shuffle 2022-11-23 13:05:16.082 Hail: INFO: Ordering unsorted dataset with network shuffle 2022-11-23 13:11:19.282 Hail: INFO: Ordering unsorted dataset with network shuffle [Stage 47:> (0 + 1) / 1]
showing top 5 rows
model_PGS002724.filter(model_PGS002724.chr_name==1).count()
[Stage 48:> (0 + 1) / 1]
99864
model_PGS002725.filter(model_PGS002725.chr_name==1).count()
[Stage 66:=============================> (1 + 1) / 2]
475080
今後のステップの実行時間を短縮するため(ここを第2のチェックポイントとするため)、一旦抽出したゲノムデータを保存しておきます。
下記のコードは、抽出したゲノムデータを outputs/chrAll.filtered.matched_PGS002724.mt
, outputs/chrAll.filtered.matched_PGS002725.mt
に保存します。
mt_match_PGS002724.write('outputs/chrAll.filtered.matched_PGS002724.mt', overwrite=True)
2022-11-21 16:25:48.517 Hail: INFO: Ordering unsorted dataset with network shuffle 2022-11-21 16:25:53.704 Hail: INFO: Ordering unsorted dataset with network shuffle 2022-11-21 16:34:47.212 Hail: INFO: Ordering unsorted dataset with network shuffle 2022-11-21 16:51:37.284 Hail: INFO: wrote matrix table with 1211967 rows and 2318 columns in 89 partitions to outputs/chrAll.filtered.matched_PGS002724.mt
mt_match_PGS002725.write('outputs/chrAll.filtered.matched_PGS002725.mt', overwrite=True)
2022-11-21 21:56:42.035 Hail: INFO: Ordering unsorted dataset with network shuffle 2022-11-21 21:56:54.561 Hail: INFO: Ordering unsorted dataset with network shuffle 2022-11-21 22:13:20.912 Hail: INFO: Ordering unsorted dataset with network shuffle 2022-11-21 22:52:14.176 Hail: INFO: wrote matrix table with 5930286 rows and 2318 columns in 89 partitions to outputs/chrAll.filtered.matched_PGS002725.mt
len(dict(mt_match_PGS002724.aggregate_rows(hl.agg.counter(mt_match_PGS002724.variantID))))
2022-11-23 13:20:39.449 Hail: INFO: Ordering unsorted dataset with network shuffle 2022-11-23 13:20:51.962 Hail: INFO: Ordering unsorted dataset with network shuffle 2022-11-23 13:26:39.054 Hail: INFO: Ordering unsorted dataset with network shuffle [Stage 57:=============================> (1 + 1) / 2]
1211916
len(dict(mt_match_PGS002725.aggregate_rows(hl.agg.counter(mt_match_PGS002725.variantID))))
2022-11-23 14:05:53.594 Hail: INFO: Ordering unsorted dataset with network shuffle 2022-11-23 14:06:12.039 Hail: INFO: Ordering unsorted dataset with network shuffle 2022-11-23 14:14:15.707 Hail: INFO: Ordering unsorted dataset with network shuffle [Stage 67:=============================> (1 + 1) / 2]
5929754
下記のコードは、ゲノムデータのアリル情報とPRSモデルのアリル情報を比較し、合致しているかをチェックします。
flip_PGS002724 = hl.case().when(
(mt_match_PGS002724.effect_allele == mt_match_PGS002724.alleles[0])
& (mt_match_PGS002724.other_allele == mt_match_PGS002724.alleles[1]), True ).when(
(mt_match_PGS002724.effect_allele == mt_match_PGS002724.alleles[1])
& (mt_match_PGS002724.other_allele == mt_match_PGS002724.alleles[0]), False ).or_missing()
flip_PGS002725 = hl.case().when(
(mt_match_PGS002725.effect_allele == mt_match_PGS002725.alleles[0])
& (mt_match_PGS002725.other_allele == mt_match_PGS002725.alleles[1]), True ).when(
(mt_match_PGS002725.effect_allele == mt_match_PGS002725.alleles[1])
& (mt_match_PGS002725.other_allele == mt_match_PGS002725.alleles[0]), False ).or_missing()
mt_match_PGS002724 = mt_match_PGS002724.annotate_rows(flip=flip_PGS002724)
mt_match_PGS002725 = mt_match_PGS002725.annotate_rows(flip=flip_PGS002725)
下記のコードは、各々のバリアントについて研究対象者の持っているアリル数(mt_match_PGS002724.DS
)とバリアントの重み(mt_match_PGS002724.effect_weight
)を掛け合わせて、ゲノムデータとPRSモデルの共通するバリアントについて足し合わせます。 これにより、PRSを計算することができます。
prs_PGS002724 = hl.agg.sum(hl.float64(mt_match_PGS002724.effect_weight) *
hl.if_else( mt_match_PGS002724.flip,
2 - mt_match_PGS002724.DS.first(),
mt_match_PGS002724.DS.first()))
mt_match_PGS002724 = mt_match_PGS002724.annotate_cols(prs=prs_PGS002724)
PGS002725 についても同様のことを行ってみましょう。
prs_PGS002725 = hl.agg.sum(hl.float64(mt_match_PGS002725.effect_weight) *
hl.if_else( mt_match_PGS002725.flip,
2 - mt_match_PGS002725.DS.first(),
mt_match_PGS002725.DS.first()))
mt_match_PGS002725 = mt_match_PGS002725.annotate_cols(prs=prs_PGS002725)
下記のコードは、PRSの値を表示します。
mt_match_PGS002724.cols().show(5)
2022-11-23 16:24:35.335 Hail: WARN: cols(): Resulting column table is sorted by 'col_key'. To preserve matrix table column order, first unkey columns with 'key_cols_by()' 2022-11-23 16:29:55.700 Hail: INFO: Ordering unsorted dataset with network shuffle 2022-11-23 16:30:04.670 Hail: INFO: Ordering unsorted dataset with network shuffle 2022-11-23 16:36:22.299 Hail: INFO: Ordering unsorted dataset with network shuffle [Stage 76:=======================================================>(88 + 1) / 89]
showing top 5 rows
mt_match_PGS002725.cols().show(5)
2022-11-23 16:53:43.532 Hail: INFO: Ordering unsorted dataset with network shuffle 2022-11-23 16:54:00.362 Hail: INFO: Ordering unsorted dataset with network shuffle 2022-11-23 17:01:50.156 Hail: INFO: Ordering unsorted dataset with network shuffle [Stage 87:=======================================================>(88 + 1) / 89]
showing top 5 rows
下記のコードは、PRSの分布を表示します。
p = hl.plot.histogram(mt_match_PGS002724.prs, title="PRS Histogram", legend="PGS002724", bins=20)
show(p)
2022-11-23 17:18:12.181 Hail: INFO: Ordering unsorted dataset with network shuffle 2022-11-23 17:18:20.554 Hail: INFO: Ordering unsorted dataset with network shuffle 2022-11-23 17:23:52.106 Hail: INFO: Ordering unsorted dataset with network shuffle 2022-11-23 17:36:20.617 Hail: INFO: Ordering unsorted dataset with network shuffle 2022-11-23 17:36:26.549 Hail: INFO: Ordering unsorted dataset with network shuffle 2022-11-23 17:42:27.683 Hail: INFO: Ordering unsorted dataset with network shuffle [Stage 109:======================================================>(88 + 1) / 89]
下記のコードは、PRSの計算結果を chrAll.filtered.PGS002724.PRS.txt ファイルに保存します。
mt_match_PGS002724.cols().export('chrAll.filtered.PGS002724.PRS.txt')
2022-11-23 17:57:54.529 Hail: INFO: Ordering unsorted dataset with network shuffle 2022-11-23 17:58:03.476 Hail: INFO: Ordering unsorted dataset with network shuffle 2022-11-23 18:03:44.519 Hail: INFO: Ordering unsorted dataset with network shuffle 2022-11-23 18:11:38.831 Hail: INFO: Coerced sorted dataset=======>(88 + 1) / 89] 2022-11-23 18:11:40.349 Hail: INFO: merging 17 files totalling 47.5K... 2022-11-23 18:11:40.442 Hail: INFO: while writing: chrAll.filtered.PGS002724.PRS.txt merge time: 91.440ms
PGS002725 についても同様のことを行ってみましょう。
p = hl.plot.histogram(mt_match_PGS002725.prs, title="PRS Histogram", legend="PGS002725", bins=20)
show(p)
2022-11-23 18:34:36.295 Hail: INFO: Ordering unsorted dataset with network shuffle 2022-11-23 18:34:57.684 Hail: INFO: Ordering unsorted dataset with network shuffle 2022-11-23 18:43:45.743 Hail: INFO: Ordering unsorted dataset with network shuffle 2022-11-23 18:59:19.136 Hail: INFO: Ordering unsorted dataset with network shuffle 2022-11-23 18:59:31.370 Hail: INFO: Ordering unsorted dataset with network shuffle 2022-11-23 19:08:03.664 Hail: INFO: Ordering unsorted dataset with network shuffle [Stage 143:======================================================>(88 + 1) / 89]
mt_match_PGS002725.cols().export('chrAll.filtered.PGS002725.PRS.txt')
2022-11-23 19:34:31.787 Hail: INFO: Ordering unsorted dataset with network shuffle 2022-11-23 19:34:49.417 Hail: INFO: Ordering unsorted dataset with network shuffle 2022-11-23 19:43:15.359 Hail: INFO: Ordering unsorted dataset with network shuffle 2022-11-23 19:54:13.539 Hail: INFO: Coerced sorted dataset========(89 + 0) / 89] 2022-11-23 19:54:13.800 Hail: INFO: merging 17 files totalling 47.5K... 2022-11-23 19:54:13.828 Hail: INFO: while writing: chrAll.filtered.PGS002725.PRS.txt merge time: 26.715ms
ここまでで 2 つの脳梗塞PRSモデルを用いて、PRSスコア値を計算しました。
このチュートリアルの最後に、そのPRSスコア値を比較してみます。
下記のコードは、PRSモデル PGS002724 を用いて計算した PRS スコア値を読み込みます。
prs_PGS002724 = hl.import_table('chrAll.filtered.PGS002724.PRS.txt', impute=True, force=True)
2022-11-23 20:19:37.892 Hail: INFO: Reading table to impute column types 2022-11-23 20:19:38.596 Hail: INFO: Finished type imputation Loading field 's' as type str (imputed) Loading field 'prs' as type float64 (imputed)
prs_PGS002725 = hl.import_table('chrAll.filtered.PGS002725.PRS.txt', impute=True, force=True)
2022-11-23 20:20:10.640 Hail: INFO: Reading table to impute column types 2022-11-23 20:20:11.192 Hail: INFO: Finished type imputation Loading field 's' as type str (imputed) Loading field 'prs' as type float64 (imputed)
下記のコードは、PRS スコア値を subject ID (変数名 s
) で検索できるようにします。
prs_PGS002724 = prs_PGS002724.key_by('s')
prs_PGS002725 = prs_PGS002725.key_by('s')
下記のコードは、2 種類の PRS スコア値を subject ID (変数名 s) で突合し、データマージを行います。
prs_merge = prs_PGS002724.rename({'s':'subjectID', 'prs':'PGS002724'})
prs_merge = prs_merge.annotate(PGS002725 = prs_PGS002725[prs_merge.subjectID].prs)
下記のコードは、データマージした結果を表示します。
prs_merge.show(5)
2022-11-23 20:24:34.655 Hail: INFO: Coerced sorted dataset 2022-11-23 20:24:35.058 Hail: INFO: Coerced sorted dataset
showing top 5 rows
下記のコードは、PGS002724
と PGS002725
のスコア値をプロットします。
p = hl.plot.scatter(prs_merge.PGS002724, prs_merge.PGS002725, xlabel="PGS002724", ylabel="PGS002725")
show(p)
2022-11-23 20:25:35.986 Hail: INFO: Coerced sorted dataset 2022-11-23 20:25:36.310 Hail: INFO: Coerced sorted dataset
下記のコードは、PRS スコア値の相関係数を計算します。
to_pandas()
関数を用いることで、hail
特有のデータタイプから、python
でよく使われる pandas
ライブラリのデータフレームへと変換することができます。
そうすることで、hail
の関数だけでなく、python
の様々な機能を使って分析することが可能になります。
prs_merge_pandas = prs_merge.to_pandas()
2022-11-23 20:26:30.449 Hail: INFO: Coerced sorted dataset 2022-11-23 20:26:30.729 Hail: INFO: Coerced sorted dataset
print(prs_merge_pandas.corr())
PGS002724 PGS002725 PGS002724 1.000000 0.168966 PGS002725 0.168966 1.000000