macでインフォマティクス

macでインフォマティクス

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

メタゲノムデータ間の類似性を計算し可視化する metafast

 

最近、コンピュータ生命科学者たちは、利用可能なショットガンメタゲノミックデータセットの量が驚異的に増加するのを目の当たりにしている。データ分析の次元性を低下させるという課題は、メタゲノムの統計分析の第一の要求である。これには、分類学的および機能的プロファイリング、 リッチさおよび類似性の評価が含まれる。ハイスループットシークエンシングの技術的進歩とコスト削減は、以前に未踏の生態学的ニッチから微生物を調べることを可能にする。前例のない方法でディテールの程度が増加した。最初のメタゲノム研究(Venter et al、2004)以来、平均カバレッジ深度は数桁増加した。単離した微生物の研究から得られた膨大なゲノムデータは、基準に基づくアプローチを開発するための基礎として役立った。続いて、参照ゲノム(例えば、参照ゲノム)との各リードの直接的アライメントよりも巧妙な技術に基づくこのような方法のブームがあった。 Kraken(Wood and Salzberg、2014)(紹介)、CLARK(Ounitら、2015)、FOCUS(Silvaら、2014)、MetaPhlAn2(Truongら、2015)(紹介)。しかし、参照に基づく方法の実際の課題の主なものは、未培養細菌の大部分を含む未知のニッチ共同体になる。参考となる微生物およびウイルスの多数のクレードの代表的なゲノムが欠如している。この問題は、何十年にもわたって徹底的に研究された環境であっても重大なものである。

 急速に蓄積されるデータ量に対処するために開発されたメタゲノム類似性を測定するための手法の1つは、adaptive subsampling(Shamsaddini et al、2014)に基づく(pubmed)。もう1つは、利用可能な参照ゲノムセットの希薄さのために、メタゲノム研究者にとって魅力的であるようであるアラインメントフリーアプローチである。そのような方法は、組成に基づく方法(k-merスペクトル分析(Chatterji et al、2008; Dubinkina et al、2016; Silva et al、2014; Wu and Ye、2011; Wu et al、2016 、Vinga and Almeida、2003)、ニューラルネットワーク(Rasheed and Rangwala、2012)、マルコフモデル(Song et al、2014))がある。しかし、これらの方法にはいくつかの制限がある。メタゲノムの2つ以上のグループ間に差異のある豊富な特徴がこの方法で隠されるか、またはほとんど情報を提供しない。

 類似性を評価するための別のアイデアは、メタゲノムの de novo assembly とそれに続くコンティグの分析(カバレッジ深度に基づく分類と存在量の差の分析)である。しかしながら、このアセンブリは、細菌種の豊富さおよび著しい種内ゲノム変動のために複雑である。これらの問題に対処するためのメタゲノムアセンブリを目的とした特別なアルゴリズムが開発されている(Boisvert et al、2012; Namiki et al、2012; Peng et al、2011; Treangen et al、2013)。特に、対形成類似性(crAss)を推定するために、メタゲノムリードの複数組み合わせたアセンブリが提案された(Dutilh et al、2012)。しかし、リードからコンティグへの完全なアセンブリは、特にpublicに入手可能なメタゲノムデータの急速な増加のために、計算上および記憶上の難度がある。

 著書らは単純化されたメタゲノム・デノボ・アセンブリに基づき、メタゲノムのde Bruijn graphの適応的セグメンテーションを使用して、メタゲノムをコンパクトに表現するためのMetaFastアルゴリズムを開発した。これは、任意の環境からのショットガンメタゲノムを、簡素化されたコンポーネントからなる修正されたDe Bruijnグラフとして表現することを可能にするアプローチである。著者らの方法は、k-merスペクトル分析とアセンブリとの間にあり、これらの2つのアライメントフリーアプローチの最高のものである:すなわち前者の速度と後者の精度のバランス。複数のメタゲノムを入力として、類似性マトリックスを出力する。このアルゴリズムで保存されたメタゲノム成分の次元構造は、微生物の固有の亜種レベルの多様性を反映している。この方法は、計算上効率的であり、特に新規な環境ニッチからのメタゲノムの分析に有望である。MetaFastのパフォーマンスをいくつかの分類学的プロファイラ(Kraken、CLARK、FOCUS、MetaPhlAn2)、クロスアセンブリベースのアルゴリズムcrAssと比較し、シミュレートされたデータと腸内マイクロバイオームと、リアルデータのNew Yorkの地下鉄と湖沼のウイルスメタゲノムで比較した。比較結果は、MetaFastが非常に効率的であり、結果が既存の方法と一致していることを示した。

 

