macでインフォマティクス

macでインフォマティクス

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

オルソログとパラログを小メモリ使用量で高速探索する SwiftOrtho

 2019 10/26 論文引用追加

 

 Gene homology type classification は、種を越えてパラログおよびオルソログを同定することからなる。オルソログは共通の先祖遺伝子から分化後に進化した遺伝子であり、パラログはduplicationのために相同な遺伝子である。遺伝子の進化の歴史は遺伝子機能と進化の理解に影響を与えるため、種を越えてオルソログとパラログを計算によって検出することは重要な問題である。
 Homology type の適切な推論は系統樹を用いた遺伝子履歴の追跡を含む[ref.1]が、いくつかの代用方法が長年にわたって開発されてきた。オルソログを推論する最も一般的な方法は、Reciprocal Best Hit (相互ベストヒット)またはRBHである[ref.2、3]。簡単に言うと、RBHは次のように述べられている。異なる2ゲノムの2つの遺伝子によってコードされている2つのタンパク質が互いに最良のスコアリングマッチとして見つかると、それらはオルソログと見なされる[ref.2、3]。

 Inparanoid(wiki)は、オルソログとパラログ内の両方を含むようにRBHオルソロジー関係を拡張する[ref.4–6]。具体的には、Inparanoidは、オルソログとIn-パラログ(link)を識別する。これらは、特定の分化イベント後に複製されたものである[ref.4–6]。その場合、当然のことながら、2種間のオルソログペアをオルソロググループに拡張することができる。ここでオルソロググループは共通の祖先から派生したと仮定される遺伝子のセットとして定義される[ref.6]。複数の種にわたるオルソログ群を同定するためにいくつかの方法が開発されてきた。これらの方法は、ツリーベースとグラフベースの2種類に分類できる。ツリーベースの方法は、異なる種における相同配列のアラインメントから遺伝子ツリーを構築し、その遺伝子ツリーをその対応する種ツリーと一致させることによってオルソロジー関係を推測する[ref.1、7、8]。正しい遺伝子ツリーと種のツリーが与えられていれば、ツリーベースの方法は正しいオルソロジー関係を推論できる[ref.9]。ツリーベースの方法の主な制限は、与えられた遺伝子ツリーと種ツリーの精度である。誤ったツリーは、不適切なオルソログとパラログの割り当てにつながる[ref.8–10]。ツリーベースの方法も計算コストが高く、それらを多数の種に適用する能力を制限している[ref.9、11–13]。グラフベースの方法は、ホモログからオルソログとパラログ内(図1)を推論し、次にそれらをオルソログ群に分類するために異なる戦略を使用する[ref.8、11、12]。オルソロググループのクラスタまたはCOGデータベースのクラスタは、3つの異なる種のRBHの三角形を検出し、三角形の共通の辺でマージする[ref.14]。Orthologous Matrix(OMA)は、類似度グラフから最大重みクリークを見つけることによって、RBHをオロソロググループに分類する[ref.15]。 Multi-ParanoidはInparanoidを拡張したもので、InParanoidを使用して3つの異なる種の三角オルソログとIn-パラログをシードとして検出し、それらのシードをより大きなグループにマージする[ref.16]。 OrthoMCLはまた、InParanoidアルゴリズムを使用して2つの種間のオルソログ、co-オルソログ(意味)、およびIn-パラログを検出し[ref.17]、次にマルコフクラスタリング(MCL)[ref.18]を使用してこれらの関係をオルソロググループに分類する。前者は遺伝子および種の進化の歴史から相同性タイプを直接推論するのではなく、代用法を用いてオルソログおよびパラログ内を同定するため、理論的には、グラフベースの方法はツリーベースの方法よりも正確性が低い。実際には、グラフベースの方法は、ツリーベースの方法と同様の精度を持つ[ref,9、10、19]。ツリーベースの方法とグラフベースの方法の両方を含むいくつかの方法を比較すると、ツリーベースの方法の方が大きなデータセットに対するグラフベースの方法よりもパフォーマンスがさらに悪いことがわかっった[ref.10]。ある研究では、RBH、グラフベースおよびツリーベースを含むいくつかの一般的な方法が比較され、ツリーベースの方法ではより高い特異性が得られるが感度が低いことが多くある[ref.20]。いくつかの研究はまた、グラフベースの方法が、ツリーベースの方法よりも特異性と感度の間でより良いトレードオフを見つけることを示している[ref.10、20、21]。それらのより良い速度および正確さのために、グラフに基づく方法は一般に大きなデータセットを分析するために好ましい。 OrthoMCLやInParanoidなどのグラフベースのメソッドは、何百ものゲノムを解析できるが、それらはすぐには利用できない可能性があるかなりの計算リソースを必要とする[ref.22、23]。

 ここではSwiftOrthoという名前の新しいオルソログ解析ツールを開発した。 SwiftOrthoは、速度、精度、メモリ効率に焦点を当てたグラフベースの方法である。著者らは、SwiftOrthoを、ゴールドスタンダードデータセットOrthobench [ref.12]とQuest for Orthologsサービス[ref.24]を使って、いくつかの既存のグラフベースのツールと比較した。両方のベンチマークを使用して、SwiftOrthoが他のグラフベースの方法よりも低いCPUおよびメモリ使用量で高い精度を提供することを示す。

 

