macでインフォマティクス

macでインフォマティクス

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

Pacbioシーケンシングリードのオーバーラップ検出感度を改善する GroupK

 

 リード長の増加により、第3世代のシークエンシングでゲノムアセンブリのギャップを埋め[ref.1, 2]、構造の変化を明らかにし[ef.13]、トランスクリプトームシークエンシングで遺伝子アイソフォームをより正確に定量できるようになった[ef.14]。さらに、ロングリードを使用することは、微生物群集構造を明らかにし、微生物群集の種内変動を解読することにおいて有望である[ref.5、6]。

 第三世代のシークエンシングデータを用いたゲノムの構築には、専用の方法とツールが必要になる。既存のゲノム構築ツールは、主に2つのタイプのグラフモデル、すなわちオーバーラップグラフおよびde Brujinグラフを利用する。エラー率が低い場合、de Bruijnグラフには、シーケンスサイズが大きくなってもグラフサイズが大きく増加しないという理論的な利点がある。これは通常、イルミナのデータセットでは高い値である。第3世代のシークエンシングデータでは、高いエラー率と低いカバレッジにより、オーバーラップグラフがゲノム構築に適した選択肢となる[ref.7]。オーバーラップグラフを構築する際の重要なステップは、オーバーラップを共有するリードペアを同定することであり、これはそれらリードが基礎となるゲノム内の同じ遺伝子座からシーケンシングされたことを示す。オーバーラップアラインメントを実施するために利用可能な多数のシーケンスアラインメントプログラムがあるが[ref.8、9]、それらの大部分は動的計画法に頼っており、そしてハイスループットシーケンシングデータにとって計算コストがかかりすぎる。エラー率が高いため、BWT(Burrows-Wheeler変換)またはハッシュテーブル[ref.10、11]を使用した既存のショートリードのオーバーラップ検出ソフトウェアは、ロングリードに直接適用することはできない。

 エラーが発生しやすいロングリードの重複を検出するために、現在2つの方法が採用されている。1つの戦略は、オーバーラップ検出の前に、Pacbio(パシフィックバイオサイエンス)およびONT(オックスフォードナノポア)データにおけるシーケンシングエラーを修正することを試みるものである。様々なシーケンシングエラー訂正ツールが存在する[ref.7、12]。それらのうちのいくつかは、少なくとも2つのシーケンシングライブラリーおよび数種類のシーケンシングラン調製を必要とし、したがって多くの用途にとって費用効果が高くないハイブリッドシーケンシングに依存している。他のものはロングリードだけを使用してエラー訂正を行う。代表的な方法の1つが、ChinらのHGAP [ref.12]に記載されており、その性能は、より高いリードカバレッジ範囲で向上する。 HGAPのようなアライメントベースのエラー訂正方法にとって、重要なステップは迅速にアライメントできるリードを識別することである。本質的に、オーバーラップ検出に使用される技術はアライメント検出にも同様に使用することができる。

 2番目のカテゴリは、エラー訂正の難しさを回避し、rawリードを使用して重複を識別する。 PacBioおよびONTデータには、さまざまな近似類似検索方法が適用されている[ref.13]。これらは一般的にシードチェーンアラインメントの手順に従う[ref.14]。シードベースのフィルタリングステップは、感度と計算効率の間のトレードオフを制御する上で重要な役割を果たす。通常、これらの方法では、フィルタリングステップとして短い文字列の(完全)一致を使用する。短い文字列またはk-merの一致には、2つのシーケンス間のk個の連続した文字の完全一致が必要である。直感的には、オーバーラップするリードは、オーバーラップしないリードよりも一般的なk-merを共有する傾向がある。したがって、共有k-mer数を迅速に見つけることができる戦略を適用することができる。このセクションでは、いくつかの最先端のオーバーラップ検出ツールの主な戦略をまとめる。次のセクションで、我々(著者ら)の方法と既存の方法との違いを強調する。

