macでインフォマティクス

macでインフォマティクス

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

FASTAアライメントからSNP distance matrixを出力する snp-dists

 

snp-distsは、Torsten Seemannさん (GIthub) がGithubで公開されている、配列間のSNP数を計算して、行列出力するツール。

 

 

インストール

mac os10.14のminiconda3-4.3.30環境でテストした。

ビルド依存

  • snp-dists is written in C to the C99 standard and only depends on zlib.

本体 Github

github.com

#Anaconda環境でcondaを使い導入
conda install -c bioconda -c conda-forge snp-dists

#またはhomebrewを使い導入
brew install brewsci/bio/snp-dists

> snp-dists -h

$ snp-dists -h

SYNOPSIS

  Pairwise SNP distance matrix from a FASTA alignment

USAGE

  snp-dists [options] alignment.fasta[.gz] > matrix.tsv

OPTIONS

  -h Show this help

  -v Print version and exit

  -q Quiet mode; do not print progress information

  -a Count all differences not just [AGTC]

  -k Keep case, don't uppercase all letters

  -c Output CSV instead of TSV

  -b Blank top left corner cell instead of 'snp-dists 0.6'

URL

  https://github.com/tseemann/snp-dists (Torsten Seemann)

——

 

テストラン

以下のような4つの配列のアライメントから、総当たりでペアワイズSNP distanceを計算し、matrixで出力する。

cat test/good.aln

$ cat test/good.aln 

>seq1

AGTCAGTC

>seq2

AGGCAGTC

>seq3

AGTGAGTA

>seq4

TGTTAGAC

 

snp-distsを実行

snp-dists test/good.aln > distances.tab

 

出力のSNP distance matrixをexcelで開いてみる。

f:id:kazumaxneo:20181119141852j:plain

数値は塩基置換の数を表す。当然その配列自身で比較すると0になる。もっとも距離があるのはseq2とseq4、及びseq3とseq4。

 

本ツールはindelによるギャップ(-)は計算に入れません。

引用

https://github.com/tseemann/snp-dists

 

Minhashを使い、genomic DNA / proteinを高速比較する sourmash

 

 sourmashは、ゲノムデータのMinHash sketchesを作成、比較、操作するためのツールボックスである。MinHash sketchは、大規模なDNAまたはRNAシーケンスコレクションの"signatures"を保存し、Jaccard indexを使用してそれらを比較または検索するための軽量な方法を提供する。 MinHash sketchは、サンプルを同定し、類似のサンプルを見出し、共有配列を有するデータセットを同定し、系統樹を構築するために使用することができる(Ondov et al、2015)。sourmashはコマンドラインスクリプトPythonライブラリ、MinHashスケッチ用のCPythonモジュールを提供する。 

 

sourmash紹介

https://sourmash.readthedocs.io/en/latest/

document

https://sourmash.readthedocs.io/en/latest/tutorials.html

Using sourmash: a practical guide

https://github.com/dib-lab/sourmash/blob/master/doc/using-sourmash-a-guide.md

A sourmash tutorial(一番説明が丁寧)

2017-dibsi-metagenomics/sourmash.md at master · dib-lab/2017-dibsi-metagenomics · GitHub

 

インストール

mac os10.14のminiconda3-4.3.30環境でテストした。

本体 Github

#Anaconda環境で導入、バージョン指定しないとversion2が入る。
conda install -y -c bioconda

> sourmash -h

sourmash -h

usage: sourmash <command> [<args>]

 

Commands can be:

 

   compute <filenames>         Compute signatures for sequences in these files.

   compare <filenames.sig>     Compute distance matrix for given signatures.

   search <query> <against>    Search for matching signatures.

   plot <matrix>               Plot a distance matrix made by 'compare'.

 

   import_csv                  Import signatures from a CSV file.

.

 

work with RNAseq signatures

 

positional arguments:

  command

 

optional arguments:

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

> sourmash search -h

$ sourmash search -h

# running sourmash subcommand: search

usage: sourmash [-h] [--threshold THRESHOLD] [-k KSIZE] [-f]

                query against [against ...]

 

positional arguments:

  query

  against

 

optional arguments:

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

  --threshold THRESHOLD

  -k KSIZE, --ksize KSIZE

  -f, --force

> sourmash compute -h

$ sourmash compute -h

# running sourmash subcommand: compute

usage: sourmash [-h] [--protein] [--input-is-protein] [-k KSIZES]

                [-n NUM_HASHES] [-f] [-o OUTPUT] [--email EMAIL]

                filenames [filenames ...]

 

positional arguments:

  filenames

 

optional arguments:

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

  --protein

  --input-is-protein

  -k KSIZES, --ksizes KSIZES

                        comma-separated list of k-mer sizes

  -n NUM_HASHES, --num-hashes NUM_HASHES

                        number of hashes to use in each sketch

  -f, --force

  -o OUTPUT, --output OUTPUT

  --email EMAIL

> sourmash compare -h

$ sourmash compare -h

# running sourmash subcommand: compare

usage: sourmash [-h] [-k KSIZE] [-o OUTPUT] signatures [signatures ...]

 

positional arguments:

  signatures

 

optional arguments:

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

  -k KSIZE, --ksize KSIZE

  -o OUTPUT, --output OUTPUT

> sourmash plot -h

$ sourmash plot -h

# running sourmash subcommand: plot

usage: sourmash [-h] [--pdf] [--labels] [--indices] [--vmax VMAX]

                [--vmin VMIN]

                distances

 

positional arguments:

  distances    output from 'sourmash compare'

 

optional arguments:

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

  --pdf

  --labels

  --indices

  --vmax VMAX

  --vmin VMIN

 

 

