macでインフォマティクス

macでインフォマティクス

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

メタゲノムのgene-targeted assembler: MegaGTA

 

 次世代シーケンシングは、近年のメタゲノミクスの研究を大きく促進してきた。これらの研究は、しばしば何百万から数十億のリードをde novoでアセンブリし、コンティグにして遺伝子アノテーションすることを含む。これは、メタゲノムのアセンブリ効率を大幅に向上させる高度なアルゴリズムの研究の開始のトリガーとなった[論文より ref.1,2,3]。一方、不均一なカバレッジと cross-genome repeatsのために、断片化された遺伝子配列になるのが一般的である。これらの欠点を克服するために、EMIRGE [ref.5]、REAGO [ref.6]、SAT-Assembler [ref.7]、Xander [ref.8]を含むいくつかの gene targeted assembly methodsが公開されている。 Xanderは、アセンブリ前に標的とされた遺伝子に由来する可能性があるリードを分類しようとする最初の3つの方法とは異なり、はじめにde Bruijn graph(DBG)を構築し、オンザフライでk-merのトラバーサルを直接探して遺伝子を探し出す。特定の遺伝子の集合は、基礎となる遺伝子ファミリーの複数のシーケンスアライメント結果を用いて構築された隠れマルコフモデル(HMM)のプロファイル[ref.9]によって導かれる。ノードから出発して、Xanderはgraphのブランチ部分で決定を行い、HMMの可能性が最も高いユニークパスを出力する。これは、ブランチで本来停止するデノボアセンブリの問題を克服する。より具体的には、Xanderアルゴリズムは、DBGとHMMとの組み合わせgraph上で操作される。そのワークフローは論文図1に示されており、次のセクションで詳しく説明する。

 Xanderは、以前の方法よりも長くよりハイクオリティの遺伝子特異的コンティグを産生するが、まだ改善の余地がある。このペーパーでは、以下の3つの点を改善した。

  •  A. 複数のサイズのk-merを使用する。 Xanderは固定kを使ってk-mersのde Bruijn graphを作成する。固定k-merは古典的なジレンマを引き起こす。小さいkは短いリピートで崩壊してde Bruijn graphに過度の分岐をもたらし、大きなkは低いカバレッジのゲノム領域の間にギャップをもたらすため、これらの領域由来の遺伝子はアセンブリされにくい[ref.10]。 HMM-guided assemblyは、HMMによって「示唆」された最良の経路を選択することによってリピートを解決することを目標としているが、異なる遺伝子の2つの部分が反復を介してキメラコンティグに結合する可能性はゼロではない。これに関して、Iterativeなde Bruijn graph [ref.10]は、複数のk-merでレバレッジをかけるものであり、いくつかのde novoアセンブラでその利点が示されている[3,11,12,13]。以下のベンチマークは、HMM誘導アセンブリがIterative(繰り返し)なde Bruijn graphから利益を得ることができることを示している。
  • B. De Bruijn graphの誤ったk-merをフィルターする。与えられたリードセットで一度だけ出現するk-mersはエラーの傾向が強い。低カバレッジな遺伝子の感度を上がるために、Xanderは、頻度の低いk-merをフィルタリングしない方針を選択した。代わりに、Xanderは誤ったk-mersを避けるためHMMに依存している。しかし、それは構造的な誤りまたは誤った塩基を伴う多くのコンティグをもたらし、特に小さいサイズのk-merが使用される場合に起こりやすくなる。この不具合を回避するために、著者らは、HMM誘導検索の間に1回だけ発生するk-mersにペナルティを課している(これらのk-mersに事前の誤った確率を設定することに等しい)。
  • C. 確率的なデータ構造によって引き起こされる偽陽性のk-merを避ける。より良いメモリ効率を達成するために、XanderはBloomフィルタ[ref.14]を使用してde Bruijn graphを表している。これは確率的なデータ構造であり、一定割合の偽陽性(ただし偽陰性はない)を引き起こす。著者らは、Bloomフィルタを、偽陽性偽陰性の両方から解放された最新のメモリ効率の高いデータ構造である簡潔なde Bruijn graph[ref.3、15]に置き換えて解決した。新しい手法との比較のため、Xanderのk-mer選択による誤ったコンティグを調べ、De Bruijn graphの正確な表現に問い合わせベンチマークした。

