macでインフォマティクス

macでインフォマティクス

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

メタゲノムのビニングされた真核生物由来コンティグの品質を調べる EukCC

 

 微生物のDNAは日常的に抽出され、配列決定され、ゲノムにアセンブリされている。回収されたゲノムの品質を推定することは、不完全なゲノムや汚染されたゲノムが公表されるのを防ぐために非常に重要である。シングルコピーマーカー遺伝子(SCMG)は、新たにアセンブルされたゲノムの品質を推定するために日常的に使用されている。これらの遺伝子はゲノム内に一度しか存在しないと予想されるため、ドラフトゲノム内で発見されたSCMGの数と予想されるマーカー遺伝子の数を比較することで、完全性の推定が可能となり、マーカー遺伝子の追加コピーの存在は汚染の指標として使用することができる。このアプローチは原核生物でも真核生物でも広く受け入れられている[1,2,3,4]。原核生物のゲノムについては、CheckM [4]が完全性と汚染度を推定するために最も広く使用されているツールだが、他のアプローチも使用されており、BUSCO [2,5]やanvi'o [6]によって原核生物のSCMGのセットも提供されている。CheckMは、ゲノムのクレードを識別するために普遍的なSCMGの初期セットを使用し、その後、品質を推定するためにクレード固有のセットを使用する。

 同様に、SCMGは真核生物の単離ゲノムの評価にも使用されている。CEGMA [1]は、6つのモデル生物から同定された240のユニバーサルシングルコピーマーカー遺伝子を用いてゲノムの完全性を評価している。BUSCOはCEGMAに取って代わり、真核生物の単一普遍マーカー遺伝子セットに加えて、複数の真核生物および原核生物のクレードのマーカー遺伝子セットをキュレーションしたことが大きな進歩であった。しかし、BUSCO(バージョン3.1)は真核生物、原生生物、植物、真菌の完全性を推定するためのセットを提供しているが、ゲノムの品質を評価する際にどのセットが最も適しているかは利用者の選択に委ねられている。これらのより普遍的なアプローチの他にも、特定の真核生物のクレードに焦点を当てたアプローチも存在する。FGMPは、SCMGと真菌ゲノム内の高度に保存された領域の両方を利用して、真菌単離ゲノムの品質を推定している[7]。原生生物については、anvi'oは、真核生物のゲノムの品質を推定するために、BUSCOの「protist」セットからのプロファイルの数を減らしたものを提供している(未発表)。

 過去20年の間に、海や土壌、ヒトの腸のような宿主関連サイトなどの環境(マイクロバイオームとして知られている)に存在する微小な生物の理解が飛躍的に進歩した。このような知見のほとんどは、メタバーコーディング(マーカー遺伝子の増幅)やメタゲノミクス(ショットガンシーケンス)などの手法を用いて、微生物の集合的な遺伝物質に最新のDNAシーケンス技術を適用したことに由来している。このような配列データの解析から、全微生物の最大99%が培養されていないと考えられている[8]。

 これまでのところ、圧倒的な数のメタバーコーディングおよびメタゲノミクス研究は、サンプル内に存在する細菌に焦点を当ててきた。しかし、ウイルスや真核生物もまた、その数と機能の両面で、生態系の重要なメンバーである[9,10,11,12]。実際、単細胞原生生物と真菌は、世界の微生物バイオマスの約17%を占めると推定されている。微生物の真核生物バイオマスの中では、原生生物として知られる遺伝的に多様な単細胞生物が25%もの割合を占めている[13]。

 利用可能な完全なゲノム配列や完全に近いゲノム配列の数が増加しているにもかかわらず、微生物真核生物の解析を含むメタバーコーディングやメタゲノムアプローチは、原生生物の真の多様性が、現在のゲノム参照データベース(RefSeqやENAなど)に反映されているものよりもはるかに大きいことを実証してきた。例えば、メタバーコーディングシークエンシングに基づく最近の推定では、海洋だけでも15万種の真核生物が存在することが示唆されている[14]が、GenBankにエントリーがある代表種はわずか4551種である(15.Nov 2019)。したがって、マイクロバイオームの機能的役割を完全に理解するためには、これらのまだ解明されていない生物が何であるか、そしてそれらが実行している可能性のある機能的役割を知る必要がある。

 現在、マイクロバイオームの機能を理解するための最良のアプローチの一つは、ショットガンリード(通常は200〜500bpの長さ)を集めて、より長いコンティグ(通常は2000〜50万bpの範囲)を得ることである。これらのコンティグは、完全なタンパク質へのアクセスを提供し、そのタンパク質は周囲の遺伝子の文脈の中で解釈される可能性がある。ここ数年、この種の解析を拡張して、メタゲノムアセンブルゲノム(MAG)と呼ばれる仮説ゲノムを回収することが一般的になってきた。MAGは、単一の生物に由来すると考えられるコンティグをグループ化することによって生成されるが、これはビニングと呼ばれるプロセスである[15]。しかし、ビニング後であっても、生物学的な理由(微生物の豊富さなど)、実験的な理由(シーケンシングの深さなど)、技術的な理由(アルゴリズムなど)の組み合わせにより、MAGの完全性にはばらつきがあり、一般的には断片化されている。さらに、コンティグをビニングするために使用される計算手法は、異なる生物に由来するコンティグの区別に失敗することがあり、キメラゲノム(コンタミネーション)につながることがある。これは、MAGを分析する際に考慮しなければならない問題である[16]。そのため、高品質のMAGの同定と選択を可能にするためには、ゲノムの品質評価がプロセスの重要な部分を占めている。

 ここでは、異なる真核生物のクレードにまたがる現在のアプローチの性能を調査し、真核生物のMAGの品質評価にこのツールを適用することを特に視野に入れ、完成度とコンタミネーションの観点から真核生物のゲノム品質を推定するための教師なし手法であるEukCCについて説明する。

 BUSCOでは、特定のクレードのSCMGセットを作成しているため、より正確な品質推定が可能である。この観察を基に、著者らは原生生物や真菌類の広い範囲をカバーする複数のSCMGセットを定義することを目指した。マーカー遺伝子ライブラリの主な用途は、特徴の乏しいゲノムへの応用であり、そのような場合、de novo予測によって遺伝子が同定される可能性が高いと予想している。そこで、メタゲノミクスのユースケースをより正確に表現するために、GeneMark-ESで生成された遺伝子予測を使用して、すべての真核生物のRefSeq種(両生類や維管束植物に属さない)に(再)アノテーションを行うことを選択した。De novoアプローチで予測可能なマーカー遺伝子の選択にも有利である。GeneMark-ESを選択して適用したのは、このツールが大規模な種の範囲でうまく機能し、一般的にRefSeqのアノテーションベンチマークに近い性能を発揮することを以前に実証したからである。さらに、UniProtKBにリファレンスゲノムとして使用されている真核生物の全種を追加した。結果として得られたタンパク質は、hmmer (version 3.2) [22]を用いて、PANTHER 14.1 [21]のファミリーレベルプロファイルHMMを用いてアノテーションを行った。PANTHERを選択したのは、テスト済みのデータベースの中で、分析したタンパク質のカバレッジが最大であることが示されていること [23] と、PANTHERのプロファイルHMMが構成グローバルドメインではなく、完全長のタンパク質ファミリーをモデル化しているからである。

 パラログの分離を高め、共通ドメインによる局所的な一致を最小限に抑えるために、プロファイル固有のビットスコアしきい値を定義することを目的とした。これを達成するために、分類学的にバランスのとれた種のセットに依存し、各プロファイルについて、シングルコピーマッチの数が最も多いビットスコアしきい値を特定した(「方法」のセクションを参照)。

 その後、クレード固有のSCMGを定義するために、まず、55の広く存在するSCMGを用いて、与えられたゲノムの参照ツリーを構築した(ここから「参照セット」と呼ぶ、「方法」の項を参照)。ツリーの各クレードにおいて、率が98%以上のSCMGをチェックした。次に、前述の閾値に一致するクラッド内に20以上のPANTHERファミリーが見つかった場合には、マーカー遺伝子のセットを定義した。この方法を用いて、木全体で477個のSCMGセットを定義することができた。BUSCOやCEGMAとは対照的に、PANTHERモデルに基づいて真核生物キングダム全体に適用可能なSCMGを特定することはできなかったが、代わりに多くのサブクレードに適用可能なセットを発見した。これは特異性を高めるためには望ましいことだが、明らかな欠点は、どのセットを使用するのが最も適切であるかを知ることである。そこで著者らは、最適なSCMGを選択し、それを用いてゲノムの品質を推定するソフトウェアパッケージであるEukCCを開発した。
 分類学的系統が不明な新規ゲノムに対して最も特異的なSCMGのセットを選択するために、EukCCは、55の広く存在するSCMGの参照セットを使用して、de novoで予測されたタンパク質をアノテーションすることにより、初期の分類学的分類を行う。次に、Pplacer [24]を適用して、参照ツリー内の各マッチを系統的に文脈化する。ツリー内の各配置を追跡すると、EukCCは、SCMGセットがデータベースに定義されている最低共通祖先(LCA)ノードを決定する。

 予想されるように、pplacerは多くの場合、すべての配列を参照ツリーの単純で狭い領域に配置するが、時折、整合性のない、遠くに関連するクレード内に配置されることがある。このような場合、SCMGの単一のセットがすべての場所を網羅しているとは限らない。このようなケースを克服するために、配置の最大の割合をカプセル化したSCMGセットが配置される。このプロセスは、不正確な配置または一貫性のない配置のために外れる配置が発生する場合を克服するが、このアプローチは、新規ゲノムからの参照SCMGとの一致がツリー内に確実に配置できない場合、不正確なSCMGセットを選択する可能性がある。これを制御するために、EukCCは常にセット内のプロファイルの数を報告し、配置位置をプロットするオプションを提供している(図5S)。したがって、参照SCMGの一部しか網羅していないセットが選択された状況では、このMAGのより詳細な分析を行うことができるので、これを行うべきである。(以下省略)

 

