De novoゲノムアセンブリでは、通常、完全なゲノムではなく、コンティグのセットが作成される。そのため、ゲノムの完全な構造を決定するためには、遺伝的連鎖地図、オプティカルマップ、Hi-Cデータなどの追加データが必要となる。従来の研究では、コンティグを整列・配置するために、このようなデータを利用していた。
本論文では、追加データを用いたゲノムアセンブリのためのフレームワークを紹介する。本アプローチは、リードをクラスタリングし、各クラスタ内のリードが追加データに従ってゲノムの近傍の位置に由来するようにすることに基づく。これらのリードはそれぞれ独立にアセンブルされ、得られたコンティグがさらに階層的にアセンブルされる。HGGAと呼ばれるツールで遺伝的連鎖地図に対する著者らのアプローチを実装した。
Pacific Biosciences社のロングリードと遺伝的連鎖地図のシミュレーションと実測を行った結果、HGGAはリードのアセンブリとコンティグアセンブリを統合した従来の手法に比べて、1.2倍から9.8倍のNGA50またはN50を達成し、より連続性の高いアセンブリと少ないコンティグを生成することが分かった。さらに、リードデータのみを用いたアセンブリと比較して、アセンブリの正しさも同等かより向上している。
我々は以前の研究で、参照配列の代わりに遺伝的連鎖地図を使ってアセンブリを誘導する手法であるKermitを紹介した。遺伝的連鎖地図は、コンティグを染色体内で方向付けして配置し、ミスアセンブリを検出するための手法である。遺伝的連鎖地図自体は、遺伝マーカーで構成されている。マーカーは染色体に分割され、各染色体内では、さらにマーカーはビンに分類される。染色体内のビンの順番は分かっているが、ビン内のマーカーの順番は分かっていない。
地図中のマーカーは、一塩基変異などの変異の集合に由来する。変異は、関連する個体の集団であるシークエンス交配から見出される。ゲノム上で互いに近い変異は、一緒に受け継がれる可能性が高い。したがって、交配個体の遺伝子型を決定し、変異が一緒に遺伝する確率を調べることにより、遺伝的連鎖地図を作成することができる。Kermit はリードセットを遺伝的連鎖地図にマッピングすることで色付けし、色と一致しないエッジをアセンブリグラフから削除する。同じ方法は、理論的にはオプティカルマップのようにリードの線形順序として表現できるあらゆるガイドデータに拡張できるが、非線形ガイドデータへの一般化には失敗している。
Kermitでは、アセンブリグラフの構築とグラフを用いたゲノムアセンブリの両方にminiasmを使用している。ここでは、比較しやすいように、同様にminiasmを使用する。miniasmは、まずminimapを用いてリードの重なりを見つける。そして、その重複部分を元に、交差する辺を削除して文字列グラフを作成する。次に、グラフの先端を削除したり、バブルをはじいたりして、グラフをきれいにする。最後に、miniasmは、得られたグラフのunitigs（分岐していないパス）をcontigsとして報告する。Miniasmはコンセンサスフェーズを実装していないため、miniasmが生成するcontigsのエラー率はreadのエラー率と同じになる。
インストール
依存
- python3
- pysam
#本体
wget -O NC_000913.fa 'http://www.ncbi.nlm.nih.gov/sviewer/viewer.cgi?sendto=on&db=nuccore&dopt=fasta&val=556503834'
git clone https://github.com/rikuu/hgga && cd hgga
#依存
#./minimap2/minimap2のようにminimap2やminiasmを参照するので、hggaのホームディレクトリにcloneしてビルドする。
git clone https://github.com/lh3/miniasm && (cd miniasm && make)
git clone https://github.com/lh3/minimap2 && (cd minimap2 && make)
git clone https://github.com/lbcb-sci/racon && cd racon && mkdir build && cd build
cmake -DCMAKE_BUILD_TYPE=Release .. && make -j8 && cd ../..
#pysamが必要
conda install -c bioconda pysam
> python hgga.py
Usage: hgga.py <reads> <alignments> [number of threads] [min leaf size]
> python closest_color.py
Usage: closest_color.py <mappings> <markers> [markers...]
テストラン
１、ゲノム、ロングリード、linkage mapをダウンロードする。
#genome
wget -O NC_000913.fa 'http://www.ncbi.nlm.nih.gov/sviewer/viewer.cgi?sendto=on&db=nuccore&dopt=fasta&val=556503834'
#PacBio reads
wget -O- http://www.cbcb.umd.edu/software/PBcR/data/selfSampleData.tar.gz | tar zxf -
ln -s selfSampleData/pacbio_filtered.fastq reads.fq
#linkage map
wget -O ecoli.txt https://raw.githubusercontent.com/rikuu/kermit/master/misc/map.txt
連鎖地図は、コンティグ上の位置とbinでマーカーを定義したタブ区切りのファイルとして提供する。大腸菌K-12の連鎖地図のシミュレーションがkermit/misc/に用意されている。
２、ロングリードをリファレンスにアラインしてPAF形式で出力、coloringする。
minimap2/minimap2 -x map-pb NC_000913.fa reads.fq > align.paf
python closest_color.py align.paf ecoli.txt > color.cf
（マップが十分に密になると、アセンブリの品質はほとんど変化しないことが分かる。Kermit では密度が 50k に達したとき、HGGA では密度が 10k に達したときに、このような現象が起こる（論文より））
ガイドアセンブルする。ロングリードとcolor.cf（PAF形式ではない）を指定する。
python hgga.py reads.fq color.cf
引用
HGGA: hierarchical guided genome assembler
Riku Walve & Leena Salmela
BMC Bioinformatics volume 23, Article number: 167 (2022)
関連
ホタペンさんのリツートでこの論文に気づきました。ありがとうございました。
