macでインフォマティクス

macでインフォマティクス

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

メタゲノムデータを種レベルで検出し割合を計算するMOCAT

 

 

公式サイト

f:id:kazumaxneo:20170904104753j:plain

チュートリアル

http://vm-lux.embl.de/~kultima/MOCAT/tutorial.html

ダウンロード

http://vm-lux.embl.de/~kultima/MOCAT/download.html

fetchMGとの違い

Taxonomic profiling using mOTUs

 

インストール

macOSでは動作しない。cent OSに導入した。

公式サイトからダウンロードしたフォルダを解凍して、以下のコマンドを打つ。

cd MOCAT
./setup.MOCAT.pl
source ~/.bash_profile #パスを通す。

 setup.MOCAT.plのインストール中にexampleデータ(3.5GBと2.3GB)のダウンロードなどの選択肢が何度か出る。yes、noで判断して進める。

 

userachのインストール。

USEARCH downloadからlinuxのversion5をダウンロードする。

mv usearch5.2.32_i86linux32 usearch #リネームする。

 

 

 

ラン

MOCAT.pl -sf my.samples -rtf 
MOCAT.pl -sf my.samples -a MOCAT.pl -sf my.samples -gp assembly
MOCAT.pl -sf my.samples -make_gene_catalog -assembly_type assembly
MOCAT.pl -sf my.samples -annotate_gene_catalog
MOCAT.pl -sf my.samples -s my.samples.padded -identity 95
MOCAT.pl -sf my.samples -f my.samples.padded -identity 95
MOCAT.pl -sf my.samples -p my.samples.padded -identity 95 -mode functional MOCAT.pl -sf my.samples -ss

 

作成途中

 

 

引用
MOCAT2: a metagenomic assembly, annotation and profiling framework
Kultima JR, Coelho LP, Forslund K, Huerta-Cepas J, Li SS, Driessen M, Voigt AY, Zeller G, Sunagawa S, Bork P

Bioinformatics. 2016 Aug 15;32(16):2520-3

 
MOCAT: a metagenomics assembly and gene prediction toolkit
Kultima JR1, Sunagawa S, Li J, Chen W, Chen H, Mende DR, Arumugam M, Pan Q, Liu B, Qin J, Wang J, Bork P

PLoS One. 2012;7(10):e47656

 

 

配列のクラスタリングツール UCLUST

 2019 9/29 help追加

2019 9/30 fastaへの変換コマンド追加

 

相同な配列をクラスタリングするツール。相同性の下限値を指定してランすると、閾値以上の相同性を持った塩基配列をまとめてくれる。CD-HIT-ESTより高速に動作するとされる。

 

 

ダウンロード (linux, mac)

http://www.drive5.com/uclust/downloads1_2_22q.html

マニュアル

http://www.drive5.com/uclust/uclust_userguide_1_1_579.pdf

 

配列が一致してクラスタリングされるかどうかは、配列とクエリーとのlocal alignment(e.g., BLAST)ではなく、global alingmentの結果で判定される。

f:id:kazumaxneo:20170903180611j:plain

公式サイトから引用。

このアライメントではidentityは66.6%となる。データベースはメモリーにキャッシュされ、それからクエリとの検索が行われる。クエリと既存のクラスターが合致すれば其のクラスタ編入され(下図 右)、マッチするクラスターがなければ新しいクラスターが作られる(下図 左)。

f:id:kazumaxneo:20170903181646j:plain

公式サイトから転載。

 

インストール

ダウンロードした実行ファイルを解凍し実行権をつける。パスの通っているディレクトリに移動しておく。

mv uclust1.2.22q_i86darwin64 uclust #リネーム
chmod u+x uclust
mv uclust /user/local/bin/

uclust

]$ uclust 

uclust v1.2.22q

(C) Copyright 2009-10 Robert C. Edgar

Licensed ONLY for use in PyNAST and QIIME.

 

 

Sort sequences by length, use --mergesort for very large files:

   uclust --sort seqs.fasta --output seqs_sorted.fasta

   uclust --mergesort seqs.fasta --output seqs_sorted.fasta [--split N] [--tmpdir dir]

      --split N is partition size in Mb (default 1000).

      --tmpdir dir is directory for partition files (default current dir).

 

Cluster de novo:

   uclust --input seqs_sorted.fasta --uc results.uc --id 0.90

 