Documentation

EukCC — EukCC documentation

 

インストール

ubuntu18.04LTSでcondaの仮想環境を作ってテストした。

依存

Github

#1 ここでは3.7の仮想環境を作成
conda create -n eukcc -y python=3.7
conda activate eukcc
conda install -c bioconda -c conda-forge eukcc -y

#GeneMarkはライセンスの関係でHPより直接ダウンロードする
#2 genemakr ESを入れる準備
apt install -y cpanminus make gcc dialog
#3 perl libray
cpanm inc::Module::Install::DSL Hash::Merge MCE::Mutex FindBin Test::Pod Logger::Simple Parallel::ForkManager.pm YAML
#cpanmで入らない場合 => *1

#4 genemakr ESのダウンロード => link

#5 プログラム本体とキーをダウンロードして、
cd GeneMark_dir
#6 Install the key
cp gm_key ~/.gm_key
#7 test run
bash check_install.bash
#8 パスを通す
export PATH=$PWD:$PATH

eukcc -h

$ eukcc -h

usage: eukcc [-h] --db DB [--outdir OUTDIR] [--config CONFIG] [--ncores int]

             [--ncorespplacer int] [--hmm HMM] [--training] [--proteins]

             [--bed file.bed] [--force] [--keeptemp] [--fplace] [--noglob]

             [--quiet] [--debug] [--HPA] [--nPlacements n] [--minGenomes n]

             [--fullineage] [--minPlacementLikelyhood float] [--mindist n]

             [--touch] [--gmes] [--pygmes] [--diamond DIAMOND] [--plot] [-v]

             fasta

 

