macでインフォマティクス

macでインフォマティクス

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

アセンブルされた微生物ゲノムのクオリティ評価を行う miComplete

2022/03/18 インストール手順追記  

 

 ハイスループットシーケンスの開発により、大規模なシーケンスプロジェクトが手頃な価格になり、可用性がますます向上している。膨大な量のメタゲノムデータが生成され、未培養微生物から数千のmetagenome-assembled genomes(MAG)が公開されている(例:Parks et al、2017)。未培養微生物のゲノムは、フローサイトメーターで細胞を選別し、DNAを増幅およびシーケンシングすることでも取得できる。結果として得られるsingle-cell amplified genomes(SAG)も、微生物の多様性に関する知識を大幅に高めることに貢献している(例えば、Rinke et al、2013)。
 MAGとSAGの品質を評価することが最も重要になった。メタゲノムデータを正しくアセンブリしてMAGにまとめることは難しく、ほとんどのMAGはコンティグが欠落している、および/または外部コンティグで汚染されている。 SAGについては、その完全性は大きく異なり、100%完全なSAGは事実上ありえない。ほとんどの方法は、単一コピーの保存されたマーカー遺伝子を識別することにより、SAGおよびMAGの完全性および汚染(または冗長性)のレベルを推定する。同定されたマーカーの割合はゲノムの完全性に対応し、追加のコピーは汚染または冗長性を表す(Rinke et al、2013)。このアプローチはCheckM(Parks et al、2015)およびBUSCO(Simãoet al、2015)に実装されている。
 これまでのところ、すべてのマーカーは、完全性または冗長性に等しく貢献していると見なされている。しかし、マーカーは原核生物のchromosomeに均一に分布しておらず、一定の量の連鎖が長い進化距離にわたって保存されている(例:Rogozin et al、2002 link)。一般的に使用されるマーカーセットには、他の非クラスター化遺伝子よりも保存されたオペロンで構成されたリボソームタンパク質遺伝子が含まれることが多いため、これは特に重要である。

 miCompleteでは、遺伝的にリンクされたマーカーの有無によって生じる潜在的なバイアスを減らす方法を実装する。 選択されたマーカーのセットは、代表的な完全な閉じたchromosomeセットで最初に識別される。 各ゲノムの各マーカーについて、上流および下流の次のマーカーまでの距離の半分が記録され、2つの合計がゲノムサイズで正規化される。 次に、クエリゲノムの完全性または冗長性を計算するときに、すべてのゲノムの中央値が重みとして使用される。 105の高度に保存された細菌マーカーのセットでは、すべての細菌(以下を参照)を代表する1175のゲノムのセットから推測される重みは、1.02e-5 (Ribosomal protein L22) から 3.17e-2 (RecR) の3桁以上の範囲だった 。 この特定のセットの重みの分布を論文補足図1に示す。重みの標準偏差に対する合計の平方根は、特定の重みのセットに付随する不確実性の尺度として使用される。

(以下略)

 

インストール

ubuntu18のpython3.8環境でテストした。

依存

External software

  • HMMER3 ( Tested with v. 3.1b2)
  • prodigal (Tested with v. 2.6.3)

Python libraries

  • Biopython (>= 1.70) ($ pip install biopython)
    Numpy (>= 1.13.1) ($ pip install numpy)
    Matplotlib (>= 2.0.2) ($ pip install matplotlib)

本体 Bitbucket

#pip
pip install micomplete==1.1.1

#virtual environment
python3 -m venv micomplete
source micomplete/bin/activate

#condaを使う(自分はこちらの手順で導入した)
mamba create -n micomplete python=3.8 -y
conda activate micomplete
#依存
mamba install -c bioconda -y prodigal
mamba install -c bioconda -y hmmer
#本体
pip install micomplete

miComplete -h

$micomplete -h

usage: miComplete [-h] [--lenient] [--format {fna,faa,gbk}] [--hlist [HLIST]]

                  [--hmms HMMS] [--weights WEIGHTS] [--linkage]

                  [--linkage-cutoff LINKAGE_CUTOFF] [--evalue EVALUE]

                  [--bias BIAS] [--domain-cutoff DOMAIN_CUTOFF]

                  [--cutoff CUTOFF] [--threads THREADS] [--log LOG] [-v]

                  [--debug] [--version] [-o OUTFILE]

                  sequence_tab

 

Quality control of metagenome assembled genomes. Able to gather relevant

statistics of completeness and redundancy of genomes given sets of marker

genes, as well as weighted versions of these statstics (including defining new

weights for any given set).

 

positional arguments:

  sequence_tab          Sequence(s) along with type (fna, faa, gbk) provided

                        in a tabular format

 