以下のようなフローを取っている(Githubより)。

  1. Assembling short genomic sequences from reads for every metagenome separately (basing on de Bruijn graph).
  2. Constructing one combined de Bruijn graph for all assembled sequences, then searching for connected components in it.
  3. Calculating a characteristic vector for every metagenome with a length equal to the number of connected components.
  4. Cross-comparing metagenomes by calculating the Bray-Curtis dissimilarity matrix based on characteristic vectors.

 

wiki

https://github.com/ctlab/metafast/wiki

 

インストール

依存

  • JRE 1.6 or higher

本体 Github

https://github.com/ctlab/metafast#installation

git clone https://github.com/ctlab/metafast.git 
cd metafast
ant
./out/metafast.sh --version

$ ./out/metafast.sh --version

Fast metagenome analysis toolkit, version 0.1.2 (revision 0863958, 01-Apr-2018)

 

 

またはGithubのリリース(リンク)から実行ファイル metafast.shをダウンロードして実行権をつける (chmod u+x)。

$ /Users/user/local/metafast/out/metafast.sh

Fast metagenome analysis toolkit, version 0.1.2 (revision 0863958, 01-Apr-2018)

 

Usage:     metafast [<Launch options>] [<Input parameters>]

Tool:           matrix-builder

Description:    Builds the distance matrix for input sequences

Input parameters (only important):

-k, --k <arg>                           k-mer size (in nucleotides, maximum 31 due to realization details) (optional, default: 31)

-i, --reads <args>                      list of reads files from single environment. FASTQ, FASTA (MANDATORY)

-b, --maximal-bad-frequency <arg>       maximal frequency for a k-mer to be assumed erroneous (optional, default: 1)

-l, --min-seq-len <arg>                 minimal sequence length to be added to a component (in nucleotides) (optional, default: 100)

-b1, --min-component-size <arg>         minimum component size in component-cutter (in k-mers) (optional, default: 1000)

-b2, --max-component-size <arg>         maximum component size in component-cutter (in k-mers) (optional, default: 10000)

-wr, --without-renumbering              don't renumber samples in the heatmap (optional)

 

Launch options (only important):

-m, --memory <arg>                      memory to use (for example: 1500M, 4G, etc.) (optional, default: 90% of free memory (currently 59.2 Gb))

-w, --work-dir <arg>                    working directory (optional, default: workDir)

-c, --continue                          continue the previous run from last succeed stage, saved in working directory (optional)

-h, --help                              print short help message (optional)

 

To see all parameters and options add --help-all.

To see full documentation visit https://github.com/ctlab/metafast/wiki

 

解析に必要なハードウエアのスペック

  • RAM: metafast requires 2-2.5 times more memory than maximum size of uncompressed FASTQ file to be processed.
  • Hard disk space: metafast requires 25-30% of total size of processed uncompressed FASTQ files.

 

実行方法

 GithubのREADME(トップページ)に載っているテストデータ( meta_test_1.fa, meta_test_2.fa、meta_test_3.fa )をダウンロードする。