Evaluate completeness and contamination of a MAG. Args that start with '--'

(eg. --db) can also be set in a config file (specified via --config). Config

file syntax allows: key=value, flag=true, stuff=[a,b,c] (for details, see

syntax at https://goo.gl/R74nmi). If an arg is specified in more than one

place, then commandline values override config file values which override

defaults.

 

positional arguments:

  fasta                 Run script on this bin (fasta file)

 

optional arguments:

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

  --db DB               Path to EukCC DB

  --outdir OUTDIR, -o OUTDIR

                        Location for the output. Names will be prefixed using

                        the bin filenames

  --config CONFIG, -c CONFIG

                        Config file to define parameters, YAML

  --ncores int, -n int  set number of cores for GeneMark-ES, pplacer and Hmmer

  --ncorespplacer int   Pplacer requires a lot of memory. If you want you can

                        set less cores for pplacer, which improves memory

                        consumption significantly

  --hmm HMM             run hmmer on all these HMMs instead

  --training            Run EukCC in training mode (needed to create a new

                        release of the DB)

  --proteins            Input fasta is proteins

  --bed file.bed, -b file.bed

                        You can pass a bedfile of the protein location to omit

                        fragmented proteins being detected twice

  --force, -f           Force rerun of computation even if output is newer

                        than input. Don't resume previous run.

  --keeptemp            Keep all temporary files, by default EukCC will remove

                        some temp files

  --fplace, -p          Force rerun of placement and subsequent steps

  --noglob, -g          Do not expand paths using glob

  --quiet, -q           Silcence most output

  --debug, -d           Debug and thus ignore safety

  --HPA                 Set placement method to HPA

  --nPlacements n       Set number of proteins to support location in tree

                        (default: 2)

  --minGenomes n        Minimal number of genomes to support a set (default:

                        3)

  --fullineage          Output full lineage for MAGs

  --minPlacementLikelyhood float

                        minimal pplacer likelyhood (default: 0.4)

  --mindist n           Distance to collapse hits (default: 2000)

  --touch               Do not run, but touch all output files

  --gmes                only run GeneMark-ES

  --pygmes              Use pygmes, will improve eukccs capability of running

                        on highly fragmented bins but will take longer

  --diamond DIAMOND     required to use pygmes option

  --plot                produce plots

  -v, --version         show program's version number and exit

 

データベースの準備

mkdir Eukcc
cd Eukcc/
wget http://ftp.ebi.ac.uk/pub/databases/metagenomics/eukcc/eukcc_db_v1.1.tar.gz
tar -xzvf eukcc_db_v1.1.tar.gz
mv eukcc_db_20191023_1 eukccdb
#remove tar ball
rm eukcc_db_v1.1.tar.gz

eukcc実行時にはこのeukccdbディレクトリを指定する。 

 

テストラン

22 Contig配列からなる小さなメタゲノムアセンブリを使用する。

wget -O testgenome.fa.gz ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/002/251/995/GCF_002251995.1_ASM225199v2/GCF_002251995.1_ASM225199v2_genomic.fna.gz
gunzip testgenome.fa.gz
eukcc --db eukccdb \
--ncores 4 \
--ncorespplacer 1 \
--outdir eukcc_testgenome \
testgenome.fa
  • --db   Path to EukCC DB
  • --outdir   Location for the output. Names will be prefixed using the bin filenames
  • --ncores    set number of cores for GeneMark-ES, pplacer and Hmmer
  • --ncorespplacer    Pplacer requires a lot of memory. If you want you can set less cores for pplacer, which improves memory consumption significantly

テストデータのランでエラーになる。修正できたら追記します。

 

実際にビニングから始める場合のチュートリアルも用意されています。上のDocumentationリンクを確認して下さい。

引用

Estimating the quality of eukaryotic genomes recovered from metagenomic analysis with EukCC

Paul Saary, Alex L. Mitchell & Robert D. Finn
Genome Biology volume 21, Article number: 244 (2020)


関連


 

*1

cpanmは標準では~/perl5/lib/perl5/に入れるが、@INCに~/perl5/lib/perl5/が含まれていないならライブラリを認識できない。対策は単純に

export PERL5LIB=~/perl5/lib/perl5/

すればよい。