実行方法

 

ここではsourmash computeのコマンドのみ紹介する。fast-GePのテストゲノムを使う(Github)。

cd fast-GeP-master/Examples/E.faecalis/input_files/

sourmashを実行

sourmash compute *.fa
sourmash compare *.sig -o e_faecalis
sourmash plot --pdf --labels e_faecalis

出力 

f:id:kazumaxneo:20181117215516p:plain

f:id:kazumaxneo:20181117215514p:plain

引用

sourmash: a library for MinHash sketching of DNA

C. Titus Brown, Luiz Irber

Journal of Open Source Software, 1(5), 27

 

関連ツール


 

移動履歴を学習し、移動をナビゲートする autojump

2018 11/17 分かりにくい文章を修正

 

autojumpは、これまでの移動結果をウエイトをつけて記憶し(学習)、補完機能によって移動を助けたり、ファイラーへの表示を助けるcdのパワーアップ版コマンド。

 

wiki

https://github.com/wting/autojump/wiki

 

autojumpに関するツイート

他にもいくつかcdの拡張コマンドはあるようですね。 

 

 インストール

mac os10.14でテストした。

依存

  • Python v2.6+ or Python v3.3+
  • Supported shells
    • bash - first class support
    • zsh - first class support
    • fish - community supported
    • tcsh - community supported
    • clink - community supported
  • Supported platforms
    • Linux - first class support
    • OS X - first class support
    • Windows - community supported
    • BSD - community supported
  • Supported installation methods
    • source code - first class support
    • Debian and derivatives - first class support
    • ArchLinux / Gentoo / openSUSE / RedHat and derivatives - community supported
    • Homebrew / MacPorts - community supported

#hokmebrew
brew install autojump

#またはソースから
git clone git://github.com/wting/autojump.git
cd autojump
./install.py or ./uninstall.py

autojump -h

$ autojump -h

usage: autojump [-h] [-a DIRECTORY] [-i [WEIGHT]] [-d [WEIGHT]] [--complete]

                [--purge] [-s] [-v]

                [DIRECTORY [DIRECTORY ...]]

 

Automatically jump to directory passed as an argument.

 

positional arguments:

  DIRECTORY             directory to jump to

 

optional arguments:

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

  -a DIRECTORY, --add DIRECTORY

                        add path

  -i [WEIGHT], --increase [WEIGHT]

                        increase current directory weight

  -d [WEIGHT], --decrease [WEIGHT]

                        decrease current directory weight

  --complete            used for tab completion

  --purge               remove non-existent paths from database

  -s, --stat            show database entries and their key weights

  -v, --version         show version information

 

Please see autojump(1) man pages for full documentation.

または"j"でも動く(jに別のaliasを割り当ててなければ)。

 

使い方

例1

一度は端末から該当ディレクトリにアクセスしておく必要がある。~/Desktopにfastqというディレクトリを作成して中に入る。

mkdir ~/Desktop/fastq
cd ~/Desktop/fastq

現在のパスを確認

$ pwd

/Users/kazuma/Desktop/fastq

テストのため、別のパス、例えばmusicに移動しておく。

cd ~/music

autojumpを試す。フルネームは手間なので、autojumpのラッパー"j"を使う (jumpするjと覚える)。

j fastq

これで~/Desktop/fastqにジャンプした。autojumpをインストール後、一度アクセスしてfastqというディレクトリ名をautojumpが記憶しているため、フルパス指定なしでジャンプすることが可能になっている。

$ pwd

/Users/kazuma/Desktop/fastq

ちゃんと移動できている。

 

例2

該当文字を含むディレクトリが複数あり、いずれにもアクセスしたことがあるなら、部分一致の状態でタブキーx2を打つと候補一覧が番号付きで表示される。その状態から該当番号を打てばジャンプできる。例えば2つのディレクトリreference_human、reference_bacteriaを作ってテストしてみる。

#作って一度アクセスする
mkdir ~/genome/reference_human/ && cd ~/genome/reference_human
mkdir ~genome/reference_bacteria && cd ~/genome/reference_bacteria

#テスト前に適当な場所に移動。ここではroot直下に移動
cd /

準備できたのでテストする。refまで打ってみる(候補は増えるが、"re"でも、"r"一文字だけでもautojumpは機能する)。

j ref

タブキーx2で候補を表示。

$ j ref__

ref__1__/Users/kazuma/genome/reference_human     

ref__2__/Users/kazuma/genome/reference_bacteria 

reference_bacteriaに飛びたいので"2"を打ち、タブキーを1度押す。2がパスに変換される。そのままエンターキーを打てばジャンプする。

 j /Users/kazuma/genome/reference_bacteria 

/Users/kazuma/genome/reference_bacteria

パス確認

$ pwd

/Users/kazuma/genome/reference_bacteria

root直下からジャンプできている。非常に便利。

 

例3

ジャンプする代わりに、ディレクトリのウィンドウを開くこともできる(mac osのopenコマンドに近いが( "open .")、こっちの方が便利)。

jo music

joコマンドには候補表示機能はありません。パス情報が不完全だとデータベースにある中でウエイトが一番高いパスに飛びます(データベースやウエイトは"j -s"で確認)。

f:id:kazumaxneo:20181117112654p:plain

 

フルパスを忘れたなら、例2のようにjコマンドでまず候補を表示、確定し、それからCtrl + aで先頭にバック、joに変えた方が早そうですね。

 

 

データベースには消去したディレクトリも残っているので、補完すると消したパスも表示され視認性が悪くなります。邪魔になってきたら"autojump --purge <directory>"で消して下さい。--increaseや--decreaseで表示のウエイトを変えることもできます。