metafast.sh -i meta_test_1.fa meta_test_2.fa meta_test_3.fa
  • -i    List of reads files from single environment. FASTQ, BINQ, FASTA files are acceptable, gzip- and bzip2-compressed files are allowed too. Files can be set by bash regexp, for example -i dir/*.fastq or -i `cat filelist.txt`. 
  • -k   K-mer size (in nucleotides, maximum 31 due to realization details). The default value is 31 nucleotide.
  • -l    Minimum sequence length to be added to a component (in nucleotides). The default value is 100 nucleotides.
  • -p   Available processors. By default metafast uses all available processors.
  • -v   Enable debug output.

 

出力ディレクトリに様々なフォルダができる。

f:id:kazumaxneo:20180401212534j:plain

matrices/内にheatmapのPNGファイルができる。

f:id:kazumaxneo:20180401212536j:plain

 Distanceがゼロ(自分自身)は白になり、それから距離が遠いほど青くなるヒートマップがdendrogram付きで出力される。上記の結果は、sample1とsample3のメタゲノムデータがsampleより近い(類似性が高い)ことを示している。

 

 

引用

MetaFast: fast reference-free graph-based comparison of shotgun metagenomic data.

Ulyantsev VI, Kazakov SV, Dubinkina VB, Tyakht AV, Alexeev DG

Bioinformatics. 2016 Sep 15;32(18):2760-7.

 

FASTA分析に使えるpythonライブラリ Goldilocks

 

 Goldilocksは基準を満たす領域のさらなる解析を行うために設計されたPythonパッケージである。パッケージをスタンドアロンスクリプトにインポートするか、コマンドラインツールを使用して使用できる。(一部略)Goldilocksはもともと、複数のサンプルにわたって明示的基準に一致する領域を返すように設計されていた。このパッケージは、以来、より柔軟性があり、GC含有量、標的モチーフの頻度、予め定義されたメトリクスおよび未同定ヌクレオチド(N)などの任意の基準に基づいて、ユーザーが関心のある領域を見つけるために使用することができる。

 

 マニュアル

Welcome to Goldilocks’s documentation! — Goldilocks 0.1.1 documentation

 ここではコマンドラインでの使い方について簡単に説明します。 

インストール

依存

  • numpy
  • matplotlib (for plotting)

To test;

  • tox
  • pytest

For coverage;

依存もこのツールも全てpipで導入できる。

pip install numpy matplotlib tox pytest nose python-coveralls 

本体 Github

GitHub - SamStudio8/goldilocks: Locating genetic regions that are "just right"

pip install goldilocks

$ Goldilocks -h

usage: Goldilocks [-h] [-t TRACKS [TRACKS ...]] [-f {bed,circos,melt,table}]

                  -l LENGTH -s STRIDE [-@ PROCESSES]

                  {gc,ref,motif,nuc} {min,max,mean,median,none} faidx

                  [faidx ...]

 

Wrapper script for Goldilocks library.

 

positional arguments:

  {gc,ref,motif,nuc}

  {min,max,mean,median,none}

  faidx

 

optional arguments:

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

  -t TRACKS [TRACKS ...], --tracks TRACKS [TRACKS ...]

  -f {bed,circos,melt,table}, --format {bed,circos,melt,table}

  -l LENGTH, --length LENGTH

  -s STRIDE, --stride STRIDE

  -@ PROCESSES, --processes PROCESSES

 

ラン

nucコマンド

ATGCをウィンドウサイズ100000bpで5000bp重複しながらカウントする。

goldilocks nuc none --tracks A C G T N -l 100000 -s 50000 -@ 8 input.fa.fai > output 
  • --tracks  TRACKS
  • -f  bed,circos,melt,table
  • -l   LENGTH
  • -s   STRIDE 
  • -@   PROCESSES

$ head output 

chr pos_start pos_end 0_A 0_C 0_G 0_N 0_T

1 1 100000 32149.0 21157.0 17026.0 0.0 29668.0

1 50001 150000 31706.0 20438.0 18182.0 0.0 29674.0

1 100001 200000 31785.0 20762.0 19049.0 0.0 28404.0

1 150001 250000 29765.0 22047.0 21919.0 0.0 26269.0

1 200001 300000 25668.0 25192.0 23868.0 0.0 25272.0

1 250001 350000 24846.0 25724.0 24445.0 0.0 24985.0

1 300001 400000 25632.0 24856.0 24349.0 0.0 25163.0

1 350001 450000 27205.0 23387.0 21910.0 0.0 27498.0

1 400001 500000 29563.0 21851.0 20667.0 0.0 27919.0

 

gcコマンド

GC含量をウィンドウサイズ100000bpで5000bp重複しながらカウント。

goldilocks gc max -l 100000 -s 50000 -@ 8 input.fa.fai > output

 

importして使う場合はオンラインドキュメントを読んでください。

Basic Package Usage — Goldilocks 0.1.1 documentation

配列の検索から、結果を可視化するplotの使い方、ggplot2やCyrcos plotとの連携など詳しく説明されています。

 

引用

Goldilocks: a tool for identifying genomic regions that are ‘just right’

Samuel M. Nicholls, Amanda Clare, and Joshua C. Randall

Bioinformatics. 2016 Jul 1; 32(13): 2047–2049.

 

コマンドライン環境のゲノムブラウザ ASCIIGenome

2019 6/17 インストール追記

 

次世代シーケンシングデータの視覚化は、研究者が結果の質を評価し仮説を生成することを可能にするゲノミクスの基本的な部分である。したがって、ゲノムデータをブラウズするためのいくつかのプログラムは、ゲノミクスコミュニティの間で広く普及しており、例えばIGV(Robinson et al、2011)またはArtemis(Carver et al、2012)が一般的な選択肢である。全てではないにしても、ほとんどのゲノムブラウザはグラフィカルユーザインタフェースGUI)を利用しており、プログラミング経験のない研究者にとって特に適した直感的な環境を作り出す。しかし、多くのバイオインフォマティシャンにとって、コマンドラインインターフェースはGUIよりも好まれ、より効率的なデータ分析が可能になる。コマンドラインインターフェース)は、実行されたコマンドをより細かく制御し、繰り返しの作業を簡単に合理化する。さらに、実行されたコマンドはプロジェクトのドキュメントの一部になることがある。さらに、ゲノムデータはリモートサーバ上に存在することが多く、そのようなデータをローカルに転送することは非実用的ではない。特に、分析者がデータを大雑把な探索する必要がある場合は特にそうである。

 ASCIIGenomeは、ターミナルから直接実行するように設計されたブラウザで、コマンドラインインターフェースを使用している。グラフィカル出力は、ターミナルウィンドウ上でASCII文字でレンダリングされるため、リモートサーバーで実行するのに適している。 UNIXコマンドラインインターフェイスに慣れている研究者は、ASCIIGenomeも同様によく知られているはずである。 SAMTools tview(Li et al、2009)は端末ウィンドウ用に設計されているが、その機能は非常に限られている。一部のブラウザでは、データをバッチ処理(IGV)したり、コマンドラインオプション(Artemis)を受け入れることができるが、GUIの表示やインタラクティブブラウジングに頼っているブラウザもある。ASCIIGenomeは一般的なGUIブラウザと同様の機能を備えているが、コマンドラインで完全に操作して端末ウィンドウにデータを視覚化する。

 

マニュアル

http://asciigenome.readthedocs.io/en/latest/

 

インストール

Github

https://github.com/dariober/ASCIIGenome

実行ファイルはリリースからもダウンロードできる。またはcondaかbrewで導入できる。

#bioconda(linux)
conda install -c bioconda -y asciigenome

#homebrew
brew install https://raw.githubusercontent.com/dariober/ASCIIGenome/master/install/brew/asciigenome.rb

 

$ asciigenome -h

usage: ASCIIGenome [-h] [--batchFile BATCHFILE] [--region REGION] [--fasta FASTA] [--exec EXEC] [--noFormat] [--nonInteractive] [--config CONFIG] [--showMem] [--showTime] [--debug {0,1,2}] [--version] [input [input ...]]

 

DESCRIPTION

Genome browser at the command line.

 

For details see http://asciigenome.readthedocs.io/

 

positional arguments:

  input                  Input files to be displayed: bam, bed, gtf, bigwig, bedgraph, etc

 

optional arguments:

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

  --batchFile BATCHFILE, -b BATCHFILE

                         Bed or gff file of regions to process in batch. Use - to read from stdin. ASCIIGenome will iterate through the regions in this file

  --region REGION, -r REGION

                         Go to region. Format 1-based as 'chrom:start-end' or 'chrom:start' or 'chrom'. E.g. chr1:1-1000

  --fasta FASTA, -fa FASTA

                         Optional reference fasta file.

                         If given, must be indexed, e.g. with `samtools faidx ref.fa`

  --exec EXEC, -x EXEC   Commands to be executed at the prompt. Either a single string, e.g. 'goto chr1 && next && seqRegex ACTG' or a file with one command per line.

  --noFormat, -nf        Do not format output with non ascii chars (colour, bold, etc.) (default: false)

  --nonInteractive, -ni  Non interactive mode: Exit after having processed cmd line args. (default: false)

  --config CONFIG, -c CONFIG

                         Source of configuration settings. It can  be  a  local  file  or  a  tag  matching  a  built-in  configuration:  'black_on_white',  'white_on_black',  'metal'.  If  null,  first  try  to read configuration from file '~/.

                         asciigenome_config'. If this file is missing use a built-in setting. For examples of configuration files see https://github.com/dariober/ASCIIGenome/blob/master/resources/config/ (default: metal)

  --showMem, -sm         Show memory usage. Typically used for debugging only. (default: false)

  --showTime, -st        Show time elapsed to process tracks. Typically used for debugging only. (default: false)

  --debug {0,1,2}        Set debugging mode. 0: off; 1: print exception stack traces; 2: print stack traces and exit. (default: 0)

  --version, -v

 

ラン

2サンプルのbamをリファレンス付きで表示する(3つ以上のbamも表示可能)

asciigenome -fa reference.fa input1.bam input2.bam

f:id:kazumaxneo:20180401170109j:plain

  • f      右向きにウィンドウサイズの1/10だけ移動。
  • ff     右向きにウィンドウサイズの半分だけ移動。
  • +10000   右向きに10000bp移動。(-10000なら左向き)。
  • 750000   ポジション750000にジャンプ(+/-はつけない)。
  • goto chr19   chr19の先頭にジャンプ
  • goto chr19: 3000   chr19のポジション3000にジャンプ
  • r     左向きにウィンドウサイズの1/10だけ移動。
  • rr    左向きにウィンドウサイズの半分だけ移動。
  • Enter   Enterだけ打つと直前のコマンドの繰り返し。
  • zo 2   2倍に縮小。
  • zi 2   2倍に拡大。

h + Enterでヘルプを表示できる。qで終了する。またfastaがなくても開くことができる。bed、gtfも表示できる。

 

オーサーが利点の1つとして挙げているように、GUIのないサーバーにいれておくと重宝すると思います。GUIが利用できる環境でも、IGVは立ち上がりが遅くクラッシュも多いので、特定の変異部位を確認するなどのタスクならASCIIGenomeのほうがはるかに高速かつ簡単だと思います。(起動してポジション手打ちすれば素早く表示できる。(5~10秒))。

 

引用

ASCIIGenome: a command line genome browser for console terminals.

Beraldi D

Bioinformatics. 2017 May 15;33(10):1568-1569.

 

倍数体のfractionation biasを視覚化する FractBias

 

 全ゲノム重複(WGD)などの倍数性事象は、単一の生物体内に2つ以上のゲノムコピーを作成する。重複(サブゲノム)に由来するホモロガスな染色体の全セットは、遺伝子が相同染色体の1つからlossするfractionationと呼ばれる過程で遺伝子欠損を受ける(Langham et al、2004)。fractionationは最終的に、重複前の元の二倍体状態に近い数まで遺伝子数を減少させる(diploidization)。興味深いことに、複製された遺伝子対のいずれかのコピーが失われる可能性も常に同じではない。シロイヌナズナ(Thomas et al、2006)、トウモロコシ(Schnable et al、2011)およびBrassica rapa(Tang et al、2012)を含むいくつかの種において、特定のサブゲノムからの遺伝子欠失のバイアスが観察されている。

 植物において、fractionationはWGDに続く主要な進化メカニズムであるが(Freeling et al、2012)、WGDの影響と結果として生じるfractionationバイアスは、主にfractionationバイアスを特徴付けるための使い易いツールがないため完全には研究されていない。この論文では、倍数体後の重複した合成領域間のバイアス分画を計算し、視覚的に評価するWebベースの自動ツール、FractBiasについて報告する。 FractBiasは、比較ゲノミクス(CoGe)プラットフォームを使用してローカルまたはオンラインで使用できる。

 

 

ラン

SynMapにアクセスする。

CoGe: SynMap

ゲノムを比較したい生物種2つを入力する。

f:id:kazumaxneo:20180401123702j:plain

Analysis optionのタブに切り替える。

f:id:kazumaxneo:20180401123816j:plain

 

Fraction Biasはdefaultでは計算されない。Fraction Biasを計算するには、パラメータを以下のように変更する。

f:id:kazumaxneo:20180401130151j:plain

  1. Synthetic depthをQuota Alignに指定。
  2. Ratio of depthを2ゲノムのpolypoidyの幅で指定。
  3. Fractionation Biasの項目のRunにチェックをつける。
  4. Sliding Window Sizeサイズを調整(defaultは100遺伝子)
  5. 計算を All genesかsyntenic genesに指定(論文のsupplementary Fig.1B)

ランを開始すると、進捗状況が表示される。

f:id:kazumaxneo:20180401122912j:plain

 

 

結果

まず比較した2ゲノムのHarr plotが表示される。

f:id:kazumaxneo:20180401123148j:plain

 

Harr plotの下の方にFraction Bias分析結果が表示される。

 

f:id:kazumaxneo:20180401123136j:plain

 

結果の解釈の仕方については、論文内で詳細が説明されている。Aは遺伝子欠失が多い領域、Bは欠失の少ない領域を表す。

f:id:kazumaxneo:20180401125313j:plain

この他の典型的な例は論文のsupplementary Fig4-6に載っている。論文に書いてあるようにどのゲノムと比較するかが大事になるが、supplementary Table 1に比較例がある。

 

引用

FractBias: a graphical tool for assessing fractionation bias following polyploidy.

Joyce BL, Haug-Baltzell A, Davey S, Bomhoff M, Schnable JC, Lyons E

Bioinformatics. 2017 Feb 15;33(4):552-554.

 

k-merを使いSimple sequence repeats (SSRs) を検索する Kmer-SSR

 

 Simple sequence repeats (SSRs) は、DNA複製、修復、または組換えに起こるミスペアリングやミスのために、少なくとも1つの塩基が何回もタンデムに繰り返されるDNAの短いリピート領域である(Levinson and Gutman、1987)。数十年間、SSRは、短いリピート配列のコピー数の増加によって引き起こされる表現型の相違を決定するために研究されている(Kashi and King、2006)。さらに、SSRは、種の適応度を低下させることなく、定量的遺伝的変異および表現型の差異を説明する(Kashi et al、1997)。 SSR濃度は、異なる種間だけでなく、同じ種内の異なる染色体間でも変化し、配列のヌクレオチド組成を評価することによって説明することができない(Katti et al、2001)。 SSRはDNA複製、組換えおよび修復の特徴的な機能を明らかにするため、生物系の相互作用の研究、次世代配列決定データを用いた再増殖ベースの疾患の研究(Kashi and King、2006)において重要である。

 SSRを識別するために、多くの異なる手法が使用されてきた。この論文では、k-merの使用を提案する。 k-merという用語は、与えられた配列から得られる長さ 'k'の部分配列を指し、一方、k-mer分解は、配列から作ることができる長さ 'k'のすべての可能な部分文字列を指す。 k-mer分解のための使用は、以前にゲノムアセンブリ機械学習(Chikhi and Medvedev、2014; Ghandi et al。、2014)などの例で概説されてきた。 (Hanら、2007)と同様の部分配列を同定するためにk-mersが使用されているが、我々の知る限り、k-mer分解によってSSR同定を試みたことはない。Kmer-SSRはほとんどのSSR識別アルゴリズムよりも速く全てのSSRを見つけることができる。Kmer-SSRのオプションは、100% recallと100% precisionを持ち、指定された長さのSSRをすべて正確に識別する。さらに生物学的に関連するSSRを特定するために、ユーザの入力に基づいてSSRのサブセットを容易に見ることを可能にするいくつかのフィルタも開発した。 Kmer-SSRは、フィルターオプションと組み合わせることで、SSRを他のSSR識別アルゴリズムよりも素早く、ユーザーフレンドリーな方法で正確かつ直感的に識別する。

 

Kmer-SSRの特徴(公式より)

  • Fast run time (linear, O(n), time complexity)
  • Memory efficient (linear, O(n), space complexity)
  • Finds all perfect repeats
  • Simple command-line interface, convenient for scripting and when running on High-Performance Computing (HPC) systems (note: no GUI provided)
  • Easily parsed, tab-delimited output
  • Runs on Linux (not Windows or Mac OS X)

 

インストール

Github

git clone https://github.com/ridgelab/Kmer-SSR.git
cd Kmer-SSR/
make
cd bin/

./kmer-ssr -h

$ ./kmer-ssr -h

./kmer-ssr: illegal option -- h

USAGE: kmer-ssr [-a a1,..,aN] [-A] [-d] [-e] [-h] [-i file] [-l int] [-L int]

                [-n int] [-N int] [-o file] [-p p1,..,pN] [-r int] [-R int]

                [-s s1,..,sN] [-t int] [-v]

 

Find SSRs in FASTA sequence data

 

Input:

 

  -i in.fasta

                        The input file in fasta format. All sequence characters

                        will be converted to uppercase. [default: stdin]

 

                        If your fasta file is compressed, do not use -i. Simply

                        use zcat, bzcat, or a similar tool and pipe it into this

                        program.

 

Output:

 

  -o out.tsv

                        The output file in tab-separated value (tsv) format.

                        Please see `README' column details. [default: stdout]

 

Algorithmic:

 

  -a a1,..,aN

                        A comma-separated list of valid, uppercase characters

                        (nucleotides). Characters not in this list will be

                        ignored. [default=A,C,G,T]

 

  -A

                        Report non-atomic SSRs (e.g., AT repeated 6 times may

                        report an ATAT repeated 3 times or an ATATAT repeated

                        2 times instead).

 

  -e

                        Disable all filters and SSR validation to report every

                        SSR. Similar to: -A -r 2 -R <big_number> -n 2 -N

                        <big_number>. This will override any options set for

                        -n, -N, -r, -R, and -s.

 

  -p p1,..,pN

                        A comma-separated list of period sizes (i.e., kmer

                        lengths). Inclusive ranges are also supported using a

                        hyphen. [default=4-8]

 

  -l int

                        Only search for SSRs in sequences with total length

                        >= l [default: 100]

 

  -L int

                        Only search for SSRs in sequences with total length

                        <= L [default: 500,000,000]

 

  -n int

                        Keep only SSRs with total length (number of

                        nucleotides) >= n [default: 16]

 

  -N int

                        Keep only SSRs with total length (number of

                        nucleotides) <= N [default: 10,000]

 

  -r int

                        Keep only SSRs that repeat >= r times [default: 2]

 

  -R int

                        Keep only SSRs that repeat <= R times [default:

                        10,000]

 

  -s s1,..,sN

                        A comma-separated list of SSRs to search for; e.g.

                        "AC,GTTA,TTCTG,CCG" or "TGA". Please note that other

                        options may prevent SSRs specified with this option

                        from appearing in the output. For example, if -p is

                        "4-6", then an SSR with a repeating "AC" will never be

                        displayed because "AC" has a period size of 2 (and, as

                        it turns out, 2 is not in the range 4-6).

 

Misc:

 

  -d

                        Disable the progress bar that normally prints to

                        stderr. Will automatically be disabled if (a) reading

                        from stdin or (b) writing to stdout without

                        redirecting it to a file.

 

  -h

                        Show this help message and exit

 

  -Q int

                        Max size of the tasks queue [default: 1,000]

 

  -t int

                        Number of threads [default: 1]

 

  -v

                        Show version number and exit

パスの通ったディレクトリに移動しておく。

 

ラン

デフォルト条件での検索。

kmer-ssr -i input.fasta -o output.tsv

 

kmer-ssr -t 4 -p 2-9 -r 3 -R 20 -i input.fasta -o output.tsv
  •  -p    Search for SSRs with a period size of 2, 3, 4, 5, 6, 7, 8, or 9. As examples: `AAAAAA...' and `ACGCAGTTGCACGCAGTTGC...' would not make the cutoff because `A' is shorter than 2 nucleotides and `ACGCAGTTGC' is longer than 9 nucleotides. However, `ACACAC...' and `AATCCTGGTAATCCTGGT...' would be included.
  • -r, -R   The min (-r) and max (-R) times the repeating unit repeats. As examples: `ACAC' and `ACACACACACACACACACACAC ACACACACACACACACACAC' (21 `AC' units) would not make the cutoff. However, `ACACAC' and `ACACACACACACACACACACACACAC ACACACACACACAC' (20 `AC' units) would be included.

 

 

引用

Kmer-SSR: a fast and exhaustive SSR search algorithm.

Pickett BD, Miller JB, Ridge PG

Bioinformatics. 2017 Dec 15;33(24):3922-3928.

 

 

 

 

 

ゲノムワイドにマイクロサテライトを高速検索する PERF

 

 Repetitive DNA はゲノムのかなりの割合を構成し、i)散在したリピートまたは転移可能なエレメントと ii)タンデムリピートの2つのカテゴリーに大別できる(Kumar et al、2010)。繰り返しモチーフの長さに依存して、タンデムリピートは、サテライト(> 100nt)、ミニサテライト(10〜30nt)およびマイクロサテライト(1-6nt)として分類される。 Simple Sequence Repeats(SSR)またはShort Tandem Repeats(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はリピートを同定するためにBernoulli-trialsに基づく確率論的モデルを使用し(Benson、1999)、いっぽうMsDetectorはヒト染色体からのリピートデータセットで訓練された隠れマルコフモデルを用いてSSRを見出す(Girgis and Sheetlin、2013)。 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というツールを紹介する。PERFは、リピートセットとの直接的な文字列比較に基づいてSSRを識別するための新しいアルゴリズムを使用する。 PERFは現在の既存の方法よりも数倍高速で、100%正確で包括的でメモリ効率が高い。他のほとんどのメソッドとは異なり、著者らのアルゴリズムは繰り返しのリピートやモチーフの途中で終わるもの見逃すことはないとされる。PERFはインタラクティブで完全なスタンドアロンのHTMLレポートを生成する機能も持ち、入力ファイルに存在するすべてのSSRのダウンストリーム分析を容易にする。

 

インストール

Github

https://github.com/rkmlab/perf

pip install perf_ssr

$ 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

uesaka-no-Air-2:~ kazumaxneo$ 

 

ラン

EnsemblからGRCh38のchr20をダウンロードして検索してみる(リンク)。

perf -i Homo_sapiens.GRCh38.dna_rm.chromosome.20.fa -a
  • -i    Input file in FASTA format
  • -a   Generate a summary HTML report

Homo_sapiens.GRCh38.dna_rm.chromosome.20_perf.tsvが出力される。BED形式に乗っ取っている。

 

$ head Homo_sapiens.GRCh38.dna_rm.chromosome.20_perf.tsv 

20 60181 60193 AATAT 12 - 2 TATAT

20 67482 67494 AAAGGC 12 - 2 GCCTTT

20 68650 68662 AAACTG 12 - 2 TTTCAG

20 82666 82679 ACTCAG 13 + 2 CTCAGA

20 83500 83514 AGAGGC 14 - 2 TGCCTC

20 84150 84162 AAGGAT 12 - 2 CCTTAT

20 84290 84302 AAGG 12 + 3 AAGG

20 84606 84621 AAAAC 15 + 3 AAAAC

20 88061 88073 AAGTG 12 - 2 TCACT

20 88562 88574 AAATAC 12 + 2 CAAATA

説明

f:id:kazumaxneo:20180401013300j:plain

公式より転載。

 

-aをつけるとhtmlレポートも出力される。

f:id:kazumaxneo:20180401012924j:plain

f:id:kazumaxneo:20180401012926j:plain

f:id:kazumaxneo:20180401013547j:plain

f:id:kazumaxneo:20180401013552j:plain

 

引用

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.

 

マイクロサテライトの高速検索を行うGUIツール Krait

 

 一般にsimple sequence repeats(SSR)またはsimple tandem repeats(STR)とも呼ばれるマイクロサテライトは、1〜6bpの単位長の短いタンデム反復DNA配列である。マッピングや集団遺伝学、法医学検査および系統解析(Ellegren 2004; Vieira et al、2016)における強力な分子マーカーとして広範囲に利用されてきた。マイクロサテライト突然変異は、多数のヒトの遺伝病に関与し、様々な調節機構および進化において重要な役割を果たす組織内および世代にわたる動的プロセスである(Kelkar et al、2011)。例えば、マイクロサテライトは主に家庭犬の形態進化とvoles(ハタネズミ)の独特の社会行動に関与していると報告されている(Merkel and Gemmell、2008)。

 ゲノムデータの利用可能性が急速に高まっているため、コストと労働集約的な従来のアプローチを使用する代わりに、コンピュータツールに基づく新しいインシリコアプローチがマイクロサテライトマーカーを抽出するために広く使用されている。MISA (Thiel et al、2003)、SciRoKo(Kofler et al、2007)、msatcommander(Faircloth、2008)、MSDB(MSDB)などの多数のアルゴリズムバイオインフォマティクスソフトウェアがマイクロサテライト検索と調査のために開発されている。それらのほとんどはゲノムワイドマイクロサテライト検索において優れた性能を示す。一般的に、ゲノム内のマイクロサテライトの探索は、それらの分布および存在量を研究するために使用されるだけでなく、その後の生物学的分析のための有用で安定なマーカーを得るために使用される。残念なことに、動植物ゲノムは非常に大きく、膨大な数のマイクロサテライトを含んでいるため、プライマー設計の要件を満たす巨大なデータセットからのマイクロサテライトマーカーをスクリーニングすることは難題になる。例えば、ヒトゲノムは、200万個のマイクロサテライト遺伝子座を含む(Subramanian et al、2003)。さらに、以前に開発されたすべての検索ツールは、生物学者が下流の分析を行い、マイクロサテライト検索結果を視覚化することを制限する独自の出力形式を持っている(論文 supplementary table S1)。

 この論文では、利用可能なツールの限界を克服しようとする、マイクロサテライトのゲノムワイドな調査のためのユーザーフレンドリーなグラフィックインターフェイスを備えた堅牢で超高速なツールKraitを紹介する。 Kraitは、完全保存、または不完全な複合マイクロサテライトの検索を可能にし、その後のプライマー設計のためにマイクロサテライトマーカーをスクリーニングすることを可能にする。

 

マニュアルPDF

krait/Manual.pdf at master · lmdu/krait · GitHub

 

インストール

Github リリースよりmac版をダウンロードする。

Krait-0.10.2-mac.tar.xzを解凍するとappファイルができるので、Applicationフォルダ(command + shift + A)に移動しておく。

f:id:kazumaxneo:20180401000408j:plain

 

ラン

起動直後のウィンドウ。perfect microsatellites (SSRs)、compound SSRs、imperfect SSRs、A variable number tandem repeat (VNTR)、primer設計、statistics解析ができる。

f:id:kazumaxneo:20180401000637j:plain

 

f:id:kazumaxneo:20180401001404j:plain

 

FASTAファイルをインポートする。

f:id:kazumaxneo:20180401000644j:plain

NCBIから直接インポートすることも可能。

f:id:kazumaxneo:20180401002239j:plain

 

perfect microsatellites (SSRs)の検索は左端のSSRsボタンをクリックする。

GRCh38のchr20検索結果。

f:id:kazumaxneo:20180401002738j:plain

 

その隣のをクリックすれば、パラメータ変更が可能になっている。

f:id:kazumaxneo:20180401001359j:plain

パラメータ変更画面。

f:id:kazumaxneo:20180401001404j:plain

 

compound SSRs

f:id:kazumaxneo:20180401003249j:plain

 

imperfect SSRs

f:id:kazumaxneo:20180401003331j:plain

 

A variable number tandem repeat (VNTR

f:id:kazumaxneo:20180401003642j:plain

 

Locateボタンでは、検出された部位がgene内のUTRか、exonかなどを色分けできる。使用するには、リファレンスのgtfファイルを指定する。

解析後のSSRs一覧。色がついている。色はLocateから確認できる。

f:id:kazumaxneo:20180401013808j:plain

 

Primerボタンでは、PCR増幅するためのprimerを設計できる。SSRs解析で見つかった1のSSRのprimerを設計する。1にチェックをつける。

f:id:kazumaxneo:20180401004653j:plain

 PrimerをクリックするとPrimerが設計される。

f:id:kazumaxneo:20180401004557j:plain

 

検索結果のまとめのレポートを作成するには、Statisticsをクリックする。

f:id:kazumaxneo:20180401004057j:plain

f:id:kazumaxneo:20180401004059j:plain

f:id:kazumaxneo:20180401004102j:plain

f:id:kazumaxneo:20180401004104j:plain

f:id:kazumaxneo:20180401004107j:plain

このように解析が終わったデータ分が表示される。 

 

引用

Krait: an ultrafast tool for genome-wide survey of microsatellites and primer design.

Du L, Zhang C, Liu Q, Zhang X, Yue B, Hancock J

Bioinformatics. 2018 Feb 15;34(4):681-683.