macでインフォマティクス

macでインフォマティクス

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

(ヒトゲノム)高速かつ精度の高いロングリードのSVコーラー cuteSV

2021 5/6 インストール手順、help、コマンド更新

 

 構造変化(SV)とは、欠失、挿入、逆位、重複、転座などのゲノムリアレンジメントで、その大きさが50 bpを超えるものを指す。ヒトゲノム上で最大のdivergencesとして、SV はヒトの疾患(遺伝性疾患やガンなど)、進化(遺伝子欠損やトランスポゾン活性など)、遺伝子制御、その他の表現型(交配や本質的な生殖隔離など)と密接に関連している。

 ショートリードベースのSVコールアプローチを開発するための努力がなされてきた。それらの多くは、リードデプス、不一致リードペア、スプリットリードアラインメント、ローカルアセンブリ、またはそれらの組み合わせなどの方法を用いており、1000 Genomes Projectのような大規模なゲノム研究において重要な役割を果たしてきた。しかし、これらのツールは比較的短いリード長であるため、高感度なSV検出には限界があり、偽陽性も存在する。

 Pacific Bioscience (PacBio)や Oxford Nanopore Technology (ONT)などのロングリードシークエンシング技術の急速な発展に伴い、ロングレンジスパニング情報は、より高分解能でより包括的に SV を検出する機会を提供している。しかし、高いシーケンシングエラー率(通常5~20%)とリード長(平均10kbp以上)に対応するためには、新たな計算アプローチが必要である。これまでの研究では、主に2種類のアプローチが採用されてきた。

 デノボアセンブリーベースのアプローチ は、より長いゲノム配列(コンティグやスキャホールド)にリードをアセンブリし、アセンブリした配列とリファレンス配列とのアラインメントからSVを発見することを目的としている。このようなアプローチは、リードアライメントに基づくアプローチに比べて、リファレンスの影響を受けにくく、特にリードアライメントのアーチファクトがない。しかし、これらのアプローチは通常計算量が多く、大規模なゲノムのハプロタイプ配列の再構成にはまだいくつかの困難がある。

 リードアライメントに基づくアプローチは、リードをリファレンスに対して直接アライメントし、アライメント結果を解析することでSVを検出する。このようなアプローチは、感度を欠くことなく計算リソースに対してより費用対効果が高く、ロングリードベースのSVコーラーで広く使用されている。いくつかのリードアライメントベースのSVコーラーが提案されており、PB-Honey [ref.31]、SMRT-SV [ref.32]、Sniffles [ref.33]、PBSV (https://github.com/PacificBiosciences/pbsv)、およびSVIM [ref.34]などがある。これらの手法では、発散性の高いアラインメントを持つ局所的なゲノム領域の同定、クリッピングされたリード部分の局所的なアセンブリと再アラインメント、SV-スパンニングシグネチャクラスタリングなど、リードアラインメントによって暗示されるSVの証拠を見つけるために様々な方法を使用している。さらに、通常、BLASR [ref.36]、NGMLR [ref.33]、Minimap2 [ref.37]、PBMM2 (https://github.com/PacificBiosciences/pbmm2) などの最先端のロングリードアライナーがリードアライメントに使用されていた。

 しかし、リードアライメントに基づくSVコールは、依然として非自明である。シーケンシングエラーが多く、複雑なSVが存在する状況下では、SVブレークポイント周辺のリードのアライメントはキメラ的で不均質であり、通常は感度や精度が低い。そのため、リードアラインメントから得られるSVシグネチャは非常に複雑であり、様々な種類のSVに対して感度の高い検出を行うためには、それらを収集・解析することが困難である。主に、最先端のツールでは、以下のような技術的な課題がある。(1) 全体的に感度がまだ満足のいくものではない。(2) いくつかのアプローチ(rMETL [ref.38]、rCANID [ref.39]、npInv [ref.40]など)は、特定の設計のために、サブセットまたは特定のクラスのSVしか検出できない。(3) いくつかのアプローチ(PBSVやSMRT-SVなど)は、まだ時間がかかり、スケーリング性能が良くないため、多くの大規模データセットには適していない可能性がある;そして(4) いくつかのアプローチ(SMRT-SVやPB-Honeyなど)は、1種類のシーケンシングデータしかサポートしていない(例えばPacBioのリードのみ)ものもある。このような欠点が、ロングリードシーケンシングデータを広く利用する上でのボトルネックとなっている。 

 ここでは、いくつかの利点を持つ多用途なSV検出アプローチであるcuteSVを紹介する。(1) cuteSVは、最先端のSVコーラーよりも高いSV検出率を持っている。特に、カバレッジの低いデータに対しても、精度を落とさずに高感度に検出することができる。(2) cuteSVは、主流のロングリードシーケンシングプラットフォームで作成された様々なエラー率のデータセットをサポートしており、様々なタイプのSV(欠失、挿入、重複、逆位、転座を含む)を検出することができる。(3) cuteSVは、RAMの使用量が少なく、最先端のアプローチと同等の速度を持っている。さらに重要なことは、CPUのスレッド数に応じてほぼ直線的に高速化できるスケーラビリティを持っている。これらの特徴により、cuteSVは大規模なデータ解析に適しており、最先端のゲノミクス研究への可能性を秘めている。

cuteSVの概要
ソートされたBAMファイルを入力として、cuteSVは、アライメントされたリード中の大きな挿入/欠失や分割アライメントをSVシグネチャクラスターとして抽出し、解析してSVをコールする。このアプローチには大きく分けて以下の3つのステップがある。

ステップ1:cuteSVは複数のシグネチャ抽出手法を用いて、様々なタイプのSVのシグネチャを網羅的に収集する。さらに、挿入と欠失をヒューリスティックに組み合わせて、脆弱なアライメントから本物のSVの証拠を回収する。

ステップ2: cuteSVは、特別に設計されたクラスタリングと精緻化のアプローチを用いて、キメラアラインメントされたリードを局所領域にクラスタリングし、さらにクラスタを精緻化して、SVのシグネチャヘテロ接合型SVと正確に区別する。

ステップ3:cuteSVは、いくつかのオーダーメイドのルールを使用して、洗練されたSVシグネチャクラスターに基づいてSVコールとジェノタイピングを実行する。

 

 

benchmark

https://github.com/tjiangHIT/sv-benchmark

 

cuteSV

 

Blog: cuteSV

https://nanoporetech.com/about-us/news/blog-cutesv-powerful-tool-uncover-full-spectrum-genomic-structural-variants

 

インストール

依存

  • python3
  • pysam
  • Biopython
  • cigar
  • numpy
  • pyvcf

Github

#bioconda (link)
mamba create -n cutesv -y
conda activate cutesv
mamba install -c bioconda cutesv -y

#pip
pip install cuteSV

#git
git clone https://github.com/tjiangHIT/cuteSV.git && cd cuteSV/ && python setup.py install

cuteSV -h  # version1.0.10

$ cuteSV -h

usage: cuteSV [-h] [--version] [-t THREADS] [-b BATCHES] [-S SAMPLE] [--retain_work_dir] [--report_readid] [-p MAX_SPLIT_PARTS] [-q MIN_MAPQ] [-r MIN_READ_LEN] [-md MERGE_DEL_THRESHOLD] [-mi MERGE_INS_THRESHOLD]

              [-s MIN_SUPPORT] [-l MIN_SIZE] [-L MAX_SIZE] [-sl MIN_SIGLENGTH] [--genotype] [--gt_round GT_ROUND] [-Ivcf IVCF] [--max_cluster_bias_INS MAX_CLUSTER_BIAS_INS] [--diff_ratio_merging_INS DIFF_RATIO_MERGING_INS]

              [--max_cluster_bias_DEL MAX_CLUSTER_BIAS_DEL] [--diff_ratio_merging_DEL DIFF_RATIO_MERGING_DEL] [--max_cluster_bias_INV MAX_CLUSTER_BIAS_INV] [--max_cluster_bias_DUP MAX_CLUSTER_BIAS_DUP]

              [--max_cluster_bias_TRA MAX_CLUSTER_BIAS_TRA] [--diff_ratio_filtering_TRA DIFF_RATIO_FILTERING_TRA]

              [BAM] reference output work_dir

 

 

Current version: v1.0.10

Author: Tao Jiang

Contact: tjiang@hit.edu.cn

 

If you use cuteSV in your work, please cite:

Jiang T et al. Long-read-based human genomic structural variation detection with cuteSV. 

Genome Biol 21,189(2020). https://doi.org/10.1186/s13059-020-02107-y

 

Suggestions:

 

For PacBio CLR data:

--max_cluster_bias_INS 100

--diff_ratio_merging_INS 0.3

--max_cluster_bias_DEL 200

--diff_ratio_merging_DEL 0.5

 

For PacBio CCS(HIFI) data:

--max_cluster_bias_INS 1000

--diff_ratio_merging_INS 0.9

--max_cluster_bias_DEL 1000

--diff_ratio_merging_DEL 0.5

 

For ONT data:

--max_cluster_bias_INS 100

--diff_ratio_merging_INS 0.3

--max_cluster_bias_DEL 100

--diff_ratio_merging_DEL 0.3

 

 

 

positional arguments:

  [BAM]                 Sorted .bam file form NGMLR or Minimap2.

  reference             The reference genome in fasta format.

  output                Output VCF format file.

  work_dir              Work-directory for distributed jobs

 

optional arguments:

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

  --version, -v         show program's version number and exit

  -t THREADS, --threads THREADS

                        Number of threads to use.[16]

  -b BATCHES, --batches BATCHES

                        Batch of genome segmentation interval.[10000000]

  -S SAMPLE, --sample SAMPLE

                        Sample name/id

  --retain_work_dir     Enable to retain temporary folder and files.

  --report_readid       Enable to report supporting read ids for each SV.

 

Collection of SV signatures:

  -p MAX_SPLIT_PARTS, --max_split_parts MAX_SPLIT_PARTS

                        Maximum number of split segments a read may be aligned before it is ignored.[7]

  -q MIN_MAPQ, --min_mapq MIN_MAPQ

                        Minimum mapping quality value of alignment to be taken into account.[20]

  -r MIN_READ_LEN, --min_read_len MIN_READ_LEN

                        Ignores reads that only report alignments with not longer than bp.[500]

  -md MERGE_DEL_THRESHOLD, --merge_del_threshold MERGE_DEL_THRESHOLD

                        Maximum distance of deletion signals to be merged. In our paper, I used -md 500 to process HG002 real human sample data.[0]

  -mi MERGE_INS_THRESHOLD, --merge_ins_threshold MERGE_INS_THRESHOLD

                        Maximum distance of insertion signals to be merged. In our paper, I used -mi 500 to process HG002 real human sample data.[100]

 

Generation of SV clusters:

  -s MIN_SUPPORT, --min_support MIN_SUPPORT

                        Minimum number of reads that support a SV to be reported.[10]

  -l MIN_SIZE, --min_size MIN_SIZE

                        Minimum size of SV to be reported.[30]

  -L MAX_SIZE, --max_size MAX_SIZE

                        Maximum size of SV to be reported.[100000]

  -sl MIN_SIGLENGTH, --min_siglength MIN_SIGLENGTH

                        Minimum length of SV signal to be extracted.[10]

 

Computing genotypes:

  --genotype            Enable to generate genotypes.

  --gt_round GT_ROUND   Maximum round of iteration for alignments searching if perform genotyping.[500]

 

Force calling:

  -Ivcf IVCF            Optional given vcf file. Enable to perform force calling. [NULL]

 

Advanced:

  --max_cluster_bias_INS MAX_CLUSTER_BIAS_INS

                        Maximum distance to cluster read together for insertion.[100]

  --diff_ratio_merging_INS DIFF_RATIO_MERGING_INS

                        Do not merge breakpoints with basepair identity more than [0.3] for insertion.

  --max_cluster_bias_DEL MAX_CLUSTER_BIAS_DEL

                        Maximum distance to cluster read together for deletion.[200]

  --diff_ratio_merging_DEL DIFF_RATIO_MERGING_DEL

                        Do not merge breakpoints with basepair identity more than [0.5] for deletion.

  --max_cluster_bias_INV MAX_CLUSTER_BIAS_INV

                        Maximum distance to cluster read together for inversion.[500]

  --max_cluster_bias_DUP MAX_CLUSTER_BIAS_DUP

                        Maximum distance to cluster read together for duplication.[500]

  --max_cluster_bias_TRA MAX_CLUSTER_BIAS_TRA

                        Maximum distance to cluster read together for translocation.[50]

  --diff_ratio_filtering_TRA DIFF_RATIO_FILTERING_TRA

                        Filter breakpoints with basepair identity less than [0.6] for translocation.

 

 

実行方法

bamファイルとリファレンスファイルを指定する。

mkdir work_dir
cuteSV -t 20 input.bam ref.fasta output.vcf work_dir
  • [BAM]    Sorted .bam file form NGMLR or Minimap2.
  • reference    The reference genome in fasta format.
  • output    Output VCF format file.
  • work_dir    Work-directory for distributed jobs

vcf4.2 formatのSVコールが出力される。

 

 

#2022/04/18 スモールゲノムアセンブリのエラーチェックに使う。

ハプロイドゲノムアセンブリで構造変異のアセンブリエラーがあれば、そのポジションはアラインされたほぼ全てのリードがエラーのシグナルをは発生する。そのため、サポートするリードの数を増やすことで、アセンブリのエラーだけを特異的に呼び出すことができる。

ロングリードは予めゲノムの100x程度のリードデプスに調整しておく。ONTのリードは長すぎるとエラー率が上がる傾向があるので、クオリティトリミングした上で、5kb-50kbくらいの長さのリードを取り出して、ゲノムのx100程度にサブサンプリングしておく。

#1 エラープロファイルを視覚化してチェック
fastp -i ONT.fq.gz

#2 サイズ選択
seqkit seq -m 5000 -M 50000 ONT.fq.gz |gzip - > ONT_5000-50000.fg.gz

#3 サンプリング
seqkit sample -p <x percent> ONT_5000-50000.fg.gz |gzip - > ONT_5000-50000_x100.fg.gz

#4 マッピング
minimap2 -ax map-ont -t 20 assembly.fasta ONT_5000-50000_x100.fg.gz |samtools sort -@ 8 - > ONT.bam

#4 任意ステップ;必要ならMAPQが0のリードを除く。ロングリードでもMAPQが0のリードはかなり出てきて、そのようなリードがSVとしてコールされることも多い
samtools view -@ 2 -bSq ONT.bam > ONT_filtered.bam
samtools index -@ 4 ONT_filtered.bam

#5 SVコール
mkdir work_dir
cuteSV -t 20 -s 50 -r 1000 -l 10 ONT.bam assembly.fasta output.vcf work_dir
  • -r     Ignores reads that only report alignments with not longer than bp.[500]
  • -s     Minimum number of reads that support a SV to be reported.[10]
  • -l      Minimum size of SV to be reported.[30]

 

閾値は不明だが、リードデプスがゲノムの100xのリードで試す限りは"-s 30-50"程度に調整するとシグナル/ノイズ比が良さそうだった。でフォルダと、サイズが30bp以上のIndelエラーが呼び出される。逆位のエラーも呼び出されるが、アラインメントを見てもわかりにくいので注意(VCFを良く見ること)。

引用

Long-read-based human genomic structural variation detection with cuteSV

Tao Jiang, Yongzhuang Liu, Yue Jiang, Junyi Li, Yan Gao, Zhe Cui, Yadong Liu, Bo Liu & Yadong Wang
Genome Biology volume 21, Article number: 189 (2020)