マイクロサテライトを高速検索する PERF

　反復DNAは複雑な生物中のゲノムのかなりの部分を構成し、i） interspersed repeats （以下、散在性反復配列）または transposable elements （以下、転移因子）とii）タンデムリピートの2つのカテゴリーに大別できる（Kumar et al、2010）。反復モチーフの長さに依存して、タンデムリピートは、サテライト（> 100nt）、ミニサテライト（10〜30nt）およびマイクロサテライト（1-6nt）に分類される。シンプルシーケンスリピート（SSR）またはショートタンデムリピート（STR）としても知られているマイクロサテライトは、ゲノムに無作為に分布しており、ポリメラーゼスリップに起因して高い突然変異率を示し、長さの増加または減少を引き起こす可能性がある（Ellegren 、2004）。その多型性が高いため、SSRはリンケージ解析（Hearne et al、1992）、ジェノタイピング（Kashi et al、1997）およびDNAフィンガープリンティング（Zietkiewicz et al、1994）に非常に有用である。コード領域におけるSSRの長さの変化は、ハンチントン病および脊髄小脳失調症（Usdin、2008）などのヒトのいくつかの神経変性疾患に関連する。マイクロサテライトは、遺伝子発現のエピジェネティックな調節（Greene et al、2007; Pietrobono et al、2005）およびゲノム構成（Kumar et al、2013; Pathak et al、2013）において重要な役割を果たすことも示されている。

