macでインフォマティクス

macでインフォマティクス

HTS (NGS) 関連のインフォマティクス情報についてまとめています。

minimap2インデックスに既知バリアント情報を組み込むことで、WGSでのSNVコールを改善する minimap2_index_modifier

 

 リファレンスゲノム配列に対するリードのアライメントは、次世代シーケンサー(NGS)技術によって得られたヒト全ゲノムシーケンスデータの解析における重要なステップの1つである。遺伝的変異の臨床的解釈の結果やゲノムワイド関連研究GWASの結果など、その後の解析ステップの質は、アライメントの結果としてリードの位置が正しく特定されるかどうかに依存している。ヒトNGS全ゲノムシーケンスデータの量は常に増加している。世界中で多くのヒトゲノム配列決定プロジェクトが行われており、その結果、配列決定されたヒトゲノムの遺伝子変異に関する大規模なデータベースが作成されている。このような既知の遺伝子変異に関する情報は、例えばゲノムグラフを作成することで、新しい個体について得られたシーケンスデータを解析する際のリードアライメント段階におけるアライメントの質を向上させるために利用することができる。リードを線形リファレンスゲノムにアライメントする既存の方法はアライメント速度が速いが、リードをゲノムグラフにアライメントする方法はゲノムの可変領域において精度が高い。線形リファレンス配列のインデックスに既知の遺伝子変異を考慮したリードアライメント法を開発することで、両手法の利点を組み合わせることができる。

 本論文では、minimap2_index_modifierツールを紹介する。このツールは、既知の一塩基バリアントとヒト集団特有の挿入/欠失(indel)を用いて、リファレンスゲノムの修正インデックスを構築することができる。modified minimap2インデックスの使用により、バイオインフォマティクスパイプラインを変更することなく、また計算オーバーヘッドを大幅に追加することなく、バリアントコールの品質が向上する。PrecisionFDA Truth Challenge V2のベンチマークデータ(GRCh38リニアリファレンス(GCA_000001405.15)にアライメントされたHG002ショートリードデータ、パラメータk = 27、w = 14)を用いて、Human Pangenome Reference Consortiumの遺伝的バリアントを用いてインデックスを修正すると、偽陰性の遺伝的バリアントの数が9500以上減少し、偽陽性の数が7000以上減少することが実証された。


~minimap2ツールの主な特徴の1つは、minimizersを使用することである。minimizersとは、与えられた配列内に位置する短い部分文字列(または、そうでなければk-measure)のことである。minimizers(最小化子)という用語は、アライメントの過程で比較する必要のある部分文字列の数を最小化するという考え方に由来する。minimizerは、与えられた配列ウィンドウ内の各k-merに与えられたハッシュ関数を適用し、ユニークな整数値を作成することで計算される。次に整数値をソートし、最小値を持つk-measureを最小化子として選択する(論文図2)。ミニマイザーを使用することで、アライメント処理が高速になる。

線形リファレンス配列にインデックスを付ける場合、ハッシュテーブルが作成され、キーは最小化子配列にハッシュ関数を適用した結果であり、値はゲノム参照配列の位置のリストである。

minimap2ツールのインデックスを変更する可能性は、ハッシュテーブルが線形リファレンス配列の所与の位置における最小化因子の数に制約を課さないという事実によって提供され、遺伝的バリアントに関する情報を追加しても、その後のアラインメントアルゴリズムには影響しない。言い換えれば、線形リファレンス配列インデックスは、ゲノムグラフと同様に、遺伝的バリアントの追加によって誘導される分岐の追加を可能にする(論文図3)。(以下省略)

 

インストール

依存

To compile from source, use this version of tools:

  • GCC/G++ 11.4.0+
  • HTSlib v1.17

Github

#HTSlibが必要(samtools解説)先にminimap2_index_modifierをcloneして、minimap2_index_modifierのrootで実行する
git clone --recurse-submodules https://github.com/samtools/htslib.git
cd htslib/
autoheader
autoreconf -i
./configure --prefix=/usr/local
make -j10
make install
cd ../

#本体(minimap2のfolk)*1(注;特別なoptionを使わない限り本家minimap2と同じものだが、パスを通す時は本家と混同しないように注意)
git clone https://github.com/ispras/minimap2_index_modifier.git
cd minimap2_index_modifier/
make -j10

#docker image(hub)
docker pull egorguga/minimap2_index_modifier:2.24

> ./minimap2

Usage: minimap2 [options] <target.fa>|<target.idx> [query.fa] [...]