Database search without clustering:

   uclust --input seqs.fasta --lib database.fasta --uc results.uc --id 0.90 --libonly

 

Database search, cluster seqs that don't match:

   uclust --input seqs_sorted.fasta --lib database.fasta --uc results.uc --id 0.90

 

File format conversions:

   uclust --uc2clstr results.uc --output results.clstr

   uclust --clstr2uc cd_hit.clstr --output results.uc

   uclust --uc2fasta results.uc --input seqs.fasta --output results.fasta [--types abc]

      --types abc...  Record types to include (default SH).

 

Create multiple alignment for each cluster:

   uclust --staralign results.fasta --output results_aligned.fasta

 

Sort results by cluster number:

   uclust --sortuc results.uc --output results_sorted.uc

 

Input:

   --usersort             Don't assume input is sorted by length (default assume sorted).

   --maxlen L             Ignore query sequences longer than L (default 10000).

   --minlen L             Ignore query sequences shorter than L (default 16).

   --amino                Input is amino acids (default infer from ACGTU frequencies).

   --nucleo               Input is nucelotides (default infer from ACGTU frequencies).

 

Output:

   --trunclabels          Truncate FASTA labels at first white space.

   --allhits              Write all accepts to .uc file (default top hit/no match only).

   --[no]output_rejects   Write rejects to .uc file (default output top hit/no match only).

   --log filename         Log file with miscellaneous information.

   --fastapairs filename  Pair-wise alignments for all hits in FASTA format.

   --blastout filename    BLAST-like alignments.

     --rowlen n             Row length for --blastout (default 64).

     --idchar c             Character for identities in --blastout (default '|').

     --diffchar c           Character for mismatches in --blastout (default ' ').

     --[no]blast_termgaps   Show terminal gaps in --blastout (default don't show).

 

Search:

   --id f                 Minimum identity for a hit (default 0.9).

   --maxaccepts n         Maximum hits before quitting search (default 1, 0=infinity).

   --maxrejects n         Maximum rejects before quitting search (default 8, 0=infinity).

   --w n                  Word length for windex (default 5 aa.s, 8 nuc.s).

   --[no]wordcountreject  Dis/enable word count rejection (default enabled).

   --bump n               Bump percent (default 50, 0=don't bump).

   --stepwords n          Target nr. of common words (default 8, 0=don't step).

   --rev                  Reverse strand matching (default plus strand only)

   --libonly              Match to --lib only (default add new seed if not matched).

   --self                 Ignore target sequence with same label as query.

   --idprefix n           First n letters of query and seed must be identical (default 0).

   --exact                Same as --maxrejects 0 --nowordcountreject.

   --optimal              Same as --maxrejects 0 --maxaccepts 0 --nowordcountreject.

 

Alignment:

   --match s              Match score (nucleotides only, default 2.0).

   --mismatch s           Mismatch score (nucleotides only, default -1.0).

   --gaopen s             Gap open penalties (see manual).

   --gapext s             Gap extend penalies (see manual).

   --[no]fastalign        [Don't] use HSPs and banding for speed (default do use).

   --hsp n                Minimum length of HSP (default 32).

   --hspscore s           Minimum score/col for HSP (default 1.0)

   --k n                  Word length for HSP-finding (default 3 amino acids, 4 nucleotides).

   --band n               Band radius for regions between HSPs (default 16, 0=don't band).

   --check_fast           Compare fast & slow alignments and report statistics in --log file.

 

Misc:

   --help                 This help.

   --quiet                Don't output progress messages to stderr.

   --version              Write version to standard output and exit.

 

 

実行方法

 Clustering

uclust --sort seqs.fasta --output seqs_sorted.fasta
uclust --input seqs_sorted.fasta --uc results.uc --id 0.90
  • --id Minimum identity for a hit (default 0.9). 

fastaへの変換

uclust --uc2fasta results.uc --input seqs_sorted.fasta --output result.fasta 

 

非冗長(nr)なデータベースを作るには以下のコマンドを打つ。

uclust --input seqs_sorted.fasta --uc2fasta results.uc --types S --output nr.fasta

 

 Database search

既存のデータベースを使ってクラスタリングを行う。

uclust --sort seqs.fasta --output seqs_sorted.fasta
uclust --input seqs_sorted.fasta --lib database.fasta --uc results.uc --id 0.90

 

Multiple alignment

uclust --input seqs_sorted.fasta --uc results.uc --id 0.90 
uclust --uc2fasta results.uc --input seqs_sorted.fasta --output results.fasta
uclust --staralign results.fasta --output aligned.fasta

 

引用

Search and clustering orders of magnitude faster than BLAST

Robert C. Edgar

Bioinformatics, Volume 26, Issue 19, 1 October 2010, Pages 2460–2461,  

メタゲノムデータを種レベルで検出し割合を計算するmOTUとfetch-MG

追記9/5;ソフト名や使い方を勘違いしておりましたので修正します。

 

 環境サンプル中の種の多様性を評価する手法として16S rRNA遺伝子を特異的に増幅する手法がよく知られているが、種によっては配列の異なるrRNA遺伝子を複数持つことがある。ここにPCR増幅のbiasもかかってくるため、16S rRNAだけでメタゲノムデータを評価すると、特に近縁種間の誤差が大きくなる可能性がある。

最近ではNGSを使った環境サンプルの全ゲノムシーケンスが選択肢として選べるようになった。しかしながら、メタゲノムシーケンスから菌叢の定量を行うと、遺伝子のコピー数の違いからbiasが非常に大きくなる。そのため、高精度な定量を行うには生物全てが保有しているようなユニバーサルでシングルコピーな(つまり滅多に水平伝達で複製しない)遺伝子セットを使ったメタゲノムデータの定量が必要になってくる。

これまでの研究で、シングルコピーでユニバーサルな遺伝子セットが見出されてきた(Ciccarelli et al., Science, 2006; Sorek et al., Science, 2007; von Mering et al., Science, 2007)。fetch-MGは先行研究で見出されたシングルコピーでユニバーサルな40遺伝子を使い、メタゲノムデータから菌の種類と割合を計算・出力する。

 

 

公式サイト

http://www.bork.embl.de/software/mOTU/index.html

チュートリアル

http://www.bork.embl.de/software/mOTU/tutorial.motu.standalone.html

 

インストール

mOTUとfetch-MGのダウンロードリンクもチュートリアルに存在している。

 fetch-MG自体はmacで動作するが、bin/の実行ファイルがmacで動作しなかったのでcent OSにインストールした。

 

追記

condaを使う

#bioconda (link)
conda install -c bioconda -y motus

 

 

実行方法

fetch-MGを使い40のMGs(protein-coding phylogenetic marker genes)を抽出する。

./fetchMG.pl -m extraction -x bin example_dataset/example_data.faa
  • -t Number of processors/threads to be used
  • -o Output directory; default = 'output'
  • -h Path to directory that contains hmm models; default = './lib'
  • -p Set if nucleotide sequences for filename.faa is not available
  • -d Fasta file with DNA sequences of the same genes; not neccesary if protein file and dna file have the same with .faa and .fna suffixes

入力はfaaファイル(ただし40のMGsが抽出できる前提で動作しているので、ゲノムの全タンパク質配列が揃っていないと正しく動作しない)。またdna配列(同じ名前で与える必要がある。e.g., example.fnaexample.faa)があれば、該当するタンパク質の塩基配列も出力してくれる。詳細はexample_datasets/を参照。

 

 

作成中。

ここから下は下書きです。

 

 

ペアードエンドfastqを指定してラン(非圧縮かgz圧縮のfastqに対応)。

./motus.pl sample.1.fq sample.2.fq
  • --verbose Be more verbose while running the analysis
  • --processors=N This should be an integer and defines the number of processors that the script will use.
  • --length-cutoff=L The minimum size per read (after quality-based trimming), default: 45.
  • --identity-cutoff=I Minimum percentage identity in alignment (default: 97) --quality-cutoff=Q Basepair quality cutoff (default: 20)
  • --fastq-format The format of the input files. Must be one of 'auto' (the default), 'sanger', or 'illumina'. Note that new Illumina machines actually use the 'sanger' format. The auto-detection should generally work well.
  • --output-directory Where to place the final results file (by default it uses a directory named ``RESULTS``).

 

RESULULTSに解析結果が出力される。 検出された菌の種名と割合が出力されている。

f:id:kazumaxneo:20170903161302j:plain

  

 

 

 複数fastqを解析する場合、fastqをそれぞれのディレクトリに入れて、其のパスを記載したファイル(1行ずつ記載)を指定してランする。

./mOTUs.pl --sample-file input.txt

 

 

 

引用

Metagenomic species profiling using universal phylogenetic marker genes.

Sunagawa S1, Mende DR, Zeller G, Izquierdo-Carrasco F, Berger SA, Kultima JR, Coelho LP, Arumugam M, Tap J, Nielsen HB, Rasmussen S, Brunak S, Pedersen O, Guarner F, de Vos WM, Wang J, Li J, Doré J, Ehrlich SD, Stamatakis A, Bork P.

Nat Methods. 2013 Dec;10(12):1196-9. doi: 10.1038/nmeth.2693. Epub 2013 Oct 20.

 

 

ユニバーサルな遺伝子に関する 先行研究

 

関連


メタゲノムデータをbinningして出力可能なGUIアプリ VizBin

2019 7/5文章修正

 

 

VizBinはメタゲノムデータをレファレンスに依存せずにbinnigする手法。テトラヌクレオタイド頻度情報を使いアセンブルデータを分類する。最終的に2次元のPCAプロットとしてビジュアル化してくれる。どこからどこまでを1つの生物として抽出するかは、人間が目で見て決定する。論文にはThe human eye-brain system for cluster identificationと書かれている。複雑なメタゲノムでは使いにくいが、比較的な単純な(low complexitiesな)メタゲノムデータであれば目で見て分類は可能であり、面白いツールである。

アプリはJAVAで構築されている。

 

公式ページ

Download appから本体をダウンロード可能。

f:id:kazumaxneo:20170903124624j:plain

Tutorial

Tutorial · claczny/VizBin Wiki · GitHub

テストデータ

example dataset AMFJ01

 

実行方法

ビルド済みのjava実行ファイルが配布されている。これを叩いて起動する。

> java -jar VizBin-dist.jar

 

File to visualiseにcontigファイルを入力する。あらかじめメタゲノム用のアセンブラでcontigは作っておく必要がある。

f:id:kazumaxneo:20170902204406j:plain

計算用のthreadsを増やして計算を高速化することもできる。

 

startをクリック。抽出されたkmerの数が表示される。

f:id:kazumaxneo:20170902204541j:plain

 ↓

しばらく待つ。 

 ↓

結果が表示された。

f:id:kazumaxneo:20170902205020j:plain

 ↓

trackpad(2本指で上下にドラッグ)やマウスのホイールで拡大縮小が可能。

 ↓

クリックして領域を囲むと、囲った領域のプロットを全て含むfastaが出力可能。

f:id:kazumaxneo:20170902205302j:plain

  ↓

exportを選択。

f:id:kazumaxneo:20170903122240j:plain

  ↓

fataが出力される。選択を間違えた時は右クリック->Selection->clear selectionでやり直しできる。

 

 

contigのサイズ、GCなども読み込める。exampleは以下のようになっている。

f:id:kazumaxneo:20170903123909j:plain

最初にヘッダー行が必要である。coverage、length、label、isMarkerなどの名前が認識に必要。contigの順番はfastaの順番と合致している必要がある。

 

入力する時はadvanxed modeに切り替える。show advanxed optionsをクリック。

f:id:kazumaxneo:20170903130235j:plain

  ↓ 

aanotton fileの項目をクリック。

f:id:kazumaxneo:20170903130408j:plain

他にもいくつかのパラメータを変更可能である。詳細は公式マニュアルを参照(一番下に説明があります)。

 

引用

VizBin - an application for reference-independent visualization and human-augmented binning of metagenomic data

Cedric C Laczny, Tomasz Sternal, Valentin Plugaru, Piotr Gawron, Arash Atashpendar, Houry Hera Margossian, Sergio Coronado, Laurens van der Maaten, Nikos Vlassis, Paul Wilmes

Microbiome 2015 3:1

 

複数のcontigをマージしアセンブリの連続性を改善する Mix

2019 6/11 追記

 

ゲノムアセンブリ構築の利点を得ることを妨げる課題の中には、未完成のアセンブリおよびその後の実験的な費用の両方がある。第一に、ゲノムデノボアセンブリのための多数のソフトウェアソリューションが利用可能であり、それぞれがその長所と短所を有し、それらの中からどのように選択するかに関する明確な指針がない。次に、これらのソリューションではドラフトアセンブリが作成され、多くの場合リソース集約型の仕上げフェーズが必要になる。

 本論文では、リファレンスゲノムに頼らずにコンティグ断片化を減らし、ゲノムの仕上げをスピードアップするという目標を持たずに、2つ以上のドラフトアセンブリをミックスするツールであるMixを開発することによってこれら2つの側面に対処する。提案したアルゴリズムは、頂点が連続体の末端を表し、辺がこれらの末端間の既存のアライメントを表す拡張グラフを構築する。これらのアライメントエッジはコンティグエクステンションに使用される。結果の出力アセンブリは、累積コンティグ長を最大化する拡張グラフ内の一連のパスに対応する。

 GAGE-B研究からの細菌NGSデータに対するMixの性能を評価し、そしてそれを新しくシーケンシングされたマイコプラズマゲノムに適用する。結果として得られる最終アセンブリは、全体的なアセンブリ品質の大幅な向上を示していた。特にMixは、de novoプロジェクトの場合のように、標準的なアセンブリ統計のみによって選択が導かれる場合でも、より良い全体的な品質結果を提供することで一貫していた。

 

 

 

 Mixはバクテリア向けに設計された、複数のconitgをマージしてより長いcontigを作る方法論。うまく使えば、細分化されたcontigをからより長いcontigを作ることができる。

  

インストール

依存

MuMmerはbrewで、Biopythonとnetworkxはpipでそれぞれ導入できる(注 nucmer v4ではエラーを起こした)。

本体 Github

 解凍してbin/のMix.pyにパスを通しておく。 

 

実行方法

3つのcontigを合体する。

preprocessing.py -o contigs.fa Assembly1.fa Assembly2.fa Assembly3.fa

#catでも動く
cat Assembly1.fa Assembly2.fa Assembly3.fa > merged.fasta

3つの contig.faがコンカテネートして1つのfastaファイル(contigs.fa)ができる。

 

あとは以下の3つのコマンドを実行する。

nucmer --maxmatch -c 30 -l 30 -banded -prefix=alignments contigs.fa contigs.fa
show-coords -rcl alignments.delta > alignments.coords
Mix.py -a alignments.coords -c contigs.fa -o output_dir/ -C 300 -A 200
Mixのflag一覧を貼っておく。



Options:

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

  -a ALN, --aln=ALN     the file containing the alignments (.coords)

  -o OUT, --out=OUT     the output directory where the scaffolds will be

                        written (it must already exist)

  -c CTG, --ctg=CTG     the multi FASTA file containing all the contigs that

                        were used in the alignment

  -A ATH, --ath=ATH     minimum length of alignment (default:0)

  -C CTH, --cth=CTH     minimum length of contig (default:500)

  -d, --dot             write the graphs in dot format

  -g, --graph           write the graphs in cytoscape format

  -r, --restrict-to-aligned

                        If on, restrict output to aligned coords

 

 output_dir/Mix_results_A200_C300/にいくつかのファイルができる。

user$ ls -lh output_dir/Mix_results_A200_C300/*

-rw-r--r--  1 user  staff   2.7M  9  1 22:41 output_dir/Mix_results_A200_C300/Mix_assembly.fasta

-rw-r--r--  1 user  staff   1.0M  9  1 22:41 output_dir/Mix_results_A200_C300/Mix_assembly_ctg.fasta

-rw-r--r--  1 user  staff   1.7M  9  1 22:41 output_dir/Mix_results_A200_C300/Mix_assembly_path.fasta

-rw-r--r--  1 user  staff     0B  9  1 22:41 output_dir/Mix_results_A200_C300/Mix_selected_isolated_contigs.csv

-rw-r--r--@ 1 user  staff   433K  9  1 22:41 output_dir/Mix_results_A200_C300/all_alignments.csv

-rw-r--r--@ 1 user  staff    28K  9  1 22:41 output_dir/Mix_results_A200_C300/all_contigs.csv

-rw-r--r--  1 user  staff    22K  9  1 22:41 output_dir/Mix_results_A200_C300/initial_assembly_graph.gml

-rw-r--r--  1 user  staff    26K  9  1 22:41 output_dir/Mix_results_A200_C300/reduced_assembly_graph.gml

 Mix_assembly.fastaがmixされて長くなったcontig.faである。

 

 

テスト

edena、SOAPdenovo2、Abyssで作ったcontigをMixにかけてみる。

preprocessing.py -o contigs.fa edena.fa abyss_k77.fa SOAPdenovo_k51.fasta
nucmer --maxmatch -c 30 -l 30 -banded -prefix=alignments contigs.fa contigs.fa
show-coords -rcl alignments.delta > alignments.coords
Mix.py -a alignments.coords -c contigs.fa -o output_dir/ -C 300 -A 200

quastで評価する。 

quast edena.fa abyss_k77.fa SOAPdenovo_k51.fasta Mix_assembly.fasta -R bacteria.fa -o results -t 20 

紫がMix.fa、赤がAbyss(k=77)

f:id:kazumaxneo:20170901230632j:plain

一番下がMix.fa。

f:id:kazumaxneo:20170901230652j:plain

contigの断片化が大きく解消されていることがわかる。一方で、図をよく見るとよくわからない部位で切れてもいる。原因は分からない。

 

 

ただし、アセンブルが止まるのはリピートだったり何かしら理由があるのが普通です。こういったツールを使うとアーティファクトができてくる恐れもあるので注意して使ってください。 

 

quastは以下で紹介しています。

 

引用

Finishing bacterial genome assemblies with Mix

Hayssam Soueidan, Florence Maurier, Alexis Groppi, Pascal Sirand-Pugnet, Florence Tardy, Christine Citti, Virginie Dupuy and Macha Nikolski

BMC Bioinformatics 2013 14(Suppl 15):S16

 

関連


de novoアセンブルしてバリアントをコールするDISCOVAR

 

DISCOVARは2014年にNature geneticsに載ったバリアントを検出する方法論。シーケンスデータをアセンブルして、バリアントをコールする。ヒトゲノムの構造変化は90%ほどは既存のツールで検出可能だが、残りの構造変化(low-complexity sequenceやsegmental duplications)の検出は困難とされる。DISCOVARはそうした難しい構造変化を検出するために構築されたと書かれている。ただし現在のバージョンではトリッキーな方法を使わないと50-Mb以上のゲノムをアセンブルすることはできない。どちらかといえばスモールゲノム向けのツールになっている。

バリアント検出の過程でアセンブルを行うため、アセンブルに用いることも可能である(DISCOVAR de novo)。

 

公式サイト

https://software.broadinstitute.org/software/discovar/blog/?page_id=98

マニュアル

https://docs.google.com/document/d/1U_o-Z0dJ0QKiJn86AV2o_YHiFzUtW9c57eh3tYjkINc/edit

 入力はilluminaでシーケンスした250-bp以上でインサートサイズが450くらいのペアリードfastq。サンプルはPCR freeのプロトコルで調整されたデータが望ましいとされる。カバレッジも条件があり、60前後がベストと書かれている。条件を箇条書きしておく。

  • Illumina MiSeq or HiSeq 2500 genome sequencers
  • PCR-free library preparation
  • 250 base paired end reads (or longer)
  • ~450 base pair fragment size
  • ~60x coverage

アセンブルできるサイズは50-Mbまでに制限されている。

 

インストール

依存

全てbrewで導入できる。

 brewでインストールする。macでもbrewでインストールできるが、原因不明の理由で失敗する。ソースからビルドするためのマニュアルもリンクが消えていたので、cent OSにbrewを使いインストールした。

brew install DISCOVAR

 

実行方法

1、リファンレンスを準備。

PrepareDiscovarGenome REF=reference.fasta

faiファイルといくつか関連するファイルができる。

 

2、bamにマッピングされたペアリードを抽出し、de novoアセンブルを実行するその後バリアントをコールする。

Discovar READS=bam OUT_HEAD=asssembly REGIONS=all REFERENCE=reference.fa TMP=temp 

コンマセパレートで複数bamを入力できる。READS=filename1,filename2,...

 

アセンブルがが終わるとassembly~がいくつかできる。(アセンブルできるサイズは50-Mbまで)。REGIONS=all にするとかなりのメモリが要求される。カバレッジx100のゲノム(4メガ)でランすると、peak memory uisageは114GBに達した。

領域を指定するにはREGIONS=1:50000-150000などと記載する(chr1の50000-150000の領域を探索)。regionは複数指定できるようになっている(e.g., REGIONS=chr:start-end,chr:start-end,...)。

 

3、アセンブルのgraphを可視化する(オプション)

dot -Tps -o assembly.final.ps assembly.final.dot -v
gv assembly.final.ps

 

 公式マニュアルでは、大きなゲノムを解析するために、少しづつ領域をオーバーラップしながらランする方法が書かれています。興味がある人は確認してみてください。アセンブルのgraphについての勉強にもなると思います。

de novo アセンブルによる構造変化検出ツールとして2015年にpublishされたfermikitなどがあります。

 

引用

Comprehensive variation discovery in single human genomes

Neil I Weisenfeld, Shuangye Yin, Ted Sharpe, Bayo Lau, Ryan Hegarty, Laurie Holmes, Brian Sogoloff, Diana Tabbaa, Louise Williams, Carsten Russ, Chad Nusbaum, Eric S Lander, Iain MacCallum & David B JaffeNeil I Weisenfeld, Shuangye Yin, Ted Sharpe, Bayo Lau, Ryan Hegarty, Laurie Holmes, Brian Sogoloff, Diana Tabbaa, Louise Williams, Carsten Russ, Chad Nusbaum, Eric S Lander, Iain MacCallum & David B Jaffe

Nature Genetics 46, 1350–1355 (2014) doi:10.1038/ng.3121

 

バクテリアやアーキアの遺伝子を予測するProdigal(メタゲノムデータセットにも対応)

2019 5/8 インストール追記

2021 7/13 help更新

2022/07/20 追記

2023/08/23 追記

 

ProdigalはDynamic Programmingの方法論により効率的にバクテリアアーキアの遺伝子を探すツール。既存の方法は様々存在するが、本手法はまずインプットゲノムを分析してモデルを構築し、それから遺伝子を予測することで、false positiveを抑えtrue callを増やすよう工夫されている。

Prodicalはドラフトゲノムやメタゲノムを扱うことが可能である。動作が非常に高速なのも特徴で、E.coliゲノムだと10秒程度で解析することが可能である。2010年に論文が発表され、2014年時点で600以上引用されている。

本手法はイントロンを含む遺伝子の予測やアノテーションは行われない。また、挿入や欠損によるフレームシフトは考慮されないので、高精度な予測を行うためには、予測に使うゲノムのエラーは入念に取り除かれていなければならない。ウィルスでの動作は検証されていないが、動作はするとされる。最新版はProdigal2で、今後Prodigal3へのアップデートが予定されている。

 

 

公式サイト

http://prodigal.ornl.gov

マニュアル

https://github.com/hyattpd/prodigal/wiki

 

インストール

Github

brewでProdigal2が導入できる。condaでも導入できる。

brew install Prodigal

conda install -c bioconda -y prodigal

prodigal -h

$ prodigal -h

 

Usage:  prodigal [-a trans_file] [-c] [-d nuc_file] [-f output_type]

                 [-g tr_table] [-h] [-i input_file] [-m] [-n] [-o output_file]

                 [-p mode] [-q] [-s start_file] [-t training_file] [-v]

 

         -a:  Write protein translations to the selected file.

         -c:  Closed ends.  Do not allow genes to run off edges.

         -d:  Write nucleotide sequences of genes to the selected file.

         -f:  Select output format (gbk, gff, or sco).  Default is gbk.

         -g:  Specify a translation table to use (default 11).

         -h:  Print help menu and exit.

         -i:  Specify FASTA/Genbank input file (default reads from stdin).

         -m:  Treat runs of N as masked sequence; don't build genes across them.

         -n:  Bypass Shine-Dalgarno trainer and force a full motif scan.

         -o:  Specify output file (default writes to stdout).

         -p:  Select procedure (single or meta).  Default is single.

         -q:  Run quietly (suppress normal stderr output).

         -s:  Write all potential genes (with scores) to the selected file.

         -t:  Write a training file (if none exists); otherwise, read and use

              the specified training file.

         -v:  Print version number and exit.

 

 

ラン

ランには3つのモードがある。詳細は公式マニュアルの説明を確認する。

 

フィニッシュしたゲノム(一本になっておりNがないゲノム)や巨大ウィルスはノーマルモードでランする。ノーマルモードはまずゲノムを分析し、それからコード領域を予測する。ハイクオリティなモデル構築のため、ゲノムは100-kb (3') ~ 500-kb (5') の長さ以上である事が推奨されている。

ノーマルモード

prodigal -i my.genome.fna -o gene.coords.gbk -a protein.faa -d nucleotide.fna
  • -i Specify FASTA/Genbank input file (fasta推奨).
  • -o Specify output file (default writes to stdout).
  • -a Write protein translations to the selected file.
  • -d Specify nucleotide sequences file.
  • -c Closed ends. Do not allow partial genes at edges of sequence.
  • -f Specify output format.

    gbk: Genbank-like format (Default)

    gff: GFF format

    sqn: Sequin feature table format

    sco: Simple coordinate output

出力はfaa、fna、genbank like、gtfなどから選択する。 genbankライクな出力は以下のようになっている。

f:id:kazumaxneo:20170901133923j:plain

アセンブルして作ったcontigの状態でもノーマルモードを利用できる。ノーマルモードでは、1本でも長いcontigがあれば、そこからモデルを構築し遺伝子を予測することが可能である。ただし、あまりにもバラバラだとエラーが出る。その場合は、近縁種の情報からモデルを構築してランできるトレーニングモードを使う。

 

 

ノーマルモードは同じプロファイルを入力に対して当てはめるので、メタゲノムのデータに対してそのまま使ってはならない。通常、メタゲノムアセンブルした全配列データはアノニマスモードでランする。アノニマスモードではプリセットのデータを使い解析される。100-kb以下のプラスミド、ファージ、ウィルスもアノニマスモード使用が推奨される。

 アノニマスモード(Prodigal2ではmetagenome)。

prodigal -i metagenome.fna -o coords.gbk -a proteins.faa -p meta
  • -p Select procedure (single or meta). Default is single.  

ノーマルモードでランすれば自分の配列で訓練するので精度が上がる。ノーマルモードを使いたいなら、contigをbinningしてMAGs配列を個々に指定する。アノニマスモードはsmall plasmidやlow qualityなゲノムにも使えるとされる。

追記

20-kbなど短い配列もアノニマスモードで実行するとノーマルモードより精度が上が流傾向が見られる。特にどこを開始コドンをするかなどは影響を受けやすい。アノニマスモードを使いたくないなら、20-kbなど短い配列に元のゲノム配列をcatなどで結合してノーマルモードで実行し、十分な訓練が行えるように工夫する(あとでアノテーションから分離する。)。このtipsは、prokkaやDFASTでも同様に使える。

 

 

フィニッシュした近縁種のゲノム情報が利用できるなら、ドラフトゲノムの予測はトレーニングモードが利用できる。トレーニングモードは、近縁種のゲノムでトレーニングデータを作成し、それを元にゲノム予測を行う方法である。ドラフトゲノムのクオリティが不十分な時に使う。

レーニングモード

1、ゲノム1を使い、トレーニングデータを出力する(-p trainは認識しないので使わない。よく似たゲノムを使わないと性能が出ない)。

prodigal -i genome1.fna -t genome1.trn 
  • -t Write a training file (if none exists); otherwise, read and use the specified training file.

2、ゲノム1のトレーニングデータを使い、ゲノム2の遺伝子を予測。

prodigal -i genome2.fna -t genome1.trn -o genome2.gbk -a genome2.faa

 

大抵のバクテリアアーキアはgenetic code11を使うが、code11で遺伝子が短すぎる場合、自動で他のcodeに切り替えるように設計されている(詳細)。

  • -g Specify a translation table to use (default 11).

-mをつけると、Nがある領域の遺伝子が予測されなくなる。

  • -m: Treat runs of N as masked sequence; don't build genes across them.

 

その他

GNU parallelsで並列化。singleモード。

mkdir gbk #.gff保存dir
mkdir Protein #.faa保存dir
#40 parallel (問題なければ--dry-runを外す)
ls *.fna | parallel -j 40 --dry-run 'prodigal -i {} -o {.}.gbk -a {.}.faa -p single && mv {.}.gbk GFF/ && mv {.}.faa Protein/'

 

Pyrodigal : Cython bindings and Python interface to Prodigal

https://pyrodigal.readthedocs.io/en/stable/index.html

 

公式マニュアルには使い方やターゲットについて丁寧に説明されています。ぜひ一度確認してみてください(リンク)。

引用 

Prodigal: prokaryotic gene recognition and translation initiation site identification

Doug HyattEmail author, Gwo-Liang Chen, Philip F LoCascio, Miriam L Land, Frank W Larimer and Loren J Hauser

BMC Bioinformatics 2010 11:119 

 

関連