SwiftOrthoに関するツイート

 

 

インストール

ubuntu18.04、python2.7.14環境でテストした。

依存

  • Python2.7 (Recommend Anaconda ) or PyPy2.7(v5.10 or greater) and Packages:
  • Networkx
  • RPython
  • numpy
  • scipy
  • Biopython
  • sklearn
  • MCL(optional)

Install packages via pip:

pip install -U rpython networkx scipy numpy biopython sklearn

本体 Github

git clone https://github.com/Rinoahu/SwiftOrtho.git
cd SwiftOrtho/

> python bin/find_hit.py -h

$ python bin/find_hit.py -h

Usage:

  search:

    python bin/find_hit.py -p blastp -i qry.fsa -d db.fsa

Parameters:

  -p: program

  -i: query sequences in fasta format

  -l: start index of query sequences

  -u: end index of query sequences

  -L: start index of reference

  -U: end index of reference

  -d: ref database

  -D: index of ref, if this parameter is specified, only this part of formatted ref will be searched against

  -o: output file

  -O: write mode of output file. w: overwrite, a: append

  -s: spaced seed in comma separated format: 1111,1110,1001

  -r: reduced amino acid alphabet. Available value is:

      1. aa9, the default setting which stands for AST,CFILMVY,DN,EQ,G,H,KR,P,W

      2. aa20, the fast setting which stands for the original alphabet of 20 amino acids

      3. user defined alphabet. The 20 amino acids in comma separated format,for example: AST,CFIL,MVY,DN,EQ,G,H,KR,P,W

  -v: number of hits to show

  -e: expect value

  -m: max ratio of pseudo hits that will trigger stop

  -j: distance between start sites of two neighbor seeds, greater will reduce the size of database

  -t: filter high frequency kmers whose counts > t

  -F: filter query sequence

  -M: bucket size of hash table, reduce this parameter will reduce memory usage but decrease sensitivity

  -c: chunck size of reference. default is 100K which means 100K sequences from reference will be used as database

  -a: number of processors to use

  -T: tmpdir to store tmp file. default ./tmpdir

 >python scripts/run_all.py 

Integrate pipeline can be used to identify orthologs, construct phylotree and get profile of pan-genomes

If the operonic information supplied, this pipeline can also cluster the operon_cluster

Usage:

  python run_all.py -i foo.pep.fsa [-r taxonomy] [-p foo.operon] [...]

Input fasta file:

 -i: protein sequences in fasta formt. The header should be like xxxx|yyyy: xxxx is taxon name and yyyy is unqiue identifier in that taxon

 

Operon information if operon clustering is needed:

 -p: operonic annotation file. The 1st column of this file should be like x0-->x1-->x2-->x3 or x0<--x1<--x2<--x3.

 

