ロングリードからphased local assemblyにより2倍体ゲノムの構造変異とindelを正確に同定する VolcanoSV

構造変異（SV）はヒトゲノムの多様性に大きく寄与し、精密医療において重要な役割を果たす。一分子ロングリードシーケンシング技術の進歩はSV検出に画期的な資源を提供するものの、SVの切断点と配列を正確かつ堅牢に特定することは依然として困難である。本著者らは、リファレンスゲノムとローカルデノボアセンブリの両方を活用してフェージングされた二倍体アセンブリを生成する革新的なハイブリッドSV検出パイプライン「VolcanoSV」を導入する。VolcanoSVはphased SNPsとユニークk-mer類似性解析を用い、精密なhaplotype-resolved SVの発見を可能にする。SNP、スモールindel、あらゆるタイプのSVを包括する包括的な遺伝子マップ構築に優れており、ヒトゲノム研究に最適である。広範な実験により、VolcanoSVが挿入・欠失SVの検出において最先端のアセンブリベースツールを凌駕し、低カバレッジ（10x）データセットを含む多様なデータセット全体で、優れた再現率、精度、F1スコア、遺伝子型正確性を示すことが実証された。VolcanoSVは、シミュレーションデータと実際のガンデータの両方において、転座、重複、逆位を含む複雑なSVの同定においてアセンブリベースのツールを上回る。さらに、VolcanoSVは様々な評価パラメータに対して頑健であり、ブレークポイントとSV配列を正確に同定する。

インストール

ubuntu24に以下の手順でインストールした。

git clone --recurse-submodules https://github.com/maiziezhoulab/VolcanoSV.git
cd VolcanoSV
sh Install.sh #注；jarファイルのダウンロードのみ（link）
mamba env create -f requirement.yaml
conda activate volcanosv

$ python bin/VolcanoSV-asm/volcanosv-asm.py -h

usage: use "python3 volcanosv-asm.py --help" for more information

optional arguments:

-h, --help show this help message and exit

--bam_file BAM_FILE, -bam BAM_FILE

--output_dir OUTPUT_DIR, -o OUTPUT_DIR

--reference REFERENCE, -ref REFERENCE

--chrnum {1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22}, -chr {1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22}

--assembler {wtdbg2,canu,miniasm,shasta,nextdenovo,hifiasm,hicanu,flye}, -asm {wtdbg2,canu,miniasm,shasta,nextdenovo,hifiasm,hicanu,flye}

--data_type {CLR,ONT,Hifi}, -dtype {CLR,ONT,Hifi}

--pacbio_subtype {CLR-rs,CLR-sq}, -pb {CLR-rs,CLR-sq}

must provide when using wtdbg2 on CLR data (default: None)

--shasta_ont_config {Nanopore-OldGuppy-Sep2020}, -shacon {Nanopore-OldGuppy-Sep2020}

must provide when using shasta (default: None)

--n_thread N_THREAD, -t N_THREAD

number of threads (default: 10)

--prefix PREFIX, -px PREFIX

file prefix in the output folder (default: Sample)

--clean, -cl

> python bin/VolcanoSV-vc/Large_INDEL/volcanosv-vc-large-indel.py -h

usage: use "python3 volcanosv-vc-large-indel.py --help" for more information

optional arguments:

-h, --help show this help message and exit

--input_dir INPUT_DIR, -i INPUT_DIR

--output_dir OUTPUT_DIR, -o OUTPUT_DIR

--data_type DATA_TYPE, -dtype DATA_TYPE

Hifi;CLR;ONT (default: None)

--bam_file BAM_FILE, -bam BAM_FILE

reads bam file for reads signature extraction (default: None)

--reference REFERENCE, -ref REFERENCE

only needed when presig is not provided (default: None)

--read_signature_dir READ_SIGNATURE_DIR, -rdsig READ_SIGNATURE_DIR

pre-extracted reads signatures; optional; if not provided, will generate one using read bamfile (default: None)

--pre_cutesig PRE_CUTESIG, -presig PRE_CUTESIG

pre-extracted cutesv signature directory;optional; if not provided, will generate a new one (default: None)

--chr_num CHR_NUM, -chr CHR_NUM

chrmosome number; Optional; if not provided, will assume input_dir contain chr1-chr22 results (default: None)

--n_thread N_THREAD, -t N_THREAD

number of threads (default: 11)

--n_thread_align N_THREAD_ALIGN, -ta N_THREAD_ALIGN

number of threads for contig alignment (default: 10)

--mem_per_thread MEM_PER_THREAD, -mempt MEM_PER_THREAD

Set maximum memory per thread for alignment; suffix K/M/G recognized; default = 768M (default: 768M)

--prefix PREFIX, -px PREFIX

file prefix in the output folder (default: Sample)

> python bin/VolcanoSV-vc/Complex_SV/volcanosv-vc-complex-sv.py -h