前述の3つの改善が実装され、MegaGTA(Xanderとの違いやMegaGTAワークフローも図1に示されている)と呼ばれる新しい遺伝子ターゲットアセンブラに統合された。最初に、HMP定義のmock communityと、大規模な根圏土壌メタゲノムサンプルの2つのデータセットを用いて、MegaGTAの有効性を実証した[ref.8]。Iterativeなde Bruijn graphで、MegaGTAはXander(よく調整されたパラメータ)よりも高い感度と精度を達成した。さらに興味深いことに、24コアサーバーでテストした場合、MegaGTAは、複数のk-merサイズを使うことで計算量が増えても、同様の量のピークメモリを使用しながら、Xanderより約2〜4倍高速だった。 MegaGTAに似た方法で複数のk-merサイズでXanderを繰り返し実行すると、MegaGTAの相対的なスピードアップはさらに大きくなる。

 根圏土壌データセットに関しては、Xanderがそれぞれ256GB、128GBおよび64GBのBloomフィルターサイズを使用するとき、Xanderが生成するコンティグの0.02、0.39、および10.52%が偽陽性k-mersを含むことがわかった。明らかにBloomフィルタは、メモリが制約されているときに精度の問題を引き起こす。 簡潔な de Bruijn graphは、Bloomフィルタの不正確さを克服するだけでなく、より速いgraph構築を可能にする。

 

 

f:id:kazumaxneo:20180610122012j:plain

ワークフロー。論文より転載。

 

インストール

cent OS6のpython2.7.14でテストした。

依存

UCHIMEはバイナリをHP(リンク)からダウンロードし、/usr/local/binなどに移動する。

本体 Github

https://github.com/HKU-BAL/megagta

git clone https://github.com/HKU-BAL/megagta.git 
cd megagta
git submodule update --init --recursive
./make.sh
cd bin/

> python megagta.py -h

$ python megagta.py -h

MegaGTA v0.1-alpha

 

Copyright (c) The University of Hong Kong

 

Usage:

  megagta.py [options] {-1 <pe1> -2 <pe2> | --12 <pe12> | -r <se>} -g gene_list.txt [-o <out_dir>]

 

  Input options that can be specified for multiple times (supporting plain text and gz/bz2 extensions)

    -1                       <pe1>          comma-separated list of fasta/q paired-end #1 files, paired with files in <pe2>

    -2                       <pe2>          comma-separated list of fasta/q paired-end #2 files, paired with files in <pe1>

    --12                     <pe12>         comma-separated list of interleaved fasta/q paired-end files

    -r/--read                <se>           comma-separated list of fasta/q single-end files

    

    -g/--gene-list           <string>       gene list

 

