Hi-CデータとロングリードDe Bruijnグラフの統合による効率的なテロメア間ゲノムアセンブリを行う Verkko2

2025/06/06 論文引用

　テロメア・ツー・テロメア・コンソーシアムは最近、ヒトゲノムの初の完全な配列を完成させた。最も複雑な繰り返しを解決するために、このプロジェクトは、長くて正確なPacBio HiFiとウルトラロングOxford Nanoporeシーケンスリードの半手動の組み合わせに依存した。その後、Verkkoアセンブラがこのプロセスを自動化し、二倍体ヒトゲノムの約半数の染色体を完全にアセンブルすることに成功した。しかし、Verkkoの最初のバージョンは計算コストが高く、典型的なヒトゲノムの全領域を解析することはできなかった。Verkko2は、より効率的なリード補正アルゴリズムを実装し、リピート解像とギャップクロージングを改善し、近接ライゲーションに基づくハプロタイプフェーズ決定とスキャフォールディングを導入し、複数のロングリードデータタイプをサポートした。これらの機能強化により、Verkkoは二倍体ヒトゲノムの全領域をアセンブルできるようになり、これには二倍体の短腕染色体や両性染色体も含まれる。これらの変更により、テロメア間のスキャフォールド数が2倍増加し、実行時間が4倍短縮され、アセンブルの正確性が向上した。19のヒトゲノムパネルにおいて、Verkko2は46本の完全な染色体のうち平均39本をスキャフォールドとしてアセンブルし、そのうち21本はギャップレス・コンティグとしてアセンブルされた。これらの改良により、テロメア間の比較やパンゲノミクスがスケールアップして可能になった。

インストール

レポジトリの指示通りcondaで環境を作って導入した。

Github

mamba create -n verkko2 -c conda-forge -c bioconda -c defaults verkko -y
conda activate verkko2

#Trioモードで実行するには、まずmerqury hapmerデータベースを生成してverkkoに渡す必要がある。

> verkko

usage: /home/kazu/mambaforge/envs/verkko/bin/verkko -d <output-directory> --hifi <hifi-reads ...> --nano <nanopore-reads ...>

MANDATORY PARAMETERS:

-d <output-directory> Directory to use for verkko intermediate and final results.

Will be created if needed.

--hifi <files ...> List of files containing PacBio HiFi reads.

Input reads can be any combination of FASTA/FASTQ,

uncompressed or gzip/bzip2/xz compressed. Any

number of files can be supplied; *.gz works.

--nano <files ...> List of files containing Oxford Nanopore reads.

Input reads can be any combination of FASTA/FASTQ/SAM/BAM,

uncompressed or gzip/bzip2/xz compressed. Any

number of files can be supplied; *.gz works.

ALGORITHM PARAMETERS:

--no-correction Do not perform Canu correction on the HiFi reads.

--no-nano Assemble without ONT data.

--hap-kmers h1 h2 type Use rukki to assign paths to haplotypes. 'h1' and 'h2

must be Meryl databases of homopolymer-compressed parental

kmers. 'type' must be 'trio', 'hic' or 'strandseq'.

--hic1 <files ...> List of files containing left hic reads.

--hic2 <files ...> List of files containing right hic reads.

Order of left and right files should be same, no interlaced files allowed.

Input reads can be any combination of FASTA/FASTQ/SAM/BAM,

uncompressed or gzip/bzip2/xz compressed. Any

number of files can be supplied; *.gz works.

--porec <files> List of files containing Pore-C reads.

Input reads can be any combination of FASTA/FASTQ/SAM/BAM,

uncompressed or gzip/bzip2/xz compressed. Any

number of files can be supplied; *.gz works.

--no-rdna-tangle Switch off option that helps to proceed large rDNA tangles which may connect multiple chromosomes.

--telomere-motif Switch to user-supplied telomere motif, automatically searched in fwd/revcomp version.

By default, the cannonical vertebrate 'CCCTAA' motif is used during Hi-C scaffolding.

--ref <path> Path to optional reference. It is used only as a guidance during scaffolding.

--uneven-depth Disable coverage-based heuristics in homozygous nodes detection for Hi-C/PoreC phasing.

--haplo-divergence Estimation on maximum divergence between haplotypes, is used only with Hi-C/PoreC data.

Should be increased for species with divergence significantly higher than in human. Default: 0.05, min 0, max 0.2

--screen <contaminant_name> </full/path/to/contaminant.fasta>