usage: use "python3 volcanosv-vc-complex-sv.py --help" for more information

optional arguments:

-h, --help show this help message and exit

--input_dir INPUT_DIR, -i INPUT_DIR

--indelvcf INDELVCF, -vcf INDELVCF

--bam_file BAM_FILE, -bam BAM_FILE

--reference REFERENCE, -ref REFERENCE

--data_type {Hifi,CLR,ONT}, -dtype {Hifi,CLR,ONT}

--output_dir OUTPUT_DIR, -o OUTPUT_DIR

--n_thread N_THREAD, -t N_THREAD

--prefix PREFIX, -px PREFIX

file prefix in the output folder (default: Sample)

> python bin/VolcanoSV-vc/Small_INDEL/volcanosv-vc-small-indel.py -h

usage: use "python3 volcanosv-vc-small-indel.py --help" for more information

optional arguments:

-h, --help show this help message and exit

--input_dir INPUT_DIR, -i INPUT_DIR

--bam_file BAM_FILE, -bam BAM_FILE

--output_dir OUTPUT_DIR, -o OUTPUT_DIR

--reference REFERENCE, -ref REFERENCE

--bedfile BEDFILE, -bed BEDFILE

optional; a high confidence bed file (default: None)

--region REGION, -r REGION

optional; exmaple: chr21:2000000-2100000 (default: None)

--n_thread N_THREAD, -t N_THREAD

number pf threads (default: 50)

--kmer_size KMER_SIZE, -k KMER_SIZE

kmer size (default: 15)

--ratio RATIO, -rt RATIO

maximum bad kmer ratio (default: 0.3)

--min_support MIN_SUPPORT, -ms MIN_SUPPORT

maximum support for bad kmer (default: 5)

--restart, -rs restart mode; assume there is kmer support file already. (default: False)

--eval, -e

--prefix PREFIX, -px PREFIX

file prefix in the output folder (default: Sample)

テストデータ

VolcanoSVアセンブリは染色体単位で動作するよう設計されている。染色体1個分のデータがテスト用にzenodoにアップされている｡具体的には､リファレンスgenome.faとONT, CLR, HiFiの3種類のbam､VCF (Single chromosomeモードのラン結果)となる｡リードはhg19リファレンスにアラインされている｡

zenodo: https://zenodo.org/records/10520476

自身のデータにVolcanoSVを適用する前に、テストデータを実行することが強く推奨されている｡zenodoに公開されている染色体サブセットのデモとは別に､全ゲノムのテストデータ(WGSモードのデモデータ)も公開されている｡

> wget https://cf.10xgenomics.com/supp/genome/refdata-hg19-2.1.0.tar.gz

> tar -xzvf refdata-hg19-2.1.0.tar.gz

1､ Single chromosomeモード - VolcanoSV Assembly

bamファイル、リファレス、リードのデータタイプのほか、任意でローカルアセンブリに使うロングリードアセンブラを指定する。アセンブラは､Hifiデータについてはhifiasmとhicanumを選ぶ｡その他のアセンブラはCLRおよびONTデータ用に用意されている。

レポジトリでは、HiFiリードの実行例として以下のコマンドが記載されている (デモデータをダウンロードして使用)｡アセンブラはHiFiasmを指定｡

cd VolcanoSV/
python3 bin/VolcanoSV-asm/volcanosv-asm.py \
-bam Hifi_L2_hg19_minimap2_chr10.bam \
-o volcanosv_asm_output -ref genome.fa \
-t 10 -chr 10 -dtype Hifi -px Hifi_L2 -asm hifiasm

-bam BAM_FILE
-o OUTPUT_DIR
-ref REFERENCE
-chr {1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22}
-asm {wtdbg2,canu,miniasm,shasta,nextdenovo,hifiasm,hicanu,flye}
-dtype {CLR, ONT, Hifi}
-pb {CLR-rs,CLR-sq} must provide when using wtdbg2 on CLR data (default: None)
-t number of threads (default: 10)

小さな染色体1個でも実行にはある程度時間がかかる｡5995WXで10スレッド指定だと1時間半ほどかかった(他の計算をバックグラウンドで実行中のため､より多くの時間がかかった可能性あり)｡

出力例

volcanosv_asm_output/chr10/

最終コンティグはchr10/assembly/final_contigs/Hifi_L2_final_contigs.fa となる。

volcanosv-asmパイプラインが正常に実行された場合、最終コンティグファイルのサイズは zenodo の Hifi_L2_contigs.fa(299,093,441 bp) とほぼ同等になる。

2､ Single chromosomeモード - VolcanoSV Assembly（ハイブリッドモード）

VolcanoSVでは､複数アセンブラを使うハイブリッドモードも提供している｡これは､アセンブラによってセグメンタル重複（SD）やその他の複雑な領域が豊富な領域をアセンブルする能力には差があり､したがって、ゲノム領域ごとに異なるアセンブラを利用することで利点が得られる可能性があるためと説明されている。