optional arguments:

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

  --lenient             By default miComplete drops hits with too high bias or

                        too low best domain score. This argument disables that

                        behavior, permitting any hit which meet the evalue

                        requirements.

  --format {fna,faa,gbk}

                        This argument should be used when a single sequence

                        file is given in place of tabulated file of sequences.

                        The argument should be followed by the format of the

                        sequence.

  --hlist [HLIST]       Write list of Present, Absent and Duplicated markers

                        for each organism to file

  --hmms HMMS           Specifies a set of HMMs to be used for completeness

                        check or linkage analysis. The default sets, "Bact105"

                        and "Arch131", can be called via their respective

                        names.

  --weights WEIGHTS     Specify a set of weights for the HMMs specified. The

                        default sets, "Bact105" and "Arch131", can be called

                        via their respective names.

  --linkage             Specifies that the provided sequences should be used

                        to calculate the weights of the provided HMMs

  --linkage-cutoff LINKAGE_CUTOFF

                        Cutoff fraction of the entire fasta which needs to be

                        contained in a single contig in order to be included

                        in linkage calculations. Disabling this is likely to

                        result in some erroneous calculations.

  --evalue EVALUE       Specify e-value cutoff to be used for completeness

                        check. Default = 4e-10

  --bias BIAS           Specify bias cutoff as fraction of score as defined by

                        hmmer. Default = 0.3

  --domain-cutoff DOMAIN_CUTOFF

                        Specify largest allowed difference between full

                        sequence- and domain evalues. Default = 1e-5

  --cutoff CUTOFF       Specify cutoff percentage of markers required to be

                        present in genome for it be included in linkage

                        calculation. Default = 0.9

  --threads THREADS     Specify number of threads to be used in parallel.

                        Default = 1

  --log LOG             Log name. Default=miComplete.log

  -v, --verbose         Enable verbose logging

  --debug

  --version             Returns miComplete version and exits

  -o OUTFILE, --outfile OUTFILE

                        Outfile can be specified. None or "-" will result in

                        printing to stdout

 

Report issues and bugs to the issue tracker at

https://bitbucket.org/evolegiolab/micomplete or directly to eric@hugoson.org

 

 

実行方法

miCompleteは複数のゲノムファイルの評価に適している。まず最初にリストファイルを作る。

find $(realpath .) -maxdepth 1 -type f -name "*.fna" | miCompletelist.sh > list.tab

出力例

f:id:kazumaxneo:20211014221802p:plain

miCompleteはこのリストファイルを使ってゲノムの品質や基本統計を計算する。

 

ゲノムの基本的な統計を出力する。

miComplete list.tab --threads 8 > output

ゲノムのあるパスに遺伝子予測して得たタンパク質配列の.faaファイルが保存されるので注意。

出力例

f:id:kazumaxneo:20211014222052p:plain

(小さなセットだと計算途中では書き込まれないので注意)

 

Completeness

基本統計に加えてゲノムアセンブリの完全性の評価も行う。

miComplete list.tab --hmms Bact105 --threads 12 > output
  • --hmms   Specifies a set of HMMs to be used for completeness  check or linkage analysis. The default sets, "Bact105" and "Arch131", can be called via their respective names.

出力例

f:id:kazumaxneo:20211014222339p:plain

 

Weighted completeness

miComplete list.tab --hmms Bact105 --threads 12 --weights Bact105 > output
  • --weights   Specify a set of weights for the HMMs specified. The default sets, "Bact105" and "Arch131", can be called via their respective names.

出力例

f:id:kazumaxneo:20211014222515p:plain

 

Creating weights

まだ重み付けがされていないマーカー遺伝子のセットを使って、独自のウエイトを作成する。

miComplete list.tab --hmms Bact105 --linkage --threads 12 > Bact105.weights
  • --linkage   Specifies that the provided sequences should be used to calculate the weights of the provided HMMs 

一般的には、できるだけ多くのよく分散した(あるいは、少なくとも重みを使おうとするデータと同じくらい広く分散した)ゲノムを用いて重みを作成するべきである(マニュアルより)。

 

大規模なメタゲノムbin配列群や、NCBIなどからダウンロードした品質のばらつきが大きいドラフトゲノムアセンブリなどをまとめて評価する際にとても有用な方法だと思います。手持ちの2900ゲノムの評価にかかった時間は20分程度でした(8スレッド使用)。

引用

miComplete: weighted quality evaluation of assembled microbial genomes

Eric Hugoson, Wai Tin Lam, Lionel Guy
Bioinformatics, Published: 22 August 2019

 

関連