Screens contaminant from provided file. Multiple screen commands are allowed and are additive.

To screen typical contaminants of human assemblies one can use --screen-human-contaminants option with no parameters, which removes rDNA, mitochondria and EBV.

--paths <paths.gaf> Option required to generate fasta file from user-provided paths through assembly graph.

Requires an existing assembly specified by '--assembly <path_to_assembly>' and input reads.

Output is specified as '-d <output_dir>' (<output_dir> should not be equal to <path_to_assembly>).

The gaf file must be formatted as follows: 'name >utig4-1<utig4-2 HAPLOTYPE1' with one path per line,

where utig4-1 is in fwd orientation and utig4-2 is in reverse complement.

COMPUTATIONAL PARAMETERS:

--python <interpreter> Path or name of a python interpreter. Default: 'python'.

--perl <interpreter> Path of name of a perl interpreter. Default: 'perl'.

--mbg <path> Path to MBG. Default for all three

--graphaligner <path> Path to GraphAligner. one packaged with verkko.

--mashmap <path> Path to MashMap.

--seqtk <path> Path to seqtk.

--winnowmap <path> Path to Winnowmap.

--minimap <path> Path to minimap2.

--bwa <path> Path to BWA.

--samtools <path> Path to Samtools.

--cleanup Remove intermediate results.

--no-cleanup Retain intermediate results (default).

--local Run on the local machine (default).

--local-memory Specify the upper limit on memory to use, in GB, default 64

--local-cpus Specify the number of CPUs to use, default 'all'

--grid Enable Grid support (SGE, Slurm, PBS, and LSF are supported).

--snakeopts <string> Append snakemake options in "string" to the

snakemake command. Options MUST be quoted.

--sto-run Set resource limits for various stages.

--ovb-run Format: number-of-cpus memory-in-gb time-in-hours

--ovs-run --cns-run 8 32 2

--red-run

--mbg-run

--utg-run

--spl-run

--ali-run

--pop-run

--utp-run

--lay-run

--sub-run

--par-run

--cns-run

--ahc-run

--fhc-run

--shc-run

ADVANCED MODULE PARAMETERS (expert users):

HiFi read correction:

--correct-k-mer-size Set the k-mer size to use for finding overlaps (425)

--correct-k-mer-window Set the window size for sketching reads when finding overlaps (19)

--correct-k-mer-coverage Set the maximum count of kmer to use for finding overlaps (1000)

--correct-min-read-length Set the overall minimum read length (4000)

--correct-min-overlap-length Set the minimum overlap length (1000)

--correct-batch-size Set the RED batch size, in Mbp (30000)

(might also need to adjust --red-run)

MBG:

--base-k

--max-k

--window

--threads

--unitig-abundance

--hifi-coverage

ONT splitting:

--split-bases

--split-reads

--min-ont-length

GraphAligner:

--seed-min-length

--seed-max-length

--align-bandwidth

--score-fraction

--min-identity

--min-score

--end-clipping

--incompatible-cutoff

--max-trace

--haploid Better ONT alignment for haploid genomes

Verkko module path: /home/kazu/mambaforge/envs/verkko/lib/verkko/

No output directory (-d) set.

No PacBio HiFi reads (--hifi) supplied.

実行方法

VerkkoはSnakemakeワークフローとして実装されている。verkko.ymlファイルを作成するラッパースクリプトによって起動される。”--hifi” オプションで指定する高精度のロングリードは必須となっている。hifi パラメーターには、PacBio HiFi/Oxford Nanopore duplexを指定できる（複数も可）。”--nano”を使用して、ウルトラロングリードデータを含めることが強く推奨されている。また、フェーシング情報の指定も推奨されている。HERRO補正リードの場合は、補正リードを"--hifi"オプションで、未補正リードを"--nano"で指定する。

#最小設定だと、作業ディレクトリとHiFi fastqを指定する。
verkko -d workdir/ --hifi hifi.reads.fq.gz

-d Directory to use for verkko intermediate and final results. Will be created if needed.　
--hifi List of files containing PacBio HiFi reads. Input reads can be any combination of FASTA/FASTQ, uncompressed or gzip/bzip2/xz compressed. Any number of files can be supplied; *.gz works.
--nano List of files containing Oxford Nanopore reads. Input reads can be any combination of FASTA/FASTQ/SAM/BAM, uncompressed or gzip/bzip2/xz compressed. Any number of files can be supplied; *.gz works.

出力例

workdir/