引用

GitHub - wting/autojump: A cd command that learns - easily navigate directories from the command line

 

参考HP


 

de bruin graphにリードをマッピングする BGREAT

 

 次世代シーケンシング技術(NGS)は、シーケンシングされたゲノムの生成を大幅に加速した。しかしながら、これらの技術は、依然として染色体当たり単一の配列を提供することができないままである。代わりに、それらは大量かつ冗長なリードセットを生成し、各リードは全ゲノムの一部である。この冗長性のために、リード間のオーバーラップを検出し、それらを一緒にアセンブルしてターゲットのゲノム配列を再構成することが可能である。

 今日でさえ、リードのアセンブリは複雑な作業のままであり、これはどのようなソフトウェアも一貫してうまく働かないためである[論文より ref.1]。アセンブリの問題は以前からずっと計算上困難なままであり、より正確にはNP困難であることが示されている[ref.2]。実用上の限界は、ゲノムの構造(正確に解読できるより長いリピート)およびシーケンシングバイアス(不均一なカバレッジおよびシークエンシングエラー)の両方から生じている。適用されたソリューションは、リードをアセンブリgraphとして表現する: graphのパスに沿ったラベルがシーケンスをエンコードする。現在のところ、大部分のアセンブラは、2種類のgraphに依存している。第2世代のシーケンシング技術によって生成されたショートリードのde Bruijn graph(DBG)[ref.3]、またはCelera Assembler [ref.4])とそのバリエーション、string graph [ref.5]のように。次に、アセンブリアルゴリズムヒューリスティックを使用してgraphを探索し、いくつかの経路を選択してそのシーケンスを出力する。これらのヒューリスティックのために、得られた配列の集合(コンティグと呼ばれる)は、シーケンシングエラー、ゲノムの変異およびリピートによって、生成されるgraphは複雑なパターンになり断片化される。コンティグセットはめったに満足のいくものではなく、通常短いコンティグを破棄するなどして後処理される。

 リードセットを分析するための最も頻繁なコンピュータタスクは、それらをリファレンスゲノム上にマッピングすることである。リファレンスゲノムが一連の配列の形の場合、リードをマップすることができるツールはたくさんある(例えば、BWA [ref.6]およびBowtie [ref.7])。Finishedゲノム配列にマッピングすることの目的は、アライメントできるかどうか、そしてできる場合はどのポジションにアライメントされるか、である。これには、たいてい、リード配列とゲノム配列の間で小さな編集距離またはハミング距離を許可するsemi-globalなアラインメント手順のヒューリスティックで行われる。リードマッピングプロセスは、低マッピング率の領域に苦しんでいる[ref.8]。リピートゲノム領域では、リードのマッピングが複数一致するため、正確にマッピングされないかもしれない。ゲノムがgraphとして表されるとき、各リピート領域出現は分解されて複数一致問題は制限されるため、マッピング問題は低減される。

 リファレンスがFinishedゲノム配列ではなく、重複したコンティグセットである場合、状況は異なる。マッピングは、リードがゲノム内に見出されるかどうかを正確に決定し得るが、複数の位置は、例えば、いくつかの真の位置が存在するかどうかを判断するには不十分であり得る。逆に、不完全なマッピングは、不完全なアセンブリまたは後処理中のいくつかのコンティグの除去に起因する可能性がある。そのような場合、コンティグセットの代わりに、アセンブリgraphを(偏っていないおよび/またはより完全な)リファレンスとして考えることは興味深いと我々(著者ら)は主張する。次いで、このグラフ経路のマッピングは、コンティグのセットに対するマッピングを補完するために必要である。これがBGREATの設計と実装を促した。

  ここでは、graph上にリードをマッピングする問題を探究する。リードをsequence graph( (a generic term meaning a graph representing sequences along its paths) )にマッピングまたはアライメントすることは、既にいくつかの文献、例えばアセンブリ、リードのエラー訂正、メタゲノミクス分野で検討されてきたことである。

 アセンブリ関連では、一旦DBGが構築されると、リードをgraphにバックマッピングすることで、サポートされていないパスを削除したり、エッジのカバレッジを計算したりすることができる。著者らが知る限り、この作業のための実用的な解決法は設計されていない(論文執筆時点)。 Ceruleanアセンブラ[ref.9]はこの可能性について言及しているが、通常のアセンブリ配列へのアラインメントしか実行しない。 Allpaths-LG [ref.10]も、第3世代シーケンシング技術のエラーの多いロングリードを使用してリピートを解決する同様のタスクを実行するが、その手順は、DBG上の任意のリードセットのマッピングに適するほど一般的ではない。理論的な観点からは、問題は、NP困難なリード通過問題(オイラーのスーパーパス問題[2,11]とも呼ばれる)であり、DBGでリードのコヒーレントパスを見つけることからなる( [ref.5]で定義されているような一連のリード配列として表現できるパス)。 SPADES [ref.12]アセンブラは、構築中に使用されるパスを追跡することによってDBGに対するリードを通す。これには相当量のメモリが必要である。ここでは、グラフの作成に使用されるリードやその他のリードのいずれかをNGSリードの任意のソースにマップすることを目指して、より一般的な問題を提案する(De Bruijn Graph Read Mapping Problem(DBGRMPと名付ける))。

 最近では、ショートリードを使ったロングリードのハイブリッドエラー訂正は、第3世代のシーケンシング技術を活用するための重要なステップとなっている。エラー訂正ツールLoRDEC [ref.13]は、ショートリードのDBGを構築し、動的プログラミングアルゴリズム(著者らの目的にとっては遅い)を用いてそれらの編集距離を計算することによってDBGの各行に対して各ロングリードをアライメントする。ショートリードの訂正ツールでは、シーケンシングエラーを訂正するために評価するためにk-merスペクトラムを使うツールが、DBGの確率的または正確な表現をリファレンスとして使ってる[ref.14,15]。

 メタゲノミクス関連では、Wang et al [ref.16]が、closely relatedな複数の細菌種からなるDBG上にリードをマッピングすることによって、サンプルの分類組成を推定した。実際、graphはこれらのゲノムの類似領域を崩壊させ、冗長マッピングを回避する。彼らのツールは、DBGのunitigsのランダムな連結配列にBWAを使用してリードをマッピングする。したがって、graphのいくつかの連続したノード(ER: il y a un pb ce n’est pas vrai))上にリードをアライメントすることはできない。同様に、いくつかの著者は、関連するゲノムを、より反復性の少ない単一のDBGに保存することを提案している[17,18,19]。しかし、これらのツールのほとんどは、hihgly closely relatedな配列に適用した場合にのみ効率的であり、結果としてフラットなgraphになる。 BlastGraphツール[19]は、graphのリードマッピングに特化しているが、実際のgraphでは使用できない(論文の結果セクションを参照)。

 ここでは、De Bruijn graph上のリードマッピングを正式化し、それがNP完全であることを示す。次に、パイプラインGGMAPを提示し、DBGの分岐パス上のリードをマップできる新しいツールBGREAT (Section GGMAP: a method to map reads on de Bruijn Graph)に入る。効率性のために、BGREATは、巨大なシーケンシングデータセットまでスケールアップするヒューリステイックスを採用している。セクションの結果では、マッピング能力と効率の面でGGMAPを評価し、それをアセンブルされたコンティグのマッピングと比較する。最後に、GGMAPの限界と利点について議論し、今後の作業の方向性について述べる(Section Discussion)。

 