ハイブリッドモードのランでは、オプションで「in-BED」アセンブラと「out-BED」アセンブラを指定する｡また､BEDファイルを指定する｡このBEDファイルと重複するフェーズブロックはin-BEDアセンブラでアセンブルされ、残りの領域はout-BEDアセンブラでアセンブルされる｡

cd VolcanoSV/
python3 bin/VolcanoSV-asm/volcanosv-asm_hybrid.py \
-bam Hifi_L2_hg19_minimap2_chr10.bam \
-o volcanosv_asm_output \
-ref genome.fa \
-bed segdups.bed \
--inbed_assembler hicanu \
--outbed_assembler hifiasm \
-t 10 \
-chr 10 \
-dtype Hifi \
-px Hifi_L2

3､Single chromosomeモード - Large Indel detection (VolcanoSV-vc)

Single chromosomeモードの結果(1 or 2)からlarge indelをコールする｡1か2の出力ディレクトリを指定する｡

cd VolcanoSV/
python3 bin/VolcanoSV-vc/Large_INDEL/volcanosv-vc-large-indel.py \
-i volcanosv_asm_output/ \
-o volcanosv_large_indel_output/ \
-dtype Hifi \
-bam Hifi_L2_hg19_minimap2_chr10.bam \
-ref genome.fa \
-chr 10 -t 10 \
-px Hifi_L2

入力ディレクトリはvolcanoSV-asmの出力ディレクトリ。このコマンドは単一染色体モードと全ゲノムシーケンスモードの両方に対応していて､オプション"-chr <NUM>"が指定された場合は単一染色体モードで実行され、指定されない場合は入力ディレクトリにchr1-chr22のコンティグが含まれていると仮定され全ゲノムシーケンスモードで実行される。プレフィックスはvolcanoSV-asmで設定したものと一致させる必要がある。

テスト時は数分で終了した｡

出力例

volcanosv_large_indel_output/chr10/

出力VCFは<ouput_folder>/volcanosv_large_indel.vcfに生成される。

4, Single chromosome mode Small Indel detection (VolcanoSV-vc)

cd VolcanoSV/
python3 bin/VolcanoSV-vc/Small_INDEL/volcanosv-vc-small-indel.py \
-i volcanosv_asm_output/ \
-o volcanosv_small_indel \
-bam Hifi_L2_hg19_minimap2_chr10.bam \
-ref genome.fa \
-r chr10 \
-t 30 \
-px Hifi_L2

LongshotのSNPコール結果が最終的なSNPコールとして採用される。フェージングされたSNP VCFファイルが生成される｡テスト時は10分ちょっとで終了した｡

5, WGSモード - VolcanoSV Assembly (VolcanoSV-asm)

VolcanoSVアセンブリは染色体単位で動作するよう設計されている。そのため、全ゲノムのSVをコールする場合、常染色体すべて順番か分散ジョブシステムで効率的にアセンブリを実行して､染色体すべてのphased assembly結果を得る(出力ディレクトリにchr-22のサブディレクトリが存在する状態にする)｡それから､上記のコマンドを-chr <NUM>なしで実行してlarge indel, small indel､染色体間の転座などの複雑なSVをコールする｡詳細はレポジトリ参照｡

レポジトリより

VolcanoSVアセンブリパイプライン(bin/VolcanoSV-asm/volcanosv-asm.py)は染色体単位での実行を想定して設計されている。hifiasm、Flye、wtdbg2、miniasm、Shasta、NextDenovo、Hicanu/Canuなど、複数の最先端アセンブラを統合していて､ユーザーはニーズに応じて適切なアセンブラを選択できる。
既に全アセンブラの実行可能版が含まれているため、個別にインストールする必要はない
VolcanoSVは常染色体専用に設計されている。性染色体のフェーシング手順は大きく異なるため
参照ファイルのコンティグ名は「chr1, chr2, chr3 ...」形式とする｡説明フィールドを付加しないこと
SVコール（特にcomplex SV や large indel）には十分な計算資源が必要
転座検出にはWGS BAMファイルが必要で、したがってSingle chromosomeモードは意味をなさない｡転座検出はWGSモードでのみ可能
VolcanoSVアセンブリは染色体単位で動作するよう設計されている。複数のジョブ提出をサポートする分散コンピューティングシステムを利用できる場合は、染色体ごとに1つのジョブを提出し、それらを並行して実行することを推奨する。ジョブスクリプトを作成するテンプレートがレポジトリに用意されている｡分散コンピューティングシステムがない場合は、異なる染色体を順次実行するためのforループを記述することを推奨している(レポジトリにコマンド例)
出力が期待値から大きくずれる場合は、アセンブリパラメータや入力データ品質（リードの質、マッピング精度、重複領域など）を見直す
HiFiデータではなくCLRデータまたはONTデータを使用する場合、引数「dtype」を対応するデータタイプ（CLR/ONT）に変更する