macでインフォマティクス

macでインフォマティクス

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

2倍体ゲノムアセンブリからHaplotigsを追い出してPrimary contigsを出力する Purge Haplotigs

 2020 7/11 図追加

2020 7/13  タイトル修正

2020 7/15 コメント追記

2021 12/23 コメント追加

2022/09/18 インストール手順修正

 

 第三世代の1分子シーケンシングにおける最近の進歩は、非常に高いレベルの連続性と完全性を持つde novoゲノムアセンブリを可能にした。さらに、最近の「diploid aware」なゲノムアセンブラの進歩により、高度にヘテロ接合した二倍体ゲノムアセンブラの品質が大幅に向上した。FALCONやCanuなどの二倍体対応アセンブラは、二倍体ゲノムのハプロタイプを融合した表現を生成することができる。また、FALCON UnzipやSupernovaなどのアセンブラの中には、さらに進化して、両方の親対立遺伝子を別々に表現した大規模なフェーズブロックを生成するものもある。FALCON Unzipアセンブリでは、アセンブリグラフ上でフェーズングが行われ、「一次コンティグ」(ハプロイドアセンブリ)と関連する「haplotigs」が生成され、これらの一次コンティグと二次ハプロチグのunionからなる二倍体アセンブリが生成される。

 理想的なハプロイド表現(プライマリコンティグ)は、2つのハプロームのすべてのヘテロ接合領域と、両方のハプロームのすべてのヘミ接合領域の1つの対立遺伝子コピーで構成されている。これにより、どちらかのハプロームに含まれる領域が、ハプロイド表現の一つの位置に完全にアラインすることが保証される。セカンダリハプロチグには、両方のハプロームに見られるヘテロ接合領域の2つの対立遺伝子コピーのうちの1つが含まれていなければならない;この点で、ハプロチグはハプロイド表現の段階的な情報として機能する。

 非常に高いヘテロ接合性を持つ領域は、de novoゲノムアセンブリではまだ問題がある。このような状況では、一組の対立遺伝子配列がヌクレオチド多様性のある閾値を超えると、ほとんどのアルゴリズムは、期待される単一のハプロタイプ融合コンティグではなく、これらの領域を別々のコンティグとしてアセンブルする 。結果、ハプロイドゲノムのサイズよりも著しく大きいアセンブリとなり、ハプロイドアセンブリにおけるこれらの対立遺伝子コンティグの存在は、下流の解析において問題となる。二倍体アセンブリを作製する場合、両方の対立遺伝子が存在する可能性があるが、対立遺伝子コンティグのペアを特定するための手順が依然として必要である。 

 この問題に対処するためにいくつかのツールが試みられてきた。HaploMerger2ツールキットとRedundansアセンブリパイプラインは、ショートリード配列からハプロタイプ融合アセンブリを生成するように設計されている。しかし、カバレッジの深さを考慮せずにコンティグ同士のアラインメントのみに基づいてコンティグを自動的に除去すると、反復的なコンティグやパラロジカルなコンティグが過剰にパージ(排除)されてしまう可能性がある。さらに、ハプロタイプ配列を分解して段階的なアセンブリを作成することが有利であることが証明されている。ロングリードのアセンブリで使用できるスクリプトには、配列のアラインメントを使用してホモログを特定し、手動でのキュレーションを支援する get_homologs.pyや、遺伝子アノテーションを使用して同調領域を一致させる HomolContigsByAnnotation などがある。それぞれに独自の長所と短所があるが、どちらもコンティグの再配置をユーザーが手動で行わなければならないという問題がある。

 本研究の目的は、1分子ロングリードシーケンシング技術を用いて作製されたアセンブリ内の対立遺伝子コンティグを迅速かつ自動的に特定し、再割り当てできる新しいパイプラインを開発することであった。Purge Haplotigsはインストールが簡単で、必要なコマンドは3つだけである。このパイプラインは、重複排除ハプロイドアセンブリを生成するためのハプロイドアセンブリ、または、改良された重複排除プライマリハプロイドアセンブリとより完全なセカンダリハプロティグアセンブリを生成するための二倍体アセンブリのいずれかで動作する。最後に、パイプラインはまた、必要に応じてアセンブリの手動検査とキュレーションを支援するように設計されたいくつかの出力を生成する。

 