Options:

  Indexing:

    -H           use homopolymer-compressed k-mer (preferrable for PacBio)

    -k INT       k-mer size (no larger than 28) [15]

    -w INT       minimizer window size [10]

    -I NUM       split index for every ~NUM input bases [4G]

    -d FILE      dump index to FILE

    --vcf-file-with-variants FILE      pass VCF FILE to modify index

    --parse-haplotype parse haplotypes from VCF to generate real SNP combinations, otherwise use all.

  Mapping:

    -f FLOAT     filter out top FLOAT fraction of repetitive minimizers [0.0002]

    -g NUM       stop chain enlongation if there are no minimizers in INT-bp [5000]

    -G NUM       max intron length (effective with -xsplice; changing -r) [200k]

    -F NUM       max fragment length (effective with -xsr or in the fragment mode) [800]

    -r NUM[,NUM] chaining/alignment bandwidth and long-join bandwidth [500,20000]

    -n INT       minimal number of minimizers on a chain [3]

    -m INT       minimal chaining score (matching bases minus log gap penalty) [40]

    -X           skip self and dual mappings (for the all-vs-all mode)

    -p FLOAT     min secondary-to-primary score ratio [0.8]

    -N INT       retain at most INT secondary alignments [5]

  Alignment:

    -A INT       matching score [2]

    -B INT       mismatch penalty (larger value for lower divergence) [4]

    -O INT[,INT] gap open penalty [4,24]

    -E INT[,INT] gap extension penalty; a k-long gap costs min{O1+k*E1,O2+k*E2} [2,1]

    -z INT[,INT] Z-drop score and inversion Z-drop score [400,200]

    -s INT       minimal peak DP alignment score [80]

    -u CHAR      how to find GT-AG. f:transcript strand, b:both strands, n:don't match GT-AG [n]

  Input/Output:

    -a           output in the SAM format (PAF by default)

    -o FILE      output alignments to FILE [stdout]

    -L           write CIGAR with >65535 ops at the CG tag

    -R STR       SAM read group line in a format like '@RG\tID:foo\tSM:bar'

    -c           output CIGAR in PAF

    --cs[=STR]   output the cs tag; STR is 'short' (if absent) or 'long' [none]

    --MD         output the MD tag

    --eqx        write =/X CIGAR operators

    -Y           use soft clipping for supplementary alignments

    -t INT       number of threads [3]

    -K NUM       minibatch size for mapping [500M]

    --version    show version number

  Preset:

    -x STR       preset (always applied before other options; see minimap2.1 for details)

                 - map-pb/map-ont - PacBio CLR/Nanopore vs reference mapping

                 - map-hifi - PacBio HiFi reads vs reference mapping

                 - ava-pb/ava-ont - PacBio/Nanopore read overlap

                 - asm5/asm10/asm20 - asm-to-ref mapping, for ~0.1/1/5% sequence divergence

                 - splice/splice:hq - long-read/Pacbio-CCS spliced alignment

                 - sr - genomic short-read mapping

 

See `man ./minimap2.1' for detailed description of these and other advanced command-line options.

 

ビルド済みindex

2つのヒトゲノムアセンブリバージョンでのmodified indexが利用できる。

Nextcloud

 

テストラ

test/README.mdの通り進めます。

1、indexing

cd minimap2_index_modifier/tests/
#比較のため通常のindexとバリアントを反映させたindexを作る
#A normal
../minimap2 -ax sr -d test.mni test.fasta

#B modified index (このminimap2 folkでないと実行できない点に注意)
../minimap2 -ax sr -d test.modified.mni --vcf-file-with-variants test_long_chr1_not_the_same.vcf.gz test.fasta

=>diffコマンドで異なるか確認
diff test.mni test.modified.mni

 

2、mapping

../minimap2 -ax sr -t 4 test.modified.mni test.1.fq test.2.fq > mapped.sam

アラインメントフォーマットには変化はないが、バリアントを含む部位でのアラインメントが改善されている可能性がある。

 

レポジトリと論文より

  • 入力 VCF ファイルにリファレンスからのvariant が含まれていない場合、修正されたインデックスは通常のものと同じになる。
  • 既知の遺伝子変異を含むデータベースの例としては、dbSNP、gnomAD、1000 Genomes Project、Human Pangenome Reference Consortiumなどがある。これらのデータベースはすべて、線形リファレンス配列インデックスを修正するために使用できる。1000 Genomes Projectのデータセットには段階的な遺伝子型が含まれており、一緒に出現する遺伝子変異のみを考慮することができる。

 

コメント

パイプラインに影響を与えないため、最小限の手間でバリアントコールパイプラインを修正することが可能になります。気になるのは、図7のグラフで、わずかな差ではありますが、デフォルトindexのminimap2はrecallの指標でbwa mem2に負けていて、ラージコホートのバリアントを使ったmodified indexを使ったminimap2はbwa mem2を上回る点でしょうか。当時minimap2がショートリードマッピングでわずかに劣る可能性が指摘されてましたが、それが見えているのかと思いました

引用

Enhancing SNV identification in whole-genome sequencing data through the incorporation of known genetic variants into the minimap2 index

Egor Guguchkin, Artem Kasianov, Maksim Belenikin, Gaukhar Zobkova, Ekaterina Kosova, Vsevolod Makeev & Evgeny Karpulevich 

BMC Bioinformatics volume 25, Article number: 238 (2024) 

 

*1

デフォルトでは"../htslib"にある HTSlibに対してビルドされるので、取ってきたminimap2_index_modifierの直下でHTSlibライブラリをビルドした。

 

関連