Ccmmutty logo
Commutty IT
7 min read

DNA変異を解析するPythonパッケージを開発し論文を出した

https://cdn.magicode.io/media/notebox/17720a29-5b3e-499e-95b8-4affb853b6c7.jpeg

はじめに

今回は私が1年ほど前にリリースしたViola-SVというPythonパッケージをご紹介します。論文はBioinformaticsという雑誌に掲載されています。
やや今更感は否めませんが、Magicodeリリース祝いということで書くことにしました笑
かなりニッチなソフトなので、紹介したところで皆さんが使うことはないと思います。しかし「記事内でプログラムを編集・実行できる」という面白い機能があるため、実験的にコードを紹介することにしました。ご容赦ください。

Viola-SVとは?

Viola-SVはDNAの構造多型(Structural Variant; SV)を解析するためのPythonパッケージです。 最近は次世代シーケンサという遺伝子配列を端から端まで読めるマシーンがあり、個体が有する配列異常をゲノム横断的に検出できる可能性を秘めています。DNAの変異にはさまざまなカテゴリーがありますが、このパッケージでは50bp以上の大きめな変異およびフィラデルフィア染色体の代表されるような転座、総称して構造多型(SV)を主に扱っています。
SVのデータはVCF(Variant Call Format)ファイルという形式に記録されることが主流です。しかしVCFファイルは人間にとっては可読性が高いのですが、プログラム処理上はかなり厄介なファイル形式でした。
#CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO	FORMAT	mouse01_N	mouse01_T
chr1	100000	M1	N	<DEL>	.	MinSomaticScore	IMPRECISE;SVTYPE=DEL;SVLEN=-100000;END=200000;CIPOS=-51,52;CIEND=-51,52;SOMATIC;SOMATICSCORE=10	PR:SR	21,0:10,0	43,4:15,3
chr1	200000	MD1	N	<DEL>	.	MinSomaticScore	IMPRECISE;SVTYPE=DEL;SVLEN=-100000;END=300000;CIPOS=-13,13;CIEND=-51,52;SOMATIC;SOMATICSCORE=10	PR:SR	21,0:10,0	43,4:15,3
chr1	300000	ML1	N	<DEL>	.	MinSomaticScore	IMPRECISE;SVTYPE=DEL;SVLEN=-100000;END=400000;CIPOS=-51,52;CIEND=-51,52;SOMATIC;SOMATICSCORE=10	PR:SR	21,0:10,0	43,4:15,3
chr1	400000	MG1	N	<DEL>	.	MinSomaticScore	IMPRECISE;SVTYPE=DEL;SVLEN=-100000;END=500000;CIPOS=-51,52;CIEND=-51,52;SOMATIC;SOMATICSCORE=10	PR:SR	21,0:10,0	43,4:15,3
chr2	100000	MDL1	N	<DUP:TANDEM>	.	MinSomaticScore	IMPRECISE;SVTYPE=DUP;SVLEN=100000;END=200000;CIPOS=-30,30;CIEND=-31,31;SOMATIC;SOMATICSCORE=10	PR:SR	21,0:10,0	43,4:15,3
chr2	100500	MDG1	N	<DUP:TANDEM>	.	MinSomaticScore	IMPRECISE;SVTYPE=DUP;SVLEN=100000;END=200500;CIPOS=-51,52;CIEND=-51,52;SOMATIC;SOMATICSCORE=10	PR:SR	21,0:10,0	43,4:15,3
chr2	200000	MLG1	N	<DUP:TANDEM>	.	MinSomaticScore	IMPRECISE;SVTYPE=DUP;SVLEN=100000;END=300000;CIPOS=-51,52;CIEND=-51,52;SOMATIC;SOMATICSCORE=10	PR:SR	21,0:10,0	43,4:15,3
chr3	100000	MDLG1o	N	N]chr4:200000]	.	MinSomaticScore	IMPRECISE;SVTYPE=BND;CIPOS=-100,100;MATEID=MDLG1h;BND_DEPTH=49;MATE_BND_DEPTH=35;SOMATIC;SOMATICSCORE=12	PR	30,1	68,9
chr4	200000	MDLG1h	N	N]chr3:100000]	.	MinSomaticScore	IMPRECISE;SVTYPE=BND;CIPOS=-100,100;MATEID=MDLG1o;BND_DEPTH=49;MATE_BND_DEPTH=35;SOMATIC;SOMATICSCORE=12	PR	30,1	68,9
なぜ厄介なのかを語れば長くなるのですが、一番は「1セルに複数の値が格納されている列がある」という部分です。VCFは上例の通りTABで仕切られたテーブルデータなのですが、行と列を指定しても値が一意に決まらないことがあります。例えば1行目のINFO列を見ると、
IMPRECISE;SVTYPE=DEL;SVLEN=-100000;END=200000;CIPOS=-51,52;CIEND=-51,52;SOMATIC;SOMATICSCORE=10
となっていますね。この中の値に到達するには、さらなる処理が必要です。なんて複雑な構造!
Viola-SVの最初のモチベーションは、VCFファイルの再構成でした。論文化するために色々な機能を載せてはいますが、本質は「VCFの再構成」につきます。