　それらの有用性を考慮すると、SSRの効率的な同定は、計算生物学における長年の目標であった。 SSRを特定するにはいくつかのツールがあるが、速度、効率、包括性、精度、使いやすさ、柔軟性という点で多くの注意点がある。既存の方法は、ヒューリスティック手法またはコンビナトリアル手法のいずれかを使用して、DNA配列からSSRを見つける（Lim et al、2013）。 TRFは、反復（Benson、1999）を同定するためにBernoulli-trialsに基づく確率論的モデルを使用するが、MsDetectorはヒト染色体からの標識反復データセット（Girgis and Sheetlin、2013）で訓練された隠れマルコフモデルを用いてSSRを見つける。 MREPSは、不一致のエッジをトリミングし、誤ったSSRをフィルタリングするために、いくつかの統計的方法を使用する（Kolpakov et al、2003）。そのようなヒューリスティックアプローチは、包括的でなく（偽陰性の問題）、正確ではない（偽陽性の問題）という犠牲を払って、高速ランタイムを有する。 SSRITとMISAはバックグラウンド正規表現を使用して、与えられたモチーフ長のすべての可能なリピートを検索する（Temnykh、2001; Thiel et al、2003）。この方法は、モチーフ配列の途中で突然終了するリピートを選ぶことができないという欠点を有する。組み合わせアプローチは網羅的かつ正確であるが、通常、非線形の時間複雑さ（典型的にはO(n log n) ）になる（Lim et al。、2013）。 SA-SSRと呼ばれる最近の包括的なアルゴリズムは、線形時間 (O(n)）でSSRを識別するためにsuffix arraysを使用する（Pickett et al、2016）。しかし、実際の実行時間は依然として大きすぎて、大きなゲノムの分析に実際に使用することはできない。

　ここでは、PERFというツールを紹介する。このツールは、リピートセットとの直接的な文字列比較に基づいてSSRを識別するための新しいアルゴリズムを使用する。 PERFは現在の既存の方法よりも数倍高速で、100％正確で包括的でメモリ効率が高い。他のほとんどのメソッドとは異なり、著者らのアルゴリズムは重複するリピートやモチーフの途中で終わるものも見逃すことはない。さらに、PERFはインタラクティブで完全なスタンドアロンのHTMLレポートを生成し、入力ファイルに存在するすべてのSSRのダウンストリーム分析を容易にする。

特徴（Githubより）

• Fast run time, despite being a single-threaded application. As an example, identification of all SSRs from the entire human genome takes less than 7 minutes. The speed can be further improved ~3 to 4 fold using PyPy (human genome finishes in less than 2 minutes using PyPy v5.8.0)
• Linear time and space complexity (O(n))
• Identifies perfect SSRs
• 100% accurate and comprehensive - Does not miss any repeats or does not pick any incorrect ones
• Easy to use - The only required argument is the input DNA sequence in FASTA format
• Flexible - Most of the parameters are customizable by the user at runtime
• Repeat cutoffs can be specified either in terms of the total repeat length or in terms of number of repeating units
• TSV output and HTML report. The default output is an easily parseable and exportable tab-separated format. Optionally, PERF also generates an interactive HTML report that depicts trends in repeat data as concise charts and tables

論文表３と表４で精度とランニングタイムが議論されている。

https://www.ncbi.nlm.nih.gov/pubmed/29121165

インストール

依存

• Compatible with both Python 2 (tested on Python 2.7) and 3 (tested on Python 3.5).

本体　Github

https://github.com/rkmlab/perf

pip install perf_ssr

 または

git clone https://github.com/RKMlab/perf.git 
cd perf 
python setup.py install 

> PERF -h

$ PERF -h

usage: PERF [-h] -i <FILE> [-o <FILE>] [-a] [-l <INT> | -u INT or FILE]

            [-rep <FILE>] [-m <INT>] [-M <INT>] [-s <INT>] [-S <FLOAT>]

            [-f <FILE> | -F <FILE>] [--version]

Required arguments:

  -i <FILE>, --input <FILE>

                        Input file in FASTA format

Optional arguments:

  -o <FILE>, --output <FILE>

                        Output file name. Default: Input file name + _perf.tsv

  -a, --analyse         Generate a summary HTML report.

  -l <INT>, --min-length <INT>

                        Minimum length cutoff of repeat

  -u INT or FILE, --min-units INT or FILE

                        Minimum number of repeating units to be considered.

                        Can be an integer or a file specifying cutoffs for

                        different motif sizes.

  -rep <FILE>, --repeats <FILE>

                        File with list of repeats (Not allowed with -m and/or

                        -M)

  -m <INT>, --min-motif-size <INT>

                        Minimum size of a repeat motif in bp (Not allowed with

                        -rep)

  -M <INT>, --max-motif-size <INT>

                        Maximum size of a repeat motif in bp (Not allowed with

                        -rep)

  -s <INT>, --min-seq-length <INT>

                        Minimum size of sequence length for consideration (in

                        bp)

  -S <FLOAT>, --max-seq-length <FLOAT>

                        Maximum size of sequence length for consideration (in

                        bp)

  -f <FILE>, --filter-seq-ids <FILE>

  -F <FILE>, --target-seq-ids <FILE>

  --version             show program's version number and exit

ラン

 マイクロサテライトを検出する。

PERF -i input.fa -a -o output

• -i      Input file in FASTA format
• -a     Generate a summary HTML report.
• -o     Output file name. Default: Input file name + _perf.tsv

検出可能なサイズなどのパラメータはGithubで詳しく説明されています。

テスト

Ensembl humanのGRCh38のランタイムを調べてみた。

time PERF -i input.fa -a -o output 

デフォルトのパラメータでは、ヒトゲノムの解析が52分で終了した(*1)。

出力

> head output

> head output

1 0 108 AACCCT 108 + 18 TAACCC

1 108 149 AACCCT 41 + 6 AACCCT

1 147 179 AACCCT 32 + 5 CCCTAA

1 172 184 AACCT 12 + 2 CCTAA

1 177 233 AACCCT 56 + 9 CCTAAC

1 231 249 AACCCT 18 + 3 CCCTAA

1 255 290 AACCCT 35 + 5 AACCCT

1 285 326 AACCCC 41 + 6 AACCCC

1 330 353 AACCCT 23 + 3 CCCTAA

1 352 392 AACCCT 40 + 6 ACCCTA

出力フォーマットはBEDに従っている（Githubより）。

"-a"をつけてランすると、htmlレポートも出力される。

 > open output.html #デフォルトブラウザで開く

グラフの表示項目はhtml上で編集できる。

*1 mac pro 2012 (x5690)で実施。tiimeコマンドの"real"で測定。

勘違いして２回紹介してしまいました。両方残しておきます。

引用

PERF: an exhaustive algorithm for ultra-fast and efficient identification of microsatellites from large DNA sequences

Avvaru AK, Sowpati DT, Mishra RK

Bioinformatics. 2018 Mar 15;34(6):943-948.