インストール

ubuntu14.04でテストした(docker使用。ホストOS: mac os10.12)。

本体 Github

git clone https://github.com/Malfoy/BGREAT.git
cd BGREAT/
make
https://github.com/Malfoy/BCALM.git

>./bgreat

$ ./bgreat 

-r read_file

-k k_value (30)

-g unitig_file (unitig.dot)

-m n_missmatch (2)

-t n_thread (1)

-e effort put in mapping (2)

-f path_file (paths)

-a not_aligned_file (notAligned.fa)

-p to align on paths instead of walks

-q for fastq read file

-c to output corrected reads

 

BCALM

#v1 (DSKも必要になる)
git clone https://github.com/Malfoy/BCALM.git
cd BCALM/
make

#またはcondaで導入(version2.2)(github)
conda install -y -c bioconda bcalm

>./bcalm

$ ./bcalm   

usage: <input> [output.dot] [minimizer length]

Note: maximum minimizer length (minimizer length = 10) requires that you type 'ulimit -n 1100' in your shell prior to running bcalm, else the software will crash

 

DBGを出力するならDSKなどのk-mer頻度カウンタも必要になる。

conda install -y -c bioconda dsk

 

実行方法

1、k-merカウント及びde Bruijn graph構築(BCALM v1)

#BCALM v1を使う場合
#k-merカウント
dsk -file test.fa -kmer-size 31 -abundance-min-threshold 3 -out out
dsk2ascii -file out.h5 -out k-mer_count
bcalm k-mer_count output.dot 10

> k-mer_counthead

$ k-mer_counthead

bash: k-mer_counthead: command not found

kazu@6a393ab918c4:/data$ head k-mer_count 

AAAAAAAGATTTGATTGCGCCGATGGTTTTT 23

AAAAAATCATTGGCACGTTTGGCAACTGTGC 9

AAAAAAGATTTGATTGCGCCGATGGTTTTTA 23

AAAAAAGTGCGCACATGTGCGCACTTTTCTG 12

AAAAACAGCCTGAACTATATGACTTTAACGG 15

AAAAACCAGCGGCGATCGTCCTTAACCTACT 7

AAAAACGTGGAGATCCTCAGGAGTTTTACCT 14

AAAAATCATTGGCACGTTTGGCAACTGTGCT 9

AAAAAGATTTGATTGCGCCGATGGTTTTTAC 23

AAAAAGTGCGCACATGTGCGCACTTTTCTGC 12

 

or

de Bruijn graph構築(BCALM v2 Github)

bcalm -in input.fq -kmer-size 21 -abundance-min 3

 .gfaファイルができる

f:id:kazumaxneo:20181116233748j:plain

 

2、 de Bruijn graphにマッピング。シーケンスデータをfasta形式、またはfastq形式("-q"をつける)で指定する。

bgreat -r reads.fa,reads2.fa -k 21 -t 12 -m 1 -g unitig.fa \
-o no_overlap -a not_aligned -p path

出力の詳細はGithubを参照してください。

 

引用
Read mapping on de Bruijn graphs
Limasset A, Cazaux B, Rivals E, Peterlongo P
BMC Bioinformatics. 2016 Jun 16;17(1):237

 

アセンブリ配列を使って全ゲノムMLST (wgMLST) を行い、アレルプロファイルから系統を比較・再構成する fast-GeP

