リファレンスゲノム配列に対するリードのアライメントは、次世代シーケンサー（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
#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
２つのヒトゲノムアセンブリバージョンでのmodified indexが利用できる。
テストラン
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
２、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のデータセットには段階的な遺伝子型が含まれており、一緒に出現する遺伝子変異のみを考慮することができる。
コメント
パイプラインに影響を与えないため、最小限の手間でバリアントコールパイプラインを修正することが可能になります。気になるのは、図７のグラフで、わずかな差ではありますが、デフォルト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ライブラリをビルドした。
関連
