今回は私が1年ほど前にリリースしたViola-SVというPythonパッケージをご紹介します。論文はBioinformaticsという雑誌に掲載されています。
やや今更感は否めませんが、Magicodeリリース祝いということで書くことにしました笑
かなりニッチなソフトなので、紹介したところで皆さんが使うことはないと思います。しかし「記事内でプログラムを編集・実行できる」という面白い機能があるため、実験的にコードを紹介することにしました。ご容赦ください。
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の再構成」につきます。
当時所属していた大学のITサークルで「Tidy Data」という考え方に出会いました。 Tidy Dataとは
という条件を満たしたテーブルデータのことです。 こうすることで、行と列の添字を指定するだけで特定の値に到達できるようになります。
これだけでもコードがシンプルになりますね。ということでTidy DataをViola-SVに採用することにしました。VCFファイルをTidy Dataに分解する、という作戦です。
私はそれぞれのTidy Dataに一意制約を課したID列(複合ID可)を必ず設け、いわゆるボイスコッド第三正規形に近いデータ構造を採用しました。
Tidy Dataとかボイスコッド正規形などについてこれ以上細かい説明をするのは面倒なので、あとは調べてみてください。
さて、Viola-SVを実際に動かしてみます。興味のない方は読み飛ばしてください。
まずはMagicode環境にViola-SVをインストールします。
インストールできたので、インポートしてみます。
さて、今回はサンプルのVCFファイルを読み込んでみましょう。
viola.read_vcf()
関数によりVCFファイルをviola.Vcf
オブジェクトとして読み込むことができます。なお、viola.read_vcf()
はファイルをurlで渡しても動作するように設計してあります。
viola.Vcf
オブジェクトはシンプルなtidy tablesの集合体となっています。viola.Vcf.table_list
でどんなテーブルが含まれているのかを確認できます。
viola.Vcf.get_table()
メソッドでテーブル自体を取得することが可能です。
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
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
の行だけを取り出したい場合は、以下のコードで簡単にフィルタリングできます。
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
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の公式ドキュメントも有用でした。
筆者の場合、まずはhello world
を出力するだけの小さなパッケージを実験的に作るところから始めました。その先は、多くのネット記事を参考にしたり、pandasやmatplotlibなど他のpackageの構造を解析したりしながら、自分のパッケージを組み立てていきました。
プログラミングはいくらでも試行錯誤ができるのが良いところです。とにかく始めてみれば、なんとかなるものです。
なので皆さんもぜひパッケージ作りに挑戦してみてください。 また、作ったら教えてください。拡散します笑
Pythonパッケージの作り方に関しては、また別の記事で書いてみようと思います。
それでは!