Optional Arguments:

  Basic assembly options:

    -c/--min-count           <int>          minimum multiplicity for filtering k-mers [1]

    -k/--k-list              <int,int,..>   comma-separated list of kmer size (in range 15-127)

                                            the last k must be a multiple of 3) [30,36,45]

    -p/--prune-len           <int>          prune the search if the score does not improve after <int> steps [20]

    -l/--low-cov-penalty     <float>        penalty for coverage one edges (in [0,1]) [0.5]

    --max-tip-len            <int>          max tip length [150]

    --no-mercy                              do not add mercy kmers

 

  Hardware options:

    -m/--memory              <float>        max memory in byte to be used in SdBG construction [0.9]

                                            (if set between 0-1, fraction of the machine's total memory)

    --mem-flag               <int>          SdBG builder memory mode [1]

                                            0: minimum; 1: moderate; others: use all memory specified by '-m/--memory'.

    --use-gpu                               use GPU

    --gpu-mem                <float>        GPU memory in byte to be used. Default: auto detect to use up all free GPU memory [0]

    -t/--num-cpu-threads     <int>          number of CPU threads, at least 2. Default: auto detect to use all CPU threads [auto]

 

  Output options:

    -o/--out-dir             <string>       output directory [./megagta_out]

    --min-contig-len         <int>          minimum length of contigs to output [450]

    --keep-tmp-files                        keep all temporary files

 

Other Arguments:

    --continue                              continue a MEGAHIT run from its last available check point.

                                            please set the output directory correctly when using this option.

    -h/--help                               print the usage message

    -v/--version                            print version

    --verbose                               verbose mode

 >./post_proc.sh

$ ./post_proc.sh

Usage: post_proc.sh -g <gene_resources> -d <workdir> -m <max_jvm_heap> -c <dist_cutoff> [-t <num_threads=1>] [-f (turn on framebot)]

 

 

ラン

ランにはXanderと同様にアセンブリターゲット遺伝子のHMMモデルが必要になる。Githubに掲載されている例を貼っておく。1行1遺伝子で、先頭カラムで遺伝子名を指定する。これは出力ディレクトリ名などに使われる。

rplB share/RDPTools/Xander_assembler/gene_resource/rplB/for_enone.hmm share/RDPTools/Xander_assembler/gene_resource/rplB/rev_enone.hmm share/RDPTools/Xander_assembler/gene_resource/rplB/ref_aligned.faa
nirK share/RDPTools/Xander_assembler/gene_resource/nirK/for_enone.hmm share/RDPTools/Xander_assembler/gene_resource/nirK/rev_enone.hmm share/RDPTools/Xander_assembler/gene_resource/nirK/ref_aligned.faa

実際には付属のXanderの中のmegagta/share/RDPTools/Xander_assembler/gene_resource/に、該当遺伝子を、たくさんの種のマルチプルアライメントから計算しHMMファイルが入っている(nifやrpoBなど)。ターゲットがこれらなら、そのまま指定すればよい。

 

該当遺伝子のHMMファイルが見つかれなければ作成する必要がある。手順はXanderのHマニュアルに書かれている。

https://github.com/rdpstaff/Xander_assembler/blob/master/README.md#per-gene-preparation-requires-biological-insight

 

ターゲットアセンブリを実行する。ターゲットのconfigファイルを"-g gene_list.txt"で指定している。

python megagta.py -r input.fq -g gene_list.txt -o output_dir

 出力

ls -lth output_dir/

$ ls -lth output_dir/

total 52K

-rw-r--r-- 1 uesaka user  20K Jun 10 18:06 log

drwxr-xr-x 4 uesaka user 4.0K Jun 10 18:06 contigs

drwxr-xr-x 7 uesaka user 4.0K Jun 10 18:06 .

drwxr-xr-x 2 uesaka user 4.0K Jun 10 18:06 k44

drwxr-xr-x 2 uesaka user 4.0K Jun 10 18:06 k35

drwxr-xr-x 2 uesaka user 4.0K Jun 10 18:05 k29

drwxr-xr-x 2 uesaka user 4.0K Jun 10 18:05 tmp

-rw-r--r-- 1 uesaka user   72 Jun 10 18:05 opts.txt

drwxr-xr-x 4 uesaka user 4.0K Jun 10 18:05 ..

contigs/の中にターゲット遺伝子のfastaが出力される。

ls -lth output_dir/contigs/nirK/

$ ls -lth output_dir/contigs/nirK/

total 0

-rw-r--r-- 1 uesaka user 0 Jun 10 18:06 prot_merged.fasta

-rw-r--r-- 1 uesaka user 0 Jun 10 18:06 nucl_merged.fasta

 

--クラスタリング--

1塩基でも異なると別出力されるため、かなり多くの配列ができる。付属スクリプトpost_proc.shを使い、例えば塩基配列同一性が99%以上の配列をまとめる。

post_proc.shラン前に、 post_proc.shの先頭部分のUCHIMEとHMMALIGNのパスを変更する。

$ head -n 11 post_proc.sh 

#!/bin/bash

 

## This script runs contigs search and post-assembly processing

## The search step requires large memory for large dataset, see Readme for instructions

## This step assumes a gene directory already exists with gene_start.txt in the directory for each gene in the genes list

## This will overwrites the previous search and post-assembly results

SCRIPT=$(readlink -f $0)

SCRIPTPATH=`dirname $SCRIPT`

JAR_DIR=${SCRIPTPATH}/../share/RDPTools/

UCHIME=/usr/local/bin/uchime

HMMALIGN=~/.linuxbrew/bin/hmmalign

準備できたらクラスタリングを実行する。

./post_proc.sh -g megagta/share/RDPTools/Xander_assembler/gene_resource/ -d output_dir/contigs -m 16G -c 0.01

 

 

引用

MegaGTA: a sensitive and accurate metagenomic gene-targeted assembler using iterative de Bruijn graphs.

Li D, Huang Y, Leung CM, Luo R, Ting HF, Lam TW.

BMC Bioinformatics. 2017 Oct 16;18(Suppl 12):408.

 

Xander: employing a novel method for efficient gene-targeted metagenomic assembly

Qiong Wang, Jordan A. Fish, Mariah Gilman, Yanni Sun, C. Titus Brown, James M. Tiedje and James R. Cole

Microbiome. 2015 Aug 5;3:32.