Optional parameters for all-vs-all homologous search:

 -s: spaced seed in format: 1111,11110,1001111 parameter

 -a: number of processors to use

 

Optional parameters for orthology inference:

 -c: min coverage of sequence [0~1]

 -y: identity [0~100]

 -n: normalization score [no|bsr|bal]. bsr: bit sore ratio; bal:  bit score over anchored length. Default: no

 

Optional parameters for clustering:

 -A: clustering algorithm. [mcl|apc]

 

Optional parameters for species tree construction:

 -r: taxonomy name used as reference [optional]

 

Optional parameters for pan-genome:

 -l: threshold to identify specific genes. For example, 0.05 means the specific gene should only exist in less than 5% of sll selected pecies

 -u: threshold to identify core genes. For example, 0.95 means core gene should exist in atleast 95% of all selected species 

 -I: Inflation parameter of mcl. default: 1.5

 

 

テストラン

1、All-to-all の相同性検索。example/ref.fsaには4生物のプロテオーム(protein.fasta)が連結して収納されている。

python bin/find_hit.py -p blastp -i example/qry.fsa -d example/ref.fsa -o input.fsa.sc -e 1e-5 -s 111111
  • -i     query sequences in fasta format
  • -d    ref database
  • -o    output file
  • -e    expect value
  • -s    spaced seed in comma separated format: 1111,1110,1001

> cat input.fsa.sc

GCF_000005825.2_ASM582v2|BPOF4_RS00005 GCF_000005825.2_ASM582v2|BPOF4_RS00005 100.00 450 0 0 1 450 1 450 2.88e-261 897 450 450 0 GCF_000005825.2_ASM582v2|BPOF4_RS00005 fake news

GCF_000005825.2_ASM582v2|BPOF4_RS00005 GCF_000006605.1_ASM660v1|JK_RS00005 53.52 340 158 0 111 450 240 579 1.60e-105 380 450 583 0 GCF_000006605.1_ASM660v1|JK_RS00005

GCF_000005825.2_ASM582v2|BPOF4_RS00005 GCF_000006645.1_ASM664v1|Y_RS21000 39.91 466 280 19 1 450 14 462 2.66e-99 359 450 462 0 GCF_000006645.1_ASM664v1|Y_RS21000

GCF_000005825.2_ASM582v2|BPOF4_RS00005 GCF_000005845.2_ASM584v2|b3702 48.20 363 188 4 94 450 106 467 3.44e-97 352 450 467 0 GCF_000005845.2_ASM584v2|b3702

GCF_000005825.2_ASM582v2|BPOF4_RS00005 GCF_000006625.1_ASM662v1|UU_RS00020 24.72 457 344 17 2 442 3 449 4.43e-33 139 450 457 0 GCF_000006625.1_ASM662v1|UU_RS00020

GCF_000005825.2_ASM582v2|BPOF4_RS00005 GCF_000005845.2_ASM584v2|b2496 25.00 232 174 10 116 346 15 232 5.34e-12 68 450 233 0 GCF_000005845.2_ASM584v2|b2496

GCF_000005825.2_ASM582v2|BPOF4_RS00005 GCF_000006645.1_ASM664v1|Y_RS07375 25.72 206 153 10 114 318 15 206 1.07e-11 67 450 235 0 GCF_000006645.1_ASM664v1|Y_RS07375

 

 2、オロソログ検索。1の出力を使う。

python bin/find_orth.py -i input.fsa.sc -c 0.5 -y 0 > input.fsa.sc.orth

 

3、クラスタリング

python bin/find_cluster.py -i input.fsa.sc.orth -a mcl -I 1.5 > input.fsa.sc.orth.mcl

 

エラーになる。修正できたら追記します。 

引用

SwiftOrtho: a Fast, Memory-Efficient, Multiple Genome Orthology Classifier

Xiao Hu, Iddo Friedberg

bioRxiv preprint first posted online Feb. 8, 2019

 

追記

SwiftOrtho: A fast, memory-efficient, multiple genome orthology classifier
Xiao Hu, Iddo Friedberg Author Notes
GigaScience, Volume 8, Issue 10, October 2019

 

関連