Tidy Dataという考え方

当時所属していた大学のITサークルで「Tidy Data」という考え方に出会いました。 Tidy Dataとは
  • 各列は一つの変数
  • 各行は一つの観測
  • 各セルには一つのデータのみが格納される という条件を満たしたテーブルデータのことです。 こうすることで、行と列の添字を指定するだけで特定の値に到達できるようになります。
これだけでもコードがシンプルになりますね。ということでTidy DataをViola-SVに採用することにしました。VCFファイルをTidy Dataに分解する、という作戦です。
私はそれぞれのTidy Dataに一意制約を課したID列(複合ID可)を必ず設け、いわゆるボイスコッド第三正規形に近いデータ構造を採用しました。
Tidy Dataとかボイスコッド正規形などについてこれ以上細かい説明をするのは面倒なので、あとは調べてみてください。

実際に動かしてみる

さて、Viola-SVを実際に動かしてみます。興味のない方は読み飛ばしてください。
まずはMagicode環境にViola-SVをインストールします。
python
!pip install setuptools==57.4.0 # 通常この行は不要です。Magicode環境ではこれがないとエラーになります。
!pip install Viola-SV

Collecting setuptools==57.4.0
Downloading setuptools-57.4.0-py3-none-any.whl (819 kB) [?25l |▍ | 10 kB 13.1 MB/s eta 0:00:01 |▉ | 20 kB 15.8 MB/s eta 0:00:01 |█▏ | 30 kB 14.5 MB/s eta 0:00:01 |█▋ | 40 kB 6.1 MB/s eta 0:00:01 |██ | 51 kB 4.6 MB/s eta 0:00:01 |██▍ | 61 kB 5.3 MB/s eta 0:00:01 |██▉ | 71 kB 6.1 MB/s eta 0:00:01 |███▏ | 81 kB 6.2 MB/s eta 0:00:01 |███▋ | 92 kB 6.9 MB/s eta 0:00:01 |████ | 102 kB 6.3 MB/s eta 0:00:01 |████▍ | 112 kB 6.3 MB/s eta 0:00:01 |████▉ | 122 kB 6.3 MB/s eta 0:00:01 |█████▏ | 133 kB 6.3 MB/s eta 0:00:01 |█████▋ | 143 kB 6.3 MB/s eta 0:00:01 |██████ | 153 kB 6.3 MB/s eta 0:00:01 |██████▍ | 163 kB 6.3 MB/s eta 0:00:01 |██████▉ | 174 kB 6.3 MB/s eta 0:00:01 |███████▏ | 184 kB 6.3 MB/s eta 0:00:01 |███████▋ | 194 kB 6.3 MB/s eta 0:00:01 |████████ | 204 kB 6.3 MB/s eta 0:00:01 |████████▍ | 215 kB 6.3 MB/s eta 0:00:01
|████████▉ | 225 kB 6.3 MB/s eta 0:00:01 |█████████▏ | 235 kB 6.3 MB/s eta 0:00:01 |█████████▋ | 245 kB 6.3 MB/s eta 0:00:01 |██████████ | 256 kB 6.3 MB/s eta 0:00:01 |██████████▍ | 266 kB 6.3 MB/s eta 0:00:01 |██████████▉ | 276 kB 6.3 MB/s eta 0:00:01 |███████████▏ | 286 kB 6.3 MB/s eta 0:00:01 |███████████▋ | 296 kB 6.3 MB/s eta 0:00:01 |████████████ | 307 kB 6.3 MB/s eta 0:00:01 |████████████▍ | 317 kB 6.3 MB/s eta 0:00:01 |████████████▉ | 327 kB 6.3 MB/s eta 0:00:01 |█████████████▏ | 337 kB 6.3 MB/s eta 0:00:01 |█████████████▋ | 348 kB 6.3 MB/s eta 0:00:01 |██████████████ | 358 kB 6.3 MB/s eta 0:00:01 |██████████████▍ | 368 kB 6.3 MB/s eta 0:00:01 |██████████████▉ | 378 kB 6.3 MB/s eta 0:00:01 |███████████████▏ | 389 kB 6.3 MB/s eta 0:00:01 |███████████████▋ | 399 kB 6.3 MB/s eta 0:00:01 |████████████████ | 409 kB 6.3 MB/s eta 0:00:01 |████████████████▍ | 419 kB 6.3 MB/s eta 0:00:01 |████████████████▉ | 430 kB 6.3 MB/s eta 0:00:01 |█████████████████▏ | 440 kB 6.3 MB/s eta 0:00:01 |█████████████████▋ | 450 kB 6.3 MB/s eta 0:00:01 |██████████████████ | 460 kB 6.3 MB/s eta 0:00:01 |██████████████████▍ | 471 kB 6.3 MB/s eta 0:00:01 |██████████████████▉ | 481 kB 6.3 MB/s eta 0:00:01 |███████████████████▏ | 491 kB 6.3 MB/s eta 0:00:01 |███████████████████▋ | 501 kB 6.3 MB/s eta 0:00:01 |████████████████████ | 512 kB 6.3 MB/s eta 0:00:01 |████████████████████▍ | 522 kB 6.3 MB/s eta 0:00:01 |████████████████████▉ | 532 kB 6.3 MB/s eta 0:00:01 |█████████████████████▏ | 542 kB 6.3 MB/s eta 0:00:01 |█████████████████████▋ | 552 kB 6.3 MB/s eta 0:00:01 |██████████████████████ | 563 kB 6.3 MB/s eta 0:00:01 |██████████████████████▍ | 573 kB 6.3 MB/s eta 0:00:01 |██████████████████████▉ | 583 kB 6.3 MB/s eta 0:00:01 |███████████████████████▏ | 593 kB 6.3 MB/s eta 0:00:01 |███████████████████████▋ | 604 kB 6.3 MB/s eta 0:00:01 |████████████████████████ | 614 kB 6.3 MB/s eta 0:00:01 |████████████████████████▍ | 624 kB 6.3 MB/s eta 0:00:01 |████████████████████████▉ | 634 kB 6.3 MB/s eta 0:00:01 |█████████████████████████▏ | 645 kB 6.3 MB/s eta 0:00:01 |█████████████████████████▋ | 655 kB 6.3 MB/s eta 0:00:01 |██████████████████████████ | 665 kB 6.3 MB/s eta 0:00:01 |██████████████████████████▍ | 675 kB 6.3 MB/s eta 0:00:01 |██████████████████████████▉ | 686 kB 6.3 MB/s eta 0:00:01 |███████████████████████████▏ | 696 kB 6.3 MB/s eta 0:00:01 |███████████████████████████▋ | 706 kB 6.3 MB/s eta 0:00:01 |████████████████████████████ | 716 kB 6.3 MB/s eta 0:00:01 |████████████████████████████▍ | 727 kB 6.3 MB/s eta 0:00:01 |████████████████████████████▉ | 737 kB 6.3 MB/s eta 0:00:01 |█████████████████████████████▏ | 747 kB 6.3 MB/s eta 0:00:01 |█████████████████████████████▋ | 757 kB 6.3 MB/s eta 0:00:01 |██████████████████████████████ | 768 kB 6.3 MB/s eta 0:00:01 |██████████████████████████████▍ | 778 kB 6.3 MB/s eta 0:00:01 |██████████████████████████████▉ | 788 kB 6.3 MB/s eta 0:00:01 |███████████████████████████████▏| 798 kB 6.3 MB/s eta 0:00:01 |███████████████████████████████▋| 808 kB 6.3 MB/s eta 0:00:01 |████████████████████████████████| 819 kB 6.3 MB/s [?25h
Installing collected packages: setuptools Attempting uninstall: setuptools Found existing installation: setuptools 60.5.0
Uninstalling setuptools-60.5.0:
Successfully uninstalled setuptools-60.5.0
Successfully installed setuptools-57.4.0
Collecting Viola-SV
Downloading Viola_SV-1.0.0.dev10-py3-none-any.whl (61 kB) [?25l |█████▎ | 10 kB 12.6 MB/s eta 0:00:01 |██████████▋ | 20 kB 15.2 MB/s eta 0:00:01 |████████████████ | 30 kB 9.4 MB/s eta 0:00:01 |█████████████████████▎ | 40 kB 5.9 MB/s eta 0:00:01 |██████████████████████████▌ | 51 kB 6.2 MB/s eta 0:00:01 |███████████████████████████████▉| 61 kB 6.9 MB/s eta 0:00:01 |████████████████████████████████| 61 kB 634 kB/s [?25h
Collecting intervaltree
Downloading intervaltree-3.1.0.tar.gz (32 kB) Preparing metadata (setup.py) ... [?25l
-
 done [?25hRequirement already satisfied: pandas in /srv/conda/envs/notebook/lib/python3.7/site-packages (from Viola-SV) (1.1.5)
Collecting biopython Downloading biopython-1.79-cp37-cp37m-manylinux_2_5_x86_64.manylinux1_x86_64.whl (2.3 MB) [?25l |▏ | 10 kB 23.6 MB/s eta 0:00:01 |▎ | 20 kB 30.6 MB/s eta 0:00:01 |▍ | 30 kB 39.2 MB/s eta 0:00:01 |▋ | 40 kB 45.4 MB/s eta 0:00:01 |▊ | 51 kB 43.6 MB/s eta 0:00:01 |▉ | 61 kB 48.1 MB/s eta 0:00:01 |█ | 71 kB 34.3 MB/s eta 0:00:01 |█▏ | 81 kB 33.2 MB/s eta 0:00:01 |█▎ | 92 kB 36.0 MB/s eta 0:00:01 |█▍ | 102 kB 13.1 MB/s eta 0:00:01 |█▋ | 112 kB 13.1 MB/s eta 0:00:01 |█▊ | 122 kB 13.1 MB/s eta 0:00:01 |█▉ | 133 kB 13.1 MB/s eta 0:00:01 |██ | 143 kB 13.1 MB/s eta 0:00:01 |██▏ | 153 kB 13.1 MB/s eta 0:00:01 |██▎ | 163 kB 13.1 MB/s eta 0:00:01 |██▍ | 174 kB 13.1 MB/s eta 0:00:01 |██▋ | 184 kB 13.1 MB/s eta 0:00:01 |██▊ | 194 kB 13.1 MB/s eta 0:00:01 |██▉ | 204 kB 13.1 MB/s eta 0:00:01 |███ | 215 kB 13.1 MB/s eta 0:00:01 |███▏ | 225 kB 13.1 MB/s eta 0:00:01 |███▎ | 235 kB 13.1 MB/s eta 0:00:01 |███▌ | 245 kB 13.1 MB/s eta 0:00:01 |███▋ | 256 kB 13.1 MB/s eta 0:00:01 |███▊ | 266 kB 13.1 MB/s eta 0:00:01 |███▉ | 276 kB 13.1 MB/s eta 0:00:01 |████ | 286 kB 13.1 MB/s eta 0:00:01 |████▏ | 296 kB 13.1 MB/s eta 0:00:01 |████▎ | 307 kB 13.1 MB/s eta 0:00:01 |████▌ | 317 kB 13.1 MB/s eta 0:00:01 |████▋ | 327 kB 13.1 MB/s eta 0:00:01 |████▊ | 337 kB 13.1 MB/s eta 0:00:01 |████▉ | 348 kB 13.1 MB/s eta 0:00:01 |█████ | 358 kB 13.1 MB/s eta 0:00:01 |█████▏ | 368 kB 13.1 MB/s eta 0:00:01 |█████▎ | 378 kB 13.1 MB/s eta 0:00:01 |█████▌ | 389 kB 13.1 MB/s eta 0:00:01 |█████▋ | 399 kB 13.1 MB/s eta 0:00:01 |█████▊ | 409 kB 13.1 MB/s eta 0:00:01 |██████ | 419 kB 13.1 MB/s eta 0:00:01 |██████ | 430 kB 13.1 MB/s eta 0:00:01 |██████▏ | 440 kB 13.1 MB/s eta 0:00:01 |██████▎ | 450 kB 13.1 MB/s eta 0:00:01 |██████▌ | 460 kB 13.1 MB/s eta 0:00:01 |██████▋ | 471 kB 13.1 MB/s eta 0:00:01 |██████▊ | 481 kB 13.1 MB/s eta 0:00:01 |███████ | 491 kB 13.1 MB/s eta 0:00:01 |███████ | 501 kB 13.1 MB/s eta 0:00:01 |███████▏ | 512 kB 13.1 MB/s eta 0:00:01 |███████▎ | 522 kB 13.1 MB/s eta 0:00:01 |███████▌ | 532 kB 13.1 MB/s eta 0:00:01 |███████▋ | 542 kB 13.1 MB/s eta 0:00:01 |███████▊ | 552 kB 13.1 MB/s eta 0:00:01 |████████ | 563 kB 13.1 MB/s eta 0:00:01 |████████ | 573 kB 13.1 MB/s eta 0:00:01 |████████▏ | 583 kB 13.1 MB/s eta 0:00:01 |████████▍ | 593 kB 13.1 MB/s eta 0:00:01 |████████▌ | 604 kB 13.1 MB/s eta 0:00:01 |████████▋ | 614 kB 13.1 MB/s eta 0:00:01 |████████▊ | 624 kB 13.1 MB/s eta 0:00:01 |█████████ | 634 kB 13.1 MB/s eta 0:00:01 |█████████ | 645 kB 13.1 MB/s eta 0:00:01 |█████████▏ | 655 kB 13.1 MB/s eta 0:00:01 |█████████▍ | 665 kB 13.1 MB/s eta 0:00:01 |█████████▌ | 675 kB 13.1 MB/s eta 0:00:01 |█████████▋ | 686 kB 13.1 MB/s eta 0:00:01 |█████████▊ | 696 kB 13.1 MB/s eta 0:00:01 |██████████ | 706 kB 13.1 MB/s eta 0:00:01 |██████████ | 716 kB 13.1 MB/s eta 0:00:01 |██████████▏ | 727 kB 13.1 MB/s eta 0:00:01 |██████████▍ | 737 kB 13.1 MB/s eta 0:00:01 |██████████▌ | 747 kB 13.1 MB/s eta 0:00:01
|██████████▋ | 757 kB 13.1 MB/s eta 0:00:01 |██████████▉ | 768 kB 13.1 MB/s eta 0:00:01 |███████████ | 778 kB 13.1 MB/s eta 0:00:01 |███████████ | 788 kB 13.1 MB/s eta 0:00:01 |███████████▏ | 798 kB 13.1 MB/s eta 0:00:01 |███████████▍ | 808 kB 13.1 MB/s eta 0:00:01 |███████████▌ | 819 kB 13.1 MB/s eta 0:00:01 |███████████▋ | 829 kB 13.1 MB/s eta 0:00:01 |███████████▉ | 839 kB 13.1 MB/s eta 0:00:01 |████████████ | 849 kB 13.1 MB/s eta 0:00:01 |████████████ | 860 kB 13.1 MB/s eta 0:00:01 |████████████▏ | 870 kB 13.1 MB/s eta 0:00:01 |████████████▍ | 880 kB 13.1 MB/s eta 0:00:01 |████████████▌ | 890 kB 13.1 MB/s eta 0:00:01 |████████████▋ | 901 kB 13.1 MB/s eta 0:00:01 |████████████▉ | 911 kB 13.1 MB/s eta 0:00:01 |█████████████ | 921 kB 13.1 MB/s eta 0:00:01 |█████████████ | 931 kB 13.1 MB/s eta 0:00:01 |█████████████▎ | 942 kB 13.1 MB/s eta 0:00:01 |█████████████▍ | 952 kB 13.1 MB/s eta 0:00:01 |█████████████▌ | 962 kB 13.1 MB/s eta 0:00:01 |█████████████▋ | 972 kB 13.1 MB/s eta 0:00:01 |█████████████▉ | 983 kB 13.1 MB/s eta 0:00:01 |██████████████ | 993 kB 13.1 MB/s eta 0:00:01 |██████████████ | 1.0 MB 13.1 MB/s eta 0:00:01 |██████████████▎ | 1.0 MB 13.1 MB/s eta 0:00:01 |██████████████▍ | 1.0 MB 13.1 MB/s eta 0:00:01 |██████████████▌ | 1.0 MB 13.1 MB/s eta 0:00:01 |██████████████▋ | 1.0 MB 13.1 MB/s eta 0:00:01 |██████████████▉ | 1.1 MB 13.1 MB/s eta 0:00:01 |███████████████ | 1.1 MB 13.1 MB/s eta 0:00:01 |███████████████ | 1.1 MB 13.1 MB/s eta 0:00:01 |███████████████▎ | 1.1 MB 13.1 MB/s eta 0:00:01 |███████████████▍ | 1.1 MB 13.1 MB/s eta 0:00:01 |███████████████▌ | 1.1 MB 13.1 MB/s eta 0:00:01 |███████████████▊ | 1.1 MB 13.1 MB/s eta 0:00:01 |███████████████▉ | 1.1 MB 13.1 MB/s eta 0:00:01 |████████████████ | 1.1 MB 13.1 MB/s eta 0:00:01 |████████████████ | 1.1 MB 13.1 MB/s eta 0:00:01 |████████████████▎ | 1.2 MB 13.1 MB/s eta 0:00:01 |████████████████▍ | 1.2 MB 13.1 MB/s eta 0:00:01 |████████████████▌ | 1.2 MB 13.1 MB/s eta 0:00:01 |████████████████▊ | 1.2 MB 13.1 MB/s eta 0:00:01 |████████████████▉ | 1.2 MB 13.1 MB/s eta 0:00:01 |█████████████████ | 1.2 MB 13.1 MB/s eta 0:00:01 |█████████████████ | 1.2 MB 13.1 MB/s eta 0:00:01 |█████████████████▎ | 1.2 MB 13.1 MB/s eta 0:00:01 |█████████████████▍ | 1.2 MB 13.1 MB/s eta 0:00:01 |█████████████████▌ | 1.2 MB 13.1 MB/s eta 0:00:01 |█████████████████▊ | 1.3 MB 13.1 MB/s eta 0:00:01 |█████████████████▉ | 1.3 MB 13.1 MB/s eta 0:00:01 |██████████████████ | 1.3 MB 13.1 MB/s eta 0:00:01 |██████████████████▏ | 1.3 MB 13.1 MB/s eta 0:00:01 |██████████████████▎ | 1.3 MB 13.1 MB/s eta 0:00:01 |██████████████████▍ | 1.3 MB 13.1 MB/s eta 0:00:01 |██████████████████▌ | 1.3 MB 13.1 MB/s eta 0:00:01 |██████████████████▊ | 1.3 MB 13.1 MB/s eta 0:00:01 |██████████████████▉ | 1.3 MB 13.1 MB/s eta 0:00:01 |███████████████████ | 1.4 MB 13.1 MB/s eta 0:00:01 |███████████████████▏ | 1.4 MB 13.1 MB/s eta 0:00:01 |███████████████████▎ | 1.4 MB 13.1 MB/s eta 0:00:01 |███████████████████▍ | 1.4 MB 13.1 MB/s eta 0:00:01 |███████████████████▌ | 1.4 MB 13.1 MB/s eta 0:00:01 |███████████████████▊ | 1.4 MB 13.1 MB/s eta 0:00:01 |███████████████████▉ | 1.4 MB 13.1 MB/s eta 0:00:01 |████████████████████ | 1.4 MB 13.1 MB/s eta 0:00:01 |████████████████████▏ | 1.4 MB 13.1 MB/s eta 0:00:01 |████████████████████▎ | 1.4 MB 13.1 MB/s eta 0:00:01 |████████████████████▍ | 1.5 MB 13.1 MB/s eta 0:00:01 |████████████████████▋ | 1.5 MB 13.1 MB/s eta 0:00:01 |████████████████████▊ | 1.5 MB 13.1 MB/s eta 0:00:01 |████████████████████▉ | 1.5 MB 13.1 MB/s eta 0:00:01 |█████████████████████ | 1.5 MB 13.1 MB/s eta 0:00:01 |█████████████████████▏ | 1.5 MB 13.1 MB/s eta 0:00:01 |█████████████████████▎ | 1.5 MB 13.1 MB/s eta 0:00:01 |█████████████████████▍ | 1.5 MB 13.1 MB/s eta 0:00:01 |█████████████████████▋ | 1.5 MB 13.1 MB/s eta 0:00:01 |█████████████████████▊ | 1.5 MB 13.1 MB/s eta 0:00:01 |█████████████████████▉ | 1.6 MB 13.1 MB/s eta 0:00:01 |██████████████████████ | 1.6 MB 13.1 MB/s eta 0:00:01 |██████████████████████▏ | 1.6 MB 13.1 MB/s eta 0:00:01 |██████████████████████▎ | 1.6 MB 13.1 MB/s eta 0:00:01 |██████████████████████▍ | 1.6 MB 13.1 MB/s eta 0:00:01 |██████████████████████▋ | 1.6 MB 13.1 MB/s eta 0:00:01 |██████████████████████▊ | 1.6 MB 13.1 MB/s eta 0:00:01 |██████████████████████▉ | 1.6 MB 13.1 MB/s eta 0:00:01 |███████████████████████ | 1.6 MB 13.1 MB/s eta 0:00:01 |███████████████████████▏ | 1.6 MB 13.1 MB/s eta 0:00:01 |███████████████████████▎ | 1.7 MB 13.1 MB/s eta 0:00:01 |███████████████████████▍ | 1.7 MB 13.1 MB/s eta 0:00:01 |███████████████████████▋ | 1.7 MB 13.1 MB/s eta 0:00:01 |███████████████████████▊ | 1.7 MB 13.1 MB/s eta 0:00:01 |███████████████████████▉ | 1.7 MB 13.1 MB/s eta 0:00:01 |████████████████████████ | 1.7 MB 13.1 MB/s eta 0:00:01 |████████████████████████▏ | 1.7 MB 13.1 MB/s eta 0:00:01 |████████████████████████▎ | 1.7 MB 13.1 MB/s eta 0:00:01 |████████████████████████▍ | 1.7 MB 13.1 MB/s eta 0:00:01 |████████████████████████▋ | 1.8 MB 13.1 MB/s eta 0:00:01 |████████████████████████▊ | 1.8 MB 13.1 MB/s eta 0:00:01 |████████████████████████▉ | 1.8 MB 13.1 MB/s eta 0:00:01 |█████████████████████████ | 1.8 MB 13.1 MB/s eta 0:00:01 |█████████████████████████▏ | 1.8 MB 13.1 MB/s eta 0:00:01 |█████████████████████████▎ | 1.8 MB 13.1 MB/s eta 0:00:01 |█████████████████████████▌ | 1.8 MB 13.1 MB/s eta 0:00:01 |█████████████████████████▋ | 1.8 MB 13.1 MB/s eta 0:00:01 |█████████████████████████▊ | 1.8 MB 13.1 MB/s eta 0:00:01 |█████████████████████████▉ | 1.8 MB 13.1 MB/s eta 0:00:01 |██████████████████████████ | 1.9 MB 13.1 MB/s eta 0:00:01 |██████████████████████████▏ | 1.9 MB 13.1 MB/s eta 0:00:01 |██████████████████████████▎ | 1.9 MB 13.1 MB/s eta 0:00:01 |██████████████████████████▌ | 1.9 MB 13.1 MB/s eta 0:00:01 |██████████████████████████▋ | 1.9 MB 13.1 MB/s eta 0:00:01 |██████████████████████████▊ | 1.9 MB 13.1 MB/s eta 0:00:01 |██████████████████████████▉ | 1.9 MB 13.1 MB/s eta 0:00:01 |███████████████████████████ | 1.9 MB 13.1 MB/s eta 0:00:01 |███████████████████████████▏ | 1.9 MB 13.1 MB/s eta 0:00:01 |███████████████████████████▎ | 1.9 MB 13.1 MB/s eta 0:00:01 |███████████████████████████▌ | 2.0 MB 13.1 MB/s eta 0:00:01 |███████████████████████████▋ | 2.0 MB 13.1 MB/s eta 0:00:01 |███ too many strings
インストールできたので、インポートしてみます。
python
import viola
viola.__version__

'1.0.0.dev10'
さて、今回はサンプルのVCFファイルを読み込んでみましょう。
viola.read_vcf()関数によりVCFファイルをviola.Vcfオブジェクトとして読み込むことができます。なお、viola.read_vcf()はファイルをurlで渡しても動作するように設計してあります。
python
vcf = viola.read_vcf('https://raw.githubusercontent.com/dermasugita/Viola-SV/master/docs/tutorial.manta.vcf')
print(vcf) # vcfはviola.Vcf classのオブジェクトである。

INFO=imprecise,svtype,svlen,end,cipos,ciend,cigar,mateid,event,homlen,homseq,svinslen,svinsseq,left_svinsseq,right_svinsseq,contig,bnd_depth,mate_bnd_depth,somatic,somaticscore,junction_somaticscore,inv3,inv5 Documentation of Vcf object ==> https://dermasugita.github.io/ViolaDocs/docs/html/reference/vcf.html id be1 be2 strand qual svtype 0 test1 chr1:82550461 chr1:82554226 +- None DEL 1 test2 chr1:22814217 chr1:92581132 -- None INV 2 test3 chr1:60567906 chr1:60675941 +- None DEL 3 test4 chr1:69583190 chr1:69590948 +- None DEL 4 test5 chr11:104534877 chr11:104536574 +- None DEL 5 test6_1 chr11:111134697 chr17:26470495 +- None BND 6 test6_2 chr17:26470495 chr11:111134697 -+ None BND
viola.Vcfオブジェクトはシンプルなtidy tablesの集合体となっています。viola.Vcf.table_listでどんなテーブルが含まれているのかを確認できます。
python
print(vcf.table_list)

['positions', 'filters', 'imprecise', 'svtype', 'svlen', 'end', 'cipos', 'ciend', 'cigar', 'mateid', 'event', 'homlen', 'homseq', 'svinslen', 'svinsseq', 'left_svinsseq', 'right_svinsseq', 'contig', 'bnd_depth', 'mate_bnd_depth', 'somatic', 'somaticscore', 'junction_somaticscore', 'inv3', 'inv5', 'formats', 'contigs_meta', 'alts_meta', 'infos_meta', 'formats_meta', 'filters_meta', 'samples_meta']
viola.Vcf.get_table()メソッドでテーブル自体を取得することが可能です。
python
vcf.get_table('positions') # pandas.Dataframeが返る

id chrom1 pos1 chrom2 pos2 ... strand2 ref alt qual svtype
0 test1 chr1 82550461 chr1 82554226 ... - G <DEL> None DEL
1 test2 chr1 22814217 chr1 92581132 ... - T <INV> None INV
2 test3 chr1 60567906 chr1 60675941 ... - T <DEL> None DEL
3 test4 chr1 69583190 chr1 69590948 ... - T <DEL> None DEL
4 test5 chr11 104534877 chr11 104536574 ... - C <DEL> None DEL
5 test6_1 chr11 111134697 chr17 26470495 ... - T T[chr17:26470495[ None BND
6 test6_2 chr17 26470495 chr11 111134697 ... + T ]chr11:111134697]T None BND

7 rows × 11 columns

python
vcf.get_table('svlen')

id value_idx svlen
0 test1 0 -3764
1 test2 0 69766915
2 test3 0 -108034
3 test4 0 -7757
4 test5 0 -1696
それぞれのテーブルはid列で紐づいており、これによりフィルタリングしやすい内部構造を実現しました。 例えばsvlen < 0の行だけを取り出したい場合は、以下のコードで簡単にフィルタリングできます。
python
vcf_filtered = vcf.filter('svlen < 0')
python
vcf_filtered.get_table('positions')

id chrom1 pos1 chrom2 pos2 ... strand2 ref alt qual svtype
0 test1 chr1 82550461 chr1 82554226 ... - G <DEL> None DEL
1 test3 chr1 60567906 chr1 60675941 ... - T <DEL> None DEL
2 test4 chr1 69583190 chr1 69590948 ... - T <DEL> None DEL
3 test5 chr11 104534877 chr11 104536574 ... - C <DEL> None DEL

4 rows × 11 columns

python
vcf_filtered.get_table('svlen')

id value_idx svlen
0 test1 0 -3764
1 test3 0 -108034
2 test4 0 -7757
3 test5 0 -1696
このように、かなりシンプルな文法で情報のフィルタリングが可能になりました。ちなみに元のVCFファイルでsvlen < 0の行だけ取り出すのは大変です。以下に再掲する行をみていただければ言うまでもないと思います。
chr1	100000	M1	N	<DEL>	.	MinSomaticScore	IMPRECISE;SVTYPE=DEL;SVLEN=-100000;END=200000;CIPOS=-51,52;CIEND=-51,52;SOMATIC;SOMATICSCORE=10	PR:SR	21,0:10,0	43,4:15,3

Pythonパッケージを作るのは簡単!

この記事を読んで、「自分もPythonパッケージを作ってみたい」と感じた方もいるかもしれません。
それならば善は急げ、躊躇せずに今すぐ作り始めてみるのが良いと思います。
実際のところPythonパッケージの開発はそこまでハードルが高くありません。筆者自身も別に情報学科出身ではなく、全く関係ない学部出身です。仕事も情報系ではありません。そのためパッケージを書き始めた頃は、プログラミングなんて素人同然でした。
幸いにもネット上にはPythonパッケージ開発に有用なページがたくさんあるので、ググっていくうちにやり方が分かっていくものです。Pythonの公式ドキュメントも有用でした。
筆者の場合、まずはhello worldを出力するだけの小さなパッケージを実験的に作るところから始めました。その先は、多くのネット記事を参考にしたり、pandasやmatplotlibなど他のpackageの構造を解析したりしながら、自分のパッケージを組み立てていきました。
プログラミングはいくらでも試行錯誤ができるのが良いところです。とにかく始めてみれば、なんとかなるものです。
なので皆さんもぜひパッケージ作りに挑戦してみてください。 また、作ったら教えてください。拡散します笑
Pythonパッケージの作り方に関しては、また別の記事で書いてみようと思います。
それでは!

Discussion

コメントにはログインが必要です。