2018 11/16 tips追記 

 

 Multilocus sequence typing(MLST)などの遺伝子ベースのタイピング法は、バクテリアpopulationsのゲノム研究のための「ゴールドスタンダード」である(Maiden et al、2013; Sheppard et al、2012)。大量の全ゲノムシーケンシング(WGS)データが入手可能になるにつれて、精査中の遺伝子の数は典型的に7 MLST(Maiden、2006; Sheppardら、2012)から53リボソームMLST(Jolley et al、 2012)、1000を超える全ゲノムMLST(wgMLST)(Cody et al、2013; Llarena et al、2017; Zhang et al、2015)に急速に拡大した。

 バクテリアWGSデータをアレルプロファイルに変換するためのいくつかのバイオインフォマティクスツールが開発されている。 PubMLSTデータベース(https://pubmlst.org)によって使用されるBIGSdbは、遺伝子ごとにアレルシーケンスデータベースを検索することによってアレルを割り当てる(Jolley and Maiden、2010)。 EnteroBase(https://github.com/EnteroBaseGroup/EnteroBase)はまた、リファレンス配列データベースを検索するためにBLASTおよびUSEARCHを使用する予め定義されたスキームに依存する。 BionumericsおよびRidom SeqSphere +のような商業的プログラムは、あらかじめ定義されたアレルデータベース(例えば、http://www.cgmlst.org/ncs)を検索し、アレル番号をゲノムに割り当てるために、遺伝子ごとに同様のアプローチを用いるようである。

 可搬性、スケーラビリティ、高解像度、明白な性質(Cody et al、2013; Maiden et al、2013; Sheppard et al、2012)などのMLSTのようなタイピング方法の利点にもかかわらず、さらなるアプリケーションを制限するボトルネックがある(Cody et al、2013; Llarena et al、2017; Sheppard et al、2012)にMLSTのようなタイピング法が用いられている。例えば、ほとんどのMLSTライクなタイピング方法は、スキームで定義されているバリエーションのみを検出することができる(Sheppard et al、2012)。既存のスキームは、いくつかのよく研究されたバクテリア種に高度に偏向されており、あまり研究されていない種は、スキームが不十分であるか、全く利用できないスキームを持っている。

 速度は、特に大量のWGSデータまたは距離が離れた分離株を処理する場合の別の制限になる(Cody et al、2013; Llarena et al、2017; Taboada et al、2017)。 GePを例にとると、distantly relatedな単離株同士の比較では、集中的なhomologous lociを探すために計算集約的なプログラムBLASTXを頻繁に呼び出すため、実行時間が劇的に増加する。対照的に、疫学調査では、疫学的に関連した単離株をアウトブレイクとは無関係のものから迅速に特定し、適時かつ適切な介入を確保することが重要である(Cody et al、2013)。

  Gene-by-gene allele calling approachは、MLST分析に由来し、元のデータ検索プロセスを反映する(Maiden 2006; Maiden et al、2013; Sheppard et al、2012)。典型的には、遺伝子配列は、PCRによって1つずつ増幅され、サンガー法によってシーケンシングされる。次いで、これらの配列を既存のアレル配列と個別に比較して完全一致を同定する(Maiden 2006; Maiden et al、2013; Sheppard et al、2012)。現在のWGS技術は、一回の操作で1株の複数遺伝子シーケンシングを可能にする。従って、現在の非効率的なアレルコール方法はリプレースされる必要がある。

方法

 手短に言えば、各単離株の全ゲノム配列を連結し、単一のクエリにして、すべての座位のすべての既知アレル配列からなるアドホックなアレルデータベースを作成してサーチする。より遠く離れたオーソロガス遺伝子を同定するため、BLASTXを使いWGSヌクレオチドを翻訳してアレルプロテインデータベースに対してサーチを行い、一次アレルコールを行った。リファレンスアミノ酸配列に完全にアライメントされる全てのオーソログ遺伝子は、後の使用のために抽出される。ゲノム中のマルチコピー遺伝子の場合、遺伝子座はパラログとしてマークされ、分析から除外された。より多くのオーソログ遺伝子を見つける機会を増やすために、フレームシフト突然変異を含むアレルの場合のように、最初の検索で省略されたアレルを見つけるために相補的なヌクレオチドデータベース検索がその後に実行された。

 新しいアルゴリズムを実装し、Fast-Genome Profiler(Fast-GeP)という新しいプログラムを作成した。さらに速度を上げるために、DIAMOND(Buchfink et al、2015)アライナーをBLAST +(Camacho et al、2009)とともにオプションとして実施した。 DIAMONDは、近年開発された配列アライナーで、BLAST +パッケージのBLASTXよりも20,000倍も高速であり、同等近い感度を備えている。Enterococcus faeciumおよびEnterococcus faecalis由来の2セットのWGSデータを用いてFast-GePをテストした。

 

GeP


インストール

ubuntu14.04のminiconda2.4.0.5環境でテストした(docker使用。ホストOS: mac os10.12)。

依存

Tsted in MacOS (OS X 10.9.5 and High Sierra v 10.13.1) and Ubuntu (14.04 and 16.04)

  • BLAST+
  • DIAMOND (v0.9.10.111 or higher)
  • MAFFT
  • SplitsTree4 (to generate a phylogenetic network from the output file Splitstree.nex)
conda install -y -c bioconda blast diamond mafft

 SplitsTree4については、HPから.dmgファイルをダウンロードし、指示に従ってインストールする。

本体 Github

git clone https://github.com/jizhang-nz/fast-GeP.git
cd fast-GeP/

perl fast-GeP.pl -h

perl fast-GeP.pl -h

 

Fast Genome Profiler (Fast-GeP) - extraction of allele profiles from genomic sequences with genome-by-genome approach.

version: 1.0

date: 4.1.2018

 

Example commands: 

fast-GeP -g list_genomes.txt -r reference.gbk

 

or,

perl fast-GeP.pl -g list_genomes.txt -r reference.gbk

 

 

Switch:

   -h: Show help

   -v: Do not show progress on the screen.

   -l: Produce extended results.

   -n: Do not produce pairwise comparison files for all the isolates. 

   -b: Using DIAMOND instead of BLAST+ as aligner in the primary search.

 

Option:

   -g: Name of the text file listing the genomic sequences need to be analyzed. Each file name should occupy a line.

   -r: Name of the reference genome sequence (GenBank format).

   -c: Minimum coverage of alignment to define an allele. [Default = 100].

   -i: Minimum identity percentage to define an allele. [Default = 80].

   -t: Minimum coverage of alignment to define a truncated allele. [Default = 50].

   -d: Number of N filled in concatinating the assemblies. [Default = 200]. 

 

テストラン

解凍

cd fast-GeP/Examples/E.faecalis/
unzip input_files.zip && cd input_files/

 

input_files(それぞれアセンブリのmuliti-fastaファイル)

f:id:kazumaxneo:20181115070809p:plain

 

ランにはリファレンスのgenebankファイル、比較したい生き物のゲノムfastaファイル (complete or draft) を使う。ゲノムfastaは場所を示したリストファイルにして与える必要がある(テストファイルでは用意されている)。(*1)

> head fast-GeP/Examples/E.faecalis/input_files/list.genomes.txt 

 $ head fast-GeP/Examples/E.faecalis/input_files/list.genomes.txt 

EN167.fas

EN182.fas

EN19.fas

EN208.fas

EN209.fas

EN215.fas

EN242.fas

EN24.fas

EN258.fas

EN28.fas

——

準備できたら実行する。

#diamondを使用
perl /fast-GeP/fast-GeP.pl -g list.genomes.txt -r EN893.gbk -b

#またはblast+を使用
perl /fast-GeP/fast-GeP.pl -g list.genomes.txt -r EN893.gbk -l

終わるまで数時間かかった。

 

出力

scheme_1541933704

f:id:kazumaxneo:20181111205135p:plain

output_1541933704

f:id:kazumaxneo:20181111205139p:plain

difference_matrix.html

f:id:kazumaxneo:20181111205105p:plain

allele_calling.txt

f:id:kazumaxneo:20181115161804j:plain

Splitstree.nexをSplitsTree4で開く。

f:id:kazumaxneo:20181111205106p:plain

 

引用

Genome-by-genome approach for fast bacterial genealogical relationship evaluation
Zhang J, Xiong Y, Rogers L, Carter GP, French N

Bioinformatics. 2018 Mar 28


Refinement of whole-genome multilocus sequence typing analysis by addressing gene paralogy
Zhang J, Halkilahti J, Hänninen ML, Rossi M

J Clin Microbiol. 2015 May;53(5):1765-7

 

参考 

*1

リストファイルを作成するなら、例えば

ls -1 *fasta > list

 

(ウィルス) コドンを考慮し、フレームシフトエラーに強いアライメントツール VIRULIGN

 

 多くのウイルス性病原体、特にRNAウイルスは、宿主内および宿主間で急速に進化しており、変化する状態への適応のマーカーがそれらのゲノムにおいて検出され得る(Lemeyら、2006)。ウイルス遺伝子型からの構造、機能および表現型予測は、ウイルス感染の薬物設計、診断および臨床管理の進歩を促進した(論文より Houldcroft et al、2017; PybusおよびRambaut、2009; Theys et al、2015)。ウイルスの遺伝子データは、進化の歴史や能動的な疫学的サーベイランスの推論にも必要不可欠である(Dellicour et al、2018; Hadfield et al、2018; Libin et al、2017)。しかしながら、遺伝子型依存性の適用は、基礎となる配列アライメントの質に強く影響される。

 ウイルス配列をアライメントさせるプロセスは、それらの広範な遺伝的多様性および頻繁な挿入および欠失によってチャレンジングになり、その結果、多数のアラインメントソフトウェアが異なる目的および用途で存在する。近年、ウイルスの集団を研究するためシーケンシングリードをマッピングおよびアセンブルするためのアライナーは、近年著しく進歩している(Posada-Cespedes et al。、2017)。ウィルスコンセンサスまたはサンガーシーケンシングリードをアライメントさせるアルゴリズムは、ペアワイズまたは複数アラインメントをもたらし、時間の経過とともに進歩が少なくなった。しかしながら、そのようなアライメントは、公衆衛生および診断の様々な側面にとって重要である。

 ウイルス遺伝子またはゲノム配列のマルチプルシーケンスアライメント(MSA)は、MAFFT、MUSCLEまたはClustal Omega(Edgar、2004; Katoh and Standley、2013; Sievers et al、2011)などの漸進的反復アプローチによって構築されることが多い。これらの発見的方法は、フレームシフトエラーを緩和する能力が低く、シーケンスデータのノイズに敏感であり得る。これは、タンパク質配列を正しいオープンリーディングフレーム(ORF)で分析する必要がある場合に有害である(例えば系統発生または薬剤耐性変異検出におけるコドン置換モデル使用時)。あるいは、リファレンス配列ガイドでのアラインメントプロセスこれらの制限を克服することができ得る(Tzou et al、2017)。しかしながら、あまりアノテーションされていないリファレンス配列の使用は、アライメントを阻害する。さらに、データセット内の劣った配列はMSA結果に大きな影響を及ぼし、自動拒絶によMSAプロセスを抑制することにより、アライメントの再現性および品質がさらに改善される。

 著書らは、closely relatedなウイルスのタンパク質コード配列のため、高速かつリファレンスガイドかつコドン訂正およびアノテーションを行うツールVIRULIGNを開発した。VIRULIGNはクロスプラットフォームGNU / LinuxUnixMacOSWindows)で使いやすいコマンドラインアプリケーションである。 VIRULIGNは、実験的に示されているように(論文のセクション5参照)、アルゴリズムの計算上の複雑さを分析して(論文のセクション4参照)、計算効率の良い方法で大規模な配列データセットを扱うことができる。単一のORFでは、VIRULIGNのアラインメントアルゴリズムは、遺伝子順序が保存されたclosely relatedなウイルスゲノムを前提に設計されているので、1塩基置換から生じるコドン異常のアライメントを訂正できる。