最終的なアセンブリ結果は、assembly.fasta。最終的なグラフは、assembly.homopolymer-compressed.gfaで、それに付属するカバレッジファイルがassembly*csv。assembly.scfmapは最終配列assembly.fastaをグラフノードに変換したもの。フェーズ情報を提供した場合は、assembly.haplotype.fasta、assembly.colors.csv、asm/assembly.paths.tsvも保存される。後者の2つのファイルは、グラフの各ノードの位相情報から得られた色と、アセンブリの位相合わせに選択されたパスに関する情報を提供する。

（複数の）Hi-C illuminaリードペアを指定

verkko -d workdir/ --hifi hifi.reads.fq.gz  --hic1 *R1*fastq.gz --hic2 *R2*fastq.gz

--hic1 List of files containing left hic reads.
--hic2 List of files containing right hic reads.

Order of left and right files should be same, no interlaced files allowed. Input reads can be any combination of FASTA/FASTQ/SAM/BAM, uncompressed or gzip/bzip2/xz compressed. Any number of files can be supplied; *.gz works.

さらに（複数の）ONTウルトラロングリードも指定

verkko -d workdir/ --hifi hifi.reads.fq.gz  --hic1 *R1*fastq.gz --hic2 *R2*fastq.gz --nano *ont.fastq.gz \

--nano List of files containing Oxford Nanopore reads.

その他

クロマチンコンフォメーションキャプチャー（3C）とONTのロングリードシーケンスを利用して、3次元クロマチン構造を解析するPore-Cリードもオプションで使用できる（wiki）。
入力リードは、FASTA/FASTQ/SAM/BAMの任意の組み合わせ、非圧縮、またはgzip/bzip2/xz圧縮に対応している。
Hi-CやPore-Cのデータが提供された場合にだけ使う、独立したScaffoldingモジュールが用意されている。scaffoldingモジュールはアセンブリ中のテロメア位置（seqtk teloで検出）を使用するので、対象の生物種が脊椎動物（CCCTAA）とは異なるテロメアリピートモチーフを持っている場合、-telomere-motifオプションで指定する。
scaffoldingモジュールはアセンブリの2倍体構造に依存しているため、--haploidオプションとは互換性がない。このような場合はスタンドアロンscaffolderのYaHSが推奨されている。
--refオプションで同じ種または近縁種の別の個体のゲノムを提供し、スキャフォールディングのガイダンスとして使用できる（アセンブリのガイドには使用されない。スキャホールディングにのみ使用される）。
-screen オプションを使用してアセンブリから汚染物質をフィルタリングできる。複数の --screen <contaminant_N_name> <contaminant_N_sequence.fasta> コマンドを使用して、任意の数のターゲットを指定可能。各汚染物質について、verkko はメインアセンブリ出力からターゲットにマッチする配列をすべて削除する。スクリーニング情報を提供した場合、出力ディレクトリにはassembly.exampleN.fastaとassembly.exampleN.exemplar.fastaファイルが作成される。このファイルには、汚染物質と判断されてアセンブリから削除されたすべての配列と、（可能であれば）環状化された大文字の配列が保存されている（開始/終了の自己相似性を指標として環状化し、カバレッジとともに判断の指標としている）。
ヒトアセンブリの典型的なコンタミナントについては、特別なオプション --screen-human-contaminants (バージョン 2.2 以前では --screen human) が用意されている。これはパラメータを必要としない。”-screen rDNA rdna.fasta --screen mitochondria mito.fasta --screen EBV ebv.fasta ”のショートカットと同義。
テストのために大腸菌K12のデータセットが用意されている（レポジトリREADMEの下の方）。

引用

Verkko2: Integrating proximity ligation data with long-read De Bruijn graphs for efficient telomere-to-telomere genome assembly, phasing, and scaffolding

Dmitry Antipov, Mikko Rautiainen, Sergey Nurk, Brian P. Walenz, Steven J. Solar, Adam M. Phillippy, Sergey Koren

bioRxiv, Posted December 22, 2024.

Genome Res. 2025 May 19:gr.280383.124. doi: 10.1101/gr.280383.124. Online ahead of print.

Verkko2 integrates proximity ligation data with long-read De Bruijn graphs for efficient telomere-to-telomere genome assembly, phasing, and scaffolding

Dmitry Antipov, Mikko Rautiainen, Sergey Nurk, Brian P Walenz, Steven J Solar, Adam M Phillippy, Sergey Koren