(3段落省略)

 この研究では、ゲノム構築にあたりオーバーラップしたロングリードを検出するため、グループシードを使用する。 GroupKと呼ばれるこの実装は、小さなオーバーラップまたは両方のリードの高いエラー率によって悪化したオーバーラップを検出するための既存の方法を補完するツールを提供する。この機能により、GroupKは、第3世代のシーケンシングプラットフォームによってシーケンスされたメタゲノムデータにおけるゲノム構築のための賢明な選択を可能にする。これらの集団サンプルは通常不均一な範囲の微生物を含んでいるので、小さなオーバーラップを識別することができることは希少な種のゲノムを再構築するために非常に重要になるだろう。

 

 

インストール

GroupKのこのpythonコードは実験のための古いコードになっている( C ++の実装は今後リリース予定と記載されている)。

  • Python 3 (You may need to modify some script if you want to use Python 2.7)
  • CMake 3.7 or higher
  • GCC 4.8 or higher
  • Biopython library for Python
  • SortedContainers library fot Python
conda install -y sortedcontainers
conda install -y matplotlib
conda install -c bioconda -y sortedcontainers

本体 Github

git clone https://github.com/Strideradu/GroupK.git 
cd GroupK/sa_filter/
g++ -std=c++11 *.cpp -o filter

python GroupK.py -h

# python GroupK.py -h

usage: GroupK.py [-h] [--k1 K1] [--threshold THRESHOLD] [--k2 K2]

                 [--accuracy ACCURACY] [--gap GAP] [--chain CHAIN]

                 [--groupc GROUPC] [--idbase IDBASE] [--ratio RATIO]

                 [--size SIZE] [--large LARGE]

                 input output

 

positional arguments:

  input                 path of input fasta file

  output                path of output file

 

optional arguments:

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

  --k1 K1               kmer size for filtration

  --threshold THRESHOLD

                        count threshold for shared k1

  --k2 K2               kmer size for group

  --accuracy ACCURACY   accuracy of the reads

  --gap GAP             gap rate

  --chain CHAIN         number of kmer in the chain to report

  --groupc GROUPC       the coefficeint c in paper to control chaining

                        threshold, must be non-zero float

  --idbase IDBASE       if the number of identical base meet this threshold

                        the output will always be reported

  --ratio RATIO         ratio of two overlap region threshold

  --size SIZE           group size threshold

  --large LARGE         release threhold for large overlap

usage: GroupK.py [-h] [--k1 K1] [--threshold THRESHOLD] [--k2 K2]

                 [--accuracy ACCURACY] [--gap GAP] [--chain CHAIN]

                 [--groupc GROUPC] [--idbase IDBASE] [--ratio RATIO]

                 [--size SIZE] [--large LARGE]

                 input output

 

positional arguments:

  input                 path of input fasta file

  output                path of output file

 

optional arguments:

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

  --k1 K1               kmer size for filtration

  --threshold THRESHOLD

                        count threshold for shared k1

  --k2 K2               kmer size for group

  --accuracy ACCURACY   accuracy of the reads

  --gap GAP             gap rate

  --chain CHAIN         number of kmer in the chain to report

  --groupc GROUPC       the coefficeint c in paper to control chaining

                        threshold, must be non-zero float

  --idbase IDBASE       if the number of identical base meet this threshold

                        the output will always be reported

  --ratio RATIO         ratio of two overlap region threshold

  --size SIZE           group size threshold

  --large LARGE         release threhold for large overlap

 

 

 

実行方法

FASTA形式のPacbioのシーケンシングリードを指定する。

python GroupK.py --k1 15 --threshold 2 --k2 9\
input.fasta output
  • --k1    kmer size for filtration (default: 15)
  • --threshold    count threshold for shared k1 (default: 2)
  • --k2    kmer size for group (default: 9)
  • --accuracy    accuracy of the reads (default: 0.85)

テスト時はエラーになった。

  

引用

Improving the sensitivity of long read overlap detection using grouped short k-mer matches

Nan Du, Jiao Chen, Yanni Sun
BMC Genomics 2019 20 (Suppl 2) :190