(以下略)

 

チュートリアル

GitHub - rega-cev/virulign-tutorial

 

f:id:kazumaxneo:20181114162751p:plain

VIRULIGN operates by aligning each target sequence (i.e., t in T) of the input file codon-correctly against the reference sequence (r). Subsequently a multiple sequence alignment MSA(r,T) is constructed based on all codon-correct (cc) pairwise aligned target sequences A_{cc}(r,t). チュートリアルより

 

VIRULIGNに関するツイート

 

インストール

mac os 10.12のanaconda2.4.3環境でテストした。

依存

Github

リリース(link)より実行ファイルをダウンロードする。

>./virulig

$ ./virulign

Usage: virulign [reference.fasta orf-description.xml] sequences.fasta

Optional parameters (first option will be the default):

  --exportKind [Mutations PairwiseAlignments GlobalAlignment PositionTable MutationTable]

  --exportAlphabet [AminoAcids Nucleotides]

  --exportWithInsertions [yes no]

  --exportReferenceSequence [no yes]

  --gapExtensionPenalty doubleValue=>3.3

  --gapOpenPenalty doubleValue=>10.0

  --maxFrameShifts intValue=>3

  --progress [no yes]

  --nt-debug directory

Output: The alignment will be printed to standard out and any progress or error messages will be printed to the standard error. This output can be redirected to files, e.g.:

   virulign ref.xml sequence.fasta > alignment.mutations 2> alignment.err

 

 