現在のところ、Purge Haplotigsは二倍体ゲノムに対してのみテストされている。

 

 

f:id:kazumaxneo:20200710173653p:plain

Flow-chart for the Purge Haplotigs pipeline. 論文より転載

 

HaplotigsとPrimary contigs の区別。Pacbioのpdf資料を転載。

f:id:kazumaxneo:20200710195219p:plain

 

チュートリアル

mroachawri / purge_haplotigs / wiki / Tutorial — Bitbucket

 

インストール

ubuntu18.04でcondaの仮想環境を作ってテストした。

 依存

  • Bash
  • BEDTools (tested with v2.26.0)
  • SAMTools (tested with v1.7)
  • Minimap2 (tested with v2.11/v2.12, https://github.com/lh3/minimap2)
  • Perl (with core modules: FindBin, Getopt::Long, Time::Piece, threads, Thread::Semaphore, Thread::Queue, List::Util)
  • Rscript (with ggplot2)

mamba create -n purge_haplotigs_env -y
conda activate purge_haplotigs_env
mamba install -c conda-forge -c bioconda purge_haplotigs -y

#インストール チェック

> purge_haplotigs test

 

[10-07-2020 17:12:19] Writing contig associations

[10-07-2020 17:12:19] Writing the reassignment table and new assembly files

[10-07-2020 17:12:19]

 

PURGE HAPLOTIGS HAS COMPLETED SUCCESSFULLY!

 

samtools faidx curated.fasta

samtools faidx curated.haplotigs.fasta

purge_haplotigs ncbiplace -p curated.fasta -h curated.haplotigs.fasta -t 4 -f

[10-07-2020 17:12:19] minimap2 OK!

[10-07-2020 17:12:19] samtools OK!

[10-07-2020 17:12:19]

Beginning Pipeline

 

PARAMETERS:

Primary contigs FASTA       curated.fasta

Haplotigs FASTA             curated.haplotigs.fasta

Out FASTA                   ncbi_placements.tsv

Threads                     4

Coverage align cutoff       50 %

 

RUNNING USING COMMAND:

purge_haplotigs place -p curated.fasta -h curated.haplotigs.fasta -t 4 -f

 

[10-07-2020 17:12:19] Reading contig lengths

[10-07-2020 17:12:19] Running minimap2 hit search

[10-07-2020 17:12:19] Minimap2 hit search done

[10-07-2020 17:12:19] Reading minimap2 alignments

[10-07-2020 17:12:19] Getting best hits

[10-07-2020 17:12:19] Getting alignments coords for best hits

[10-07-2020 17:12:19] Renaming contigs and haplotigs in the FALCON Unzip format

[10-07-2020 17:12:19] Writing placement file

[10-07-2020 17:12:19] Done!

md5sum -c validate.md5

make: md5sum: No such file or directory

make: *** [Makefile:75: chk_out_files] Error 127

 

ONE OR MORE TESTS FAILED

 

#help 

purge_haplotigs -h

$ purge_haplotigs -h

 

ERROR: no such sub-command -h

 

USAGE:

purge_haplotigs  <command>  [options]

 

COMMANDS:

-- Purge Haplotigs pipeline:

    hist            First step, generate a read-depth histogram for the genome

    cov             Second step, get contig coverage stats and flag 'suspect' contigs

    purge           Third step, identify and reassign haplotigs

 

-- Other scripts

    clip            EXPERIMENTAL; Find and clip overlapping contig ends

    place           Generate a placement file for submission to NCBI

    test            Test everything!

 

purge_haplotigs cov -h

$ purge_haplotigs cov -h

Option high requires an argument

 

USAGE:

purge_haplotigs  cov  -i aligned.bam.genecov  -l <integer>  -m <integer>  -h <integer>  [-o coverage_stats.csv -j 80  -s 80 ]

 

REQUIRED:

-i / -in        The bedtools genomecov output that was produced from 'purge_haplotigs readhist'

-l / -low       The read depth low cutoff (use the histogram to eyeball these cutoffs)

-h / -high      The read depth high cutoff

-m / -mid       The low point between the haploid and diploid peaks

 

OPTIONAL:

-o / -out       Choose an output file name (CSV format, DEFAULT = coverage_stats.csv)

-j / -junk      Auto-assign contig as "j" (junk) if this percentage or greater of the contig is 

                low/high coverage (DEFAULT = 80, > 100 = don't junk anything)

-s / -suspect   Auto-assign contig as "s" (suspected haplotig) if this percentage or less of the

                contig is diploid level of coverage (DEFAULT = 80)

> purge_haplotigs hist

USAGE:

purge_haplotigs  hist  -b aligned.bam  -g genome.fasta  [ -t threads ]

 

REQUIRED:

-b / -bam       BAM file of aligned and sorted reads/subreads to the reference

-g / -genome    Reference FASTA for the BAM file.

 

OPTIONAL:

-t / -threads   Number of worker threads to use, DEFAULT = 4, MINIMUM = 2

-d / -depth     Maximum cutoff for depth. DEFAULT = 200, increase if needed,

                set much higher than your expected average coverage.

> purge_haplotigs purge -h

USAGE:

purge_haplotigs  purge  -g genome.fasta  -c coverage_stats.csv

 

REQUIRED:

-g / -genome        Genome assembly in fasta format. Needs to be indexed with samtools faidx.

-c / -coverage      Contig by contig coverage stats csv file from the previous step.

 

OPTIONAL:

-t / -threads       Number of worker threads to use. DEFAULT = 4

-o / -outprefix     Prefix for the curated assembly. DEFAULT = "curated"

-r / -repeats       BED-format file of repeats to ignore during analysis.

-d / -dotplots      Generate dotplots for manual inspection.

-b / -bam           Samtools-indexed bam file of aligned and sorted reads/subreads to the

                    reference, required for generating dotplots.

-f / -falconNaming  Rename contigs in the style used by FALCON/FALCON-unzip

 

ADVANCED:

-a / -align_cov     Percent cutoff for identifying a contig as a haplotig. DEFAULT = 70

-m / -max_match     Percent cutoff for identifying repetitive contigs. Ignored when 

                    using repeat annotations (-repeats). DEFAULT = 250

-I                  Minimap2 indexing, drop minimisers every N bases, DEFAULT = 4G

-v / -verbose       Print EVERYTHING.

-limit_io           Limit for I/O intensive jobs. DEFAULT = -threads

-wind_min           Min window size for BED coverage plots (for dotplots). DEFAULT = 5000

-wind_nmax          Max windows per contig for BED coverage plots (for dotplots). DEFAULT = 200

 

 

 

実行方法

ランにはドラフトアセンブリ配列が必要。入力ドラフトアセンブリ配列として、ハプロイドアセンブリ(FALCONやCANUなど)、または二倍体アセンブリ(FALCON Unzipなど)のいずれかを使用できる。また、Pacbioのロングリードシークエンシングデータが必要。ショートリードでも動作するが、ロングリードに最適化されている。

 

1、マッピング

PacBioのサブリード、またはショートリードをハプロイドまたは二倍体ゲノムアセンブリマッピングしてbamを作成する。

minimap2 -t 12 -ax map-pb genome.fa subreads.fasta.gz --secondary=no \
| samtools sort -m 1G -@ 12 -o aligned.bam -T tmp.ali

 

2、カバレッジヒストグラムの生成 (スレッド増やしすぎると終わらないので注意

purge_haplotigs hist -b aligned.bam -g genome.fasta -t 12
  • -b    BAM file of aligned and sorted reads/subreads to the reference
  • -g    Reference FASTA for the BAM file
  • -t    Number of worker threads to use, DEFAULT = 4, MINIMUM = 2

BEDToolsのgenomecovファイルと、カバレッジプロファイルを示した画像ファイルが出力される。Collapsed haplotype のコンティグの場合は、両方の対立遺伝子からのリードがマップされるが、対立遺伝子が別々のコンティグとしてアセンブルされた場合は、リードが2つのコンティグに分割され、リードデプスが半分になる。これを利用してハプロチグである可能性の高いコンティグを識別する。

f:id:kazumaxneo:20200711231701p:plain

ハプロイドアセンブリでは、重複が発生している場合は二峰性の分布が観察される。0.5×の深さのピークは重複している領域から、1×の深さのピークはハプロタイプが適切に融合している領域からのものである。二倍体ではアセンブリ全体が複製されている必要があるため、1×のピークは非常に小さいか、または全く見えないかもしれない(論文より)。

 

マニュアルには、phased assemblyである場合は二峰性のピークは二倍体のピークは非常に小さい可能性があるとの記載がある。マニュアルに出力例あり。


3、コンティグごとのカバレッジの分析

コンティグごとにカバレッジを分析する。 前のステップで得られたカバレッジプロファイルからカバレッジカットオフ値を指定する。

purge_haplotigs cov -i aligned.bam.genecov -l <integer> -m <integer> -h <integer> \
-o coverage_stats.csv -j 80 -s 80
  • -i     The bedtools genomecov output that was produced from 'purge_haplotigs readhist'
  • -l      The read depth low cutoff (use the histogram to eyeball these cutoffs)
  • -h     The read depth high cutoff
  • -m    The low point between the haploid and diploid peaks
  • -o     Choose an output file name (CSV format, DEFAULT = coverage_stats.csv)
  • -j      Auto-assign contig as "j" (junk) if this percentage or greater of the contig is  low/high coverage (DEFAULT = 80, > 100 = don't junk anything)
  • -s    Auto-assign contig as "s" (suspected haplotig) if this percentage or less of the contig is diploid level of coverage (DEFAULT = 80)

コンティグのカバレッジ統計のcsvファイルが作成され、疑わしいコンティグ(0.5x)には更なる分析や除去のためのフラグが立てられる (1xのプライマリなコンティグにはフラグは立てられない(論文図2B))。二倍体アセンブリの場合、両方のハプロタイプが存在しているはずなので、ほとんどのコンティグは解析のためにフラグが立てられ、以後の解析でプライマリなコンティグに選抜されると推測される。

数値の例は公式チュートリアルが参考になる。

 

4、パージパイプラインの実行

フラグが立てられたコンティグは、コンパニオンコンティグとのシンテニーの同定を試みるために、Minimap2による配列のアラインメントが行われ、どのコンティグを維持するかが評価される(論文図2C)。アライメントスコアがカットオフ(デフォルトでは70%以上)を超えるコンティグは、ハプロティグとして再割り当ての対象となる。 

purge_haplotigs purge -g genome.fasta -c coverage_stats.csv
  • -g    Genome assembly in fasta format. Needs to be indexed with samtools faidx.
  • -c    Contig by contig coverage stats csv file from the previous step.

ハプロチグが入れ子になっていたり、重複していたり、あるいはほとんどが反復配列で構成されている場合には、コンティグの競合が発生する可能性がある。そのような場合、コンティグの同定、アラインメントスコアリング、コンフリクトの解決、コンティグの再割り当てのステップは条件を満たすコンティグがなくなるまで繰り返し行われる。

 

3つのFASTA形式のファイル;キュレーションされたコンティグ、ハプロティグとして再割り当てされたコンティグ、アーティファクトとして再割り当てされた異常なカバレッジのコンティグ,が出力される。

f:id:kazumaxneo:20200710193908p:plain

  1. <prefix>.fasta.キュレーションされたプライマリコンティグ
  2. <prefix>.haplotigs.fasta. 最初の入力アセンブリで識別されたすべてのハプロチグ
  3. <prefix>.artefacts.fasta. カバレッジが非常に低い/高いコンティグ(注意: ミトコンドリア/葉緑体/その他のコンティグもこの中に存在する可能性が高い)。

 

 

論文にもあるように、結果はBUSCOのduplicated BUSCOなどを見ることで評価できます。

 

追記

purge_haplotigs clip(ベータ)を走らせることで、末端に残ってしまっているオーバーラップをクリッピングできる()。

引用

Purge Haplotigs: allelic contig reassignment for third-gen diploid genome assemblies

Michael J. Roach, Simon A. Schmidt & Anthony R. Borneman
BMC Bioinformatics volume 19, Article number: 460 (2018)

 

参考

https://www.pacb.com/wp-content/uploads/4May2017_JasonChin_ChalengesDiploidAssembly.pdf

 

関連