macでインフォマティクス

macでインフォマティクス

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

ロングリードのself error correctionやcontigのポリッシングを行う CONSENT

2019 4/16 マッピングの画像追加

2019 7/22 インストール、help追記、エラー修正

2019 9/8 コメント追加

2019 11/11 Segmentation faultのリンク追記

  

 第3世代のシークエンシング技術Pacific BiosciencesとOxford Nanoporeは、2011年の創業以来広く使用されてきた。 このロングリードは、コンティグおよびハプロタイプアセンブリ(Patterson et al、2015; Kamath et al、2017)、スキャフォールド(Cao et al、2017)、および構造変異のコール(Sedlazeck et al、1995)などの様々な問題を解決すると期待される。しかし、それらは非常にノイズが多く、10〜30%のエラー率に達する。ショートリードのエラーはほとんど置換で構成されているが、ロングリードのエラープロファイルは主に挿入と欠失で構成されているため、ショートリードよりもはるかに複雑である。その結果、ロングリードを扱うプロジェクトの最初のステップとして、エラー訂正がしばしば必要になる。ロングリードのエラープロファイルとエラーレートはショートリードのものとはかなり異なるため、ロングリードを訂正するには特定のアルゴリズム開発が必要になる。
 ロングリードのエラー訂正は、主に2つの方法で対処されてきた。第1のアプローチ、ハイブリッドエラー訂正は、訂正を実行するために追加のショートリードデータを利用する。 2番目のアプローチである自己訂正は、ロングリードをそのシーケンスに含まれる情報のみに基づいて修正することを目的としている。ハイブリッドの方法は、次のようなさまざまな技術に依存している。
•ショートリードとロングリードのアライメント(CoLoRMAP(Haghshenas et al、2016)(紹介)、HECiL(Choudhury et al、2018)(紹介))。
•de Bruijnグラフ探索(LoRDEC(Salmela and Rivals、2014)(紹介)、Jabba(Miclotte et al、2016)(紹介)、FMLRC(Wang et al、2018)(紹介))。
•ショートリードから構築されたコンティグへのロングリードのアライメント(MiRCA(Kchouk and Elloumi、2018)、HALC(Bao and Lan、2017)(紹介))。
隠れマルコフモデル(Hercules(Firtina et al、2018)(紹介))
•異なる戦略の組み合わせ(NaS(Madoui et al、2015)、HG-
CoLoR(Morisse et al、2018)(紹介))
自己訂正法は通常、ロングリードの互いのアライメントに基づいて構築されている(PBDAGCon(Chin et al、2013)(紹介)、PBcR(Koren et al、2013))。
最初のロングリードシーケンシング実験の結果、非常にエラーの多いロングリード(平均15〜30%のエラー率)が得られたため、ほとんどの方法はショートリードデータの追加使用に頼っていた。その結果、ハイブリッド訂正は以前よりもずっと研究され、そしてより高度に開発された。実際、2014年には、5つのハイブリッド訂正ツールと2つの自己訂正ツールしか使用できなかった。ただし、第3世代のシーケンシング技術は急速に進化しており、現在では10〜12%のエラー率に達するロングリードを生み出すことに成功している。さらに、ロングリードシーケンシングテクノロジの進化により、低コストでより高いスループットのデータを生成することも可能になり、その結果、そのようなデータはより広く利用可能になった。したがって、現在、自己訂正はロングリードを扱うデータ分析プロジェクトの最初のステップとして頻繁に使用されるようになった。

 

関連
 第三世代のシーケンシング技術の急速な発展、およびそれらが現在到達しているより低いエラー率のために、様々な効率的な自己訂正法が開発された。それらのほとんどは、ロングリード間のオーバーラップというコンピューティングの最初のステップを共有している。これはマッピングアプローチを介して行うことができ、これはロングリードの類似領域の位置のみを提供する(Canu(Koren et al、2017)、MECAT(Xiao et al、2017)、FLAS(Bao et al、2018))、またはアラインメントにより、類似領域の位置、およびそれらの実際の塩基対塩基の一致、ミスマッチ、挿入、および欠失に関する一致を提供する(PBDAGCon、PBcR、Daccord(Tischler and Myers、2017))。次いで、必要ならば、マッピングされた領域の実際のアラインメントを再計算した後に、アラインメントを要約し、コンセンサスを計算するために、directed acyclic graph(DAG)が通常構築される。

(一部略)
最先端技術とは異なる効率的なアプローチを組み合わせた新しいself error correction方法であるCONSENTを紹介する。コンセンサス配列を計算するために、一致は実際にはロングリードの重複領域間の複数のシーケンスアラインメントを計算することから始まる。残りのエラーを修正し、最終的なエラー率を下げるために、これらのコンセンサスシーケンスはローカルのde Bruijn graphの助けを借りてさらに洗練される。さらに、他の現在の最先端の方法とは異なり、CONSENTは partial order graphsに基づく方法を使用して実際のマルティプルシーケンスアラインメントを計算する(Lee et al、2002)。加えて、 k-mers chainingに基づく効率的なセグメンテーション戦略を導入し、それは複数の配列シーケンスアラインメントの時間フットプリントを減らす。したがって、このセグメンテーション戦略は、スケーラブルな複数のシーケンスアライメント計算を可能にする。特に、CONSENTをオックスフォードナノポアのウルトラロングリードに効率的に拡張することができる。
 本著者らの実験は、CONSENTが最新の最先端のself error correction法とうまく比較でき、そしてOxford Nanoporeのリアルデータセットでそれらを凌駕することさえ示している。特に、CONSENTはオックスフォードナノポアのウルトラロングリードを含むヒトのシーケンシングデータセットに拡張できる唯一の方法であり、最大340 kbpの長さに達することを示している。

 

f:id:kazumaxneo:20190327204618p:plain

Fig. 1: Overview of CONSENT’s worfklow for long read error correction. 論文より転載

 

CONSENTに関するツイート


インストール

ubuntu16.04のPython 3.6.7環境でテストした(docker使用、ホストOS macos 10.14)。

依存

CONSENT has been developed and tested on x86-64 GNU/Linux.

#fpaの導入
#1 bioconda (link) *注1
conda install -c bioconda -y fpa

#Rustのcargoを使う
cargo install fpa_lr

本体 Github

git clone --recursive https://github.com/morispi/CONSENT
cd CONSENT/
./install.sh
#注*2

> ./CONSENT-correct

# ./CONSENT-correct 

Usage: ./CONSENT-correct [options] --in longReads.fasta --out result.fasta --type readsTechnology

 

Input:

longReads.fasta:               fasta file of long reads to correct, with one sequence per line.

result.fasta:                  fasta file where to output the corrected long reads.

readsTechnology:               Indicate whether the long reads are from PacBio (--type PB) or Oxford Nanopore (--type ONT)

 

Options:

--windowSize INT, -l INT:      Size of the windows to process. (default: 500)

--minSupport INT, -s INT:      Minimum support to consider a window for correction. (default: 4)

--maxSupport INT, -S INT:      Maximum number of sequences to include into a window. (default: 1,000)

--maxMSA INT, -M:              Maximum number of sequences to include into the MSA. (default: 150)

--merSize INT, -k INT:         k-mer size for chaining and polishing. (default: 9)

--solid INT, -f INT:           Minimum number of occurrences to consider a k-mer as solid during polishing. (default: 4)

--anchorSupport INT, -c INT:   Minimum number of sequences supporting (Ai) - (Ai+1) to keep the two anchors in the chaining. (default: 8)

--minAnchors INT, -a INT:      Minimum number of anchors in a window to allow consensus computation. (default: 10)

--windowOverlap INT, -o INT:   Overlap size between consecutive windows. (default: 50)

--nproc INT, -j INT:           Number of processes to run in parallel (default: number of cores).

  --minimapIndex INT, -m INT:    Split minimap2 index every INT input bases (default: 500M).

--tmpdir STRING, -t STRING:    Path where to store the temporary overlaps file (default: working directory, as Alignments_dateTimeStamp.paf).

--help, -h:                    Print this help message.

> ./CONSENT-polish

# ./CONSENT-polish

Usage: ./CONSENT-polish [options] --contigs contigs.fasta --reads longReads.fasta --out result.fasta

 

Input:

contigs.fasta:                 fasta file of contigs to polish, with one sequence per line.

longReads.fasta:               fasta file of long reads to use for polishing, with one sequence per line.

result.fasta:                  fasta file where to output the polished contigs.

 

Options:

--windowSize INT, -l INT:      Size of the windows to process. (default: 500)

--minSupport INT, -s INT:      Minimum support to consider a window for correction. (default: 4)

--maxSupport INT, -S INT:      Maximum number of sequences to include into a window. (default: 1,000)

--maxMSA INT, -M:              Maximum number of sequences to include into the MSA. (default: 150)

--merSize INT, -k INT:         k-mer size for chaining and polishing. (default: 9)

--solid INT, -f INT:           Minimum number of occurrences to consider a k-mer as solid during polishing. (default: 4)

--anchorSupport INT, -c INT:   Minimum number of sequences supporting (Ai) - (Ai+1) to keep the two anchors in the chaining. (default: 8)

--minAnchors INT, -a INT:      Minimum number of anchors in a window to allow consensus computation. (default: 10)

--windowOverlap INT, -o INT:   Overlap size between consecutive windows. (default: 50)

--nproc INT, -j INT:           Number of processes to run in parallel (default: number of cores).

  --minimapIndex INT, -m INT:    Split minimap2 index every INT input bases (default: 500M).

--tmpdir STRING, -t STRING:    Path where to store the temporary overlaps file (default: working directory, as Alignments_dateTimeStamp.paf).

--help, -h:                    Print this help message.

> bin/CONSENT

# ./CONSENT    

Usage: ./CONSENT [-a alignmentFile.paf] [-s minSupportForGoodRegions] [-l minLengthForGoodRegions] [-j threadsNb] 

 

 

 

実行方法

Oxford nanoporeのロングリードをエラー訂正する。 ロングリードのfasta or fastqを指定する(*注3)。

CONSENT-correct --in longReads.fasta --out CONSENT_result.fasta --type ONT
  • readsTechnology: Indicate whether the long reads are from PacBio (--type PB) or Oxford Nanopore (--type ONT)

厳しくするなら--anchorSupportの値を増やすのが効果的。

ディープシーケンスなら--maxSupportや--maxMSAを増やす。同時に--anchorSupportをかなり増やす。

 

contigのpolishingを行う。contigのfastaと、アセンブルに使用したロングリードのfastaを指定する。

CONSENT-polish --contigs contigs.fasta --reads longReads.fasta --out result.fasta

 

CONSENT-correctでエラー訂正してマッピングしてみた。見事にエラーが消えている。

f:id:kazumaxneo:20190416222346j:plain

上がエラー訂正前、下がエラー訂正後のリファレンスへのマッピング。不完全なindel部位はホモのindelに変化している。

 

CONSENTはマッピングのデプスやエラー率のパラメータ設定によって、defaultの条件では全然修正できない場合があります。特定の領域のリードだけ取り出し、小さなデータでパラメータ検討することをお勧めします。

引用

CONSENT: Scalable self-correction of long reads with multiple sequence alignment
Pierre Morisse, Camille Marchet, Antoine Limasset, Thierry Lecroq, Arnaud Lefebvre

bioRxiv preprint first posted online Feb. 11, 2019

 

関連 

self error correction


 

 

パフォーマンスを比較したプレプリントが出ています。どのツール/アルゴリズムを使うべきか考える指針になります。

https://www.biorxiv.org/content/biorxiv/early/2019/05/29/519330.full.pdf 

 

*1

2019/7/22に別のマシンに導入する際、fpaのversion0.5.0を入れるとCONSENTのランの途中にfpaのエラがー起きた。その時はpyenvでcondaのバージョン管理しているサンプルだったので、それも関係している可能性があるが、*2の流れで応急処置した。

 

*2

/usr/local/bin/binにコマンドがないとしてエラーが出た。スクリプト自体は修正せず、/usr/local/bin/binと/usr/local/bin/BMEAN/BOAをそれぞれ作成し、そこにバイナリとスクリプトを入れて対応した。

mkdir /usr/local/bin/bin
mkdir /usr/local/bin/BMEAN && /usr/local/bin/BMEAN/BOA

cd CONSENT/
cp CONSENT-* /usr/local/bin/bin/
cp bin/CONSENT /usr/local/bin/bin/
cp bin/reformatPAF.py /usr/local/bin/bin/
cp BMEAN/blosum80.mat /usr/local/bin/BMEAN/BOA/

 

*3

後々試した際(*1のエラーと同じ日付)、fastqを使うとエラがー出たのでfastaに変換してから使用した。

 

追記

ヘッダーにスペースがあるとSegmentation faultが起きる

https://github.com/morispi/CONSENT/issues/4