テストラン

データの準備

git clone https://github.com/rega-cev/virulign-tutorial.git
cd virulign-tutorial/examples-alignments/DENV/

デングウイルス、ジカウィルス、HIVのデータセットになっている。ここではデングウイルスのデータセットを使い、virulignを実行する。 デングウイルスは1つだけコード領域を持つ(NCBI nucletide)。チュートアリルのデータはDengue Virus Variation Database (link) [Hatcher et al, 2017]の3539のゲノムデータ(全領域)がmulti-fastaで集められている。

 

 virulignのランにはリファレンスfastaとターゲットシーケンシングのfastaを指定する。ユーザー指定の方法(--exportKind)で、リファレンスとデータセットのアライメントが行われる。例えばglobal alignmentモードで実行 (wiki)。

virulign NC_001477.fasta denv-1.fasta --exportKind GlobalAlignment > alignment
  • --exportKind [Mutations | PairwiseAlignments | GlobalAlignment | PositionTable | MutationTable]

アライメントパラメータはオプションをつけることで変更できる。また--exportAlphabet Nucleotides をつけることで、defaultではアミノ酸に翻訳して出力だが、塩基出力に変更できる。

 

またはPairwise Alignmentsモードで実行(それぞれの2配列間でアライメント)

virulign NC_001477.fasta denv-1.fasta --exportKind PairwiseAlignments > alignment

 

変異情報のテーブルを出力

virulign NC_001477.fasta denv-1.fasta --exportKind MutationTable > mutation

#またはmutations
virulign NC_001477.fasta denv-1.fasta --exportKind Mutations > mutation

 

 

XMLを使うとアノテーションができるようですが、genebankのXMLは直接扱えず、カスタムXML形式に変換する必要があるようです。そのコマンドも配布されていますが、詳細はチュートリアルで確認して下さい(GitHub - rega-cev/virulign-tutorial)。HIVxmlはreferenceに準備されています。

https://github.com/rega-cev/virulign/tree/master/references

引用

VIRULIGN: fast codon-correct alignment and annotation of viral genomes

Libin P, Deforche K, Abecasis AB, Theys K

Bioinformatics. 2018 Oct 8

参考

チュートリアルでは、別のアライメントツールとも比較して、実行時間、フレームシフトエラーに対する訂正精度を簡単に議論しています。使っているのは、本ツール、mafft、muscle、Clustal Omegaで、以下のようにコマンドを打っています。

#本ツール
virulign NC_001477.fasta denv-1.fasta \
--exportKind GlobalAlignment \
--exportAlphabet Nucleotides > denv-1-virulign.fasta

#mafftの場合
mafft --auto denv-1.fasta > denv-1-mafft.fasta

#muscleの場合
muscle -maxiters 1 -diags -in denv-1.fasta -out denv-1-muscle.fasta

#Clustal Omega の場合
clustalo --auto -i denv-1.fasta -o denv-1-clustalo.fasta

 

 

 

 

Freiburg RNA tools

 

 RNA生物学は分子生物学および生物医学研究における重要なtopicである。biological systemsにおけるRNAの機能は e.g.,  病気のプロセスに関するイノベーション(1)からCRISPR-Casに基づく最近の遺伝子編集のイノベーション(2,3)に至るまで、複雑で範囲が広い。それらの配列および他の核酸またはタンパク質との相互作用を含む、核酸の分子特性を調べる広範囲のバイオインフォマティクスツールが開発されている。

ここでは、RNA-RNA相互作用の予測、配列解析(デザイン、スプライシング、多型)を含むRNA解析用のツール群、CRISPRサイトの分類を含む単一のプラットフォームであるFreiburg RNAツールwebserver(4)の大幅なアップデートを報告する 。現在、毎月約2500のジョブが処理されている。 RNA特有のツールに加えて、ウェブサーバは、education and teachingのためのインタラクティブアルゴリズム実装を提供する。

論文の以下では、一般的なWebサーバアーキテクチャの概要を示し、続いて機能サービスとその成功したアプリケーションの概要を示す。

 
Freiburg RNA Tools

http://rna.informatik.uni-freiburg.de

f:id:kazumaxneo:20181002165507j:plain

 

 1、IntaRNA

 IntaRNAウェブサーバー(ref.4,8)は、提供されるサービスの中で最も人気がある。2RNA分子間の相互作用を調べる機能を提供している。最近、論文でベンチマークされた最先端の方法を使っており (paper link) 、クエリ/ターゲット配列間の予測に加えて、2010年の最初のリリース(ref.4)では利用できなかった原核生物のsmall RNA(sRNA)ターゲットのゲノムワイドスクリーンも提供している。

http://rna.informatik.uni-freiburg.de/IntaRNA/Input.jsp

f:id:kazumaxneo:20181002170501j:plain

ncRNA-mRNA prediction example

output

f:id:kazumaxneo:20181002183802p:plain

f:id:kazumaxneo:20181002183813p:plain

 

2、CopraRNA

CopraRNA (link)は、原核生物のsmall RNAターゲット予測のアルゴリズム(ref.7,9)(link)。偽陽性予測を有意に減少しながら系統で保存されたsRNAを特徴付ける(ref.10)。多様な細菌種を用いた多くの研究に適用されている(ref.11-18)。

http://rna.informatik.uni-freiburg.de/CopraRNA/Input.jsp

f:id:kazumaxneo:20181002170420j:plain

使うには、最低3生物のホモロガスなsRNA配列をFASTA形式で与える必要がある。 名前はNCBI RefseqのNZ_* または NC_XXXXXXになっていなければならない。

IsaR1 target predictions on 19 cyanobacterial genomes. 

CopraRNA - Results

f:id:kazumaxneo:20181003103140j:plain

  

3、GLASSgo

GLobal Automatic Small RNA Search goはsRNAのホモログを検出する(ref.20 pubmed)。反復的なBLASTn(21)検索と自動適応グラフベースのクラスタリングアルゴリズムが組み合わせられている。同定された同種のsRNA(FASTA)は、相互作用する系統樹で視覚化される。 RNAlien(ref.22)とInfernal(ref.23)の組み合わせやBLASTnのみのような関連ワークフローとは対照的に、GLASSgoは完全に自動化されており、技術的障壁が低くなっている。さらに、RNAlien + Infernalより約2桁高速である。GLASSgoを用いてsRNA OxySのホモログを徹底的に検出した研究が最近報告された(ref.24 pubmed)。 

http://rna.informatik.uni-freiburg.de/GLASSgo/Input.jsp

f:id:kazumaxneo:20181002170536j:plain

IsaR1 from PCC6803

http://rna.informatik.uni-freiburg.de/GLASSgo/Result.jsp?toolName=GLASSgo&jobID=IsaR1_PCC6803

f:id:kazumaxneo:20181003103715j:plain

  

4、metaMIR

metaMIR(ref.26)はmiRNA予測ツール。利用可能な多数のmiRNA予測ツールのデータを統合するために開発された。各アルゴリズムは、特定のmiRNAがメッセンジャーRNA(mRNA)を標的とする可能性を予測するための一連の特性に基づいている。 metaMIRは、これらの予測と、ゲノミクス、プロテオミクス、およびキュレーションされたデータベースから、検証されたmiRNAターゲッティング結果(規制の正と負の証拠)の独自に生成されたコレクションと組み合わせている。これは、陽性標的化の可能性(例えば、標的のmiRNAダウンレギュレーション)のスコアリングを可能にするだけでなく、他の標的遺伝子を標的としない選択された遺伝子を標的とするmiRNAのリストを改良するための非標的化を明示的に予測する。

http://rna.informatik.uni-freiburg.de/metaMIR/Input.jsp

f:id:kazumaxneo:20181002170605j:plain

 

RNA構造予測ツール

5、LocARNA

は、複数のRNAのmultiple alignmentsおよび2時構造予測を行い、未知の構造を有する非コードRNAの比較分析を行う。

 

6、MARNA

複数のグローバルRNAアライメントを計算する(ref.32)。 

http://rna.informatik.uni-freiburg.de/MARNA/Input.jsp

論文参照。

7、ExpaRNA

ExpaRNAは、グローバルではなくローカルに高いシーケンスアライメントを示す部位を探すことによって、RNAの構造的なモチーフを予測する(ref.35)。

http://rna.informatik.uni-freiburg.de/ExpaRNA/Input.jsp

他にも様々なツールがある。

 

education

RNA関連アルゴリズム理解をサポートするために、インタラクティブJavascriptベースの基本的なメソッドと拡張メソッドのインタフェースが提供されている(ref.48)。

http://rna.informatik.uni-freiburg.de/Teaching/

f:id:kazumaxneo:20181002171228j:plain

 

引用

Freiburg RNA tools: a central online resource for RNA-focused research and teaching
Raden M, Ali SM, Alkhnbashi OS, Busch A, Costa F, Davis JA, Eggenhofer F, Gelhausen R, Georg J, Heyne S, Hiller M, Kundu K, Kleinkauf R, Lott SC, Mohamed MM, Mattheis A1, Miladi M, Richter AS, Will S, Wolff J, Wright PR, Backofen R

Nucleic Acids Res. 2018 Jul 2;46(W1):W25-W29