Skip to content

Instantly share code, notes, and snippets.

@hanzou666
Last active December 18, 2018 06:16
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save hanzou666/2b28f1f107b367d5af1983ce96b7f90c to your computer and use it in GitHub Desktop.
Save hanzou666/2b28f1f107b367d5af1983ce96b7f90c to your computer and use it in GitHub Desktop.
(メタ)ゲノム解析に用いる配列処理ツールの説明

CheckM

MAG(Metagenome-Assembled Genome) や SAG(Single Amplified Genome) のクオリティチェックをバクテリアがもつべき遺伝子が揃っているかどうかでチェックするツール

インストール

3rd party ツール

必要なツールをあらかじめインストールし、パスを通しておく必要がある。

  • HMMER (>=3.1b1)
  • prodigal (2.60 or >=2.6.1)
  • pplacer (>=1.1)

本体

2系のPythonが入っている状態で、

pip install numpy  # 入っていない場合のみ
pip install checkm-genome

リファレンスDBのセットアップ

リファレンスのフォルダをダウンロードしてパスを通しておく必要がある。 ここから最新版のcheckm_data*.tar.gzをダウンロードし、解凍する。そのあと、

checkm data setRoot <解凍したディレクトリのパス>

でCheckMにリファレンスDBの居場所を教えておく。

基本的な使い方

checkm lineage_wf -t 16 --tab_table -f evaluation.tsv -x fa input_file_dir output_dir
  • input_file_dir は、ビニングされたfastaファイルが置いてあるディレクトリを指定する。
    • 拡張子で探すので、 fna 以外の拡張子が使われている場合は、 -x で明示する。
  • evaluation.tsv に各ビンの評価が書かれている。
    • デフォルトでは複数のスペース区切りとなっていて使いにくいので、 --tab_table でタブ区切りにしておく。
    • 1列目: Bin Id
      • fastaファイルの名前から拡張子をとったもの。
    • 2列目: Marker lineage
      • リファレンスツリー上の位置。このクレードでuniversalなマーカー遺伝子セットを用いる。
    • 3列目: # genomes
      • このクレードに位置するリファレンスゲノムの数。
    • 4列目: # markers
      • このクレードのマーカー遺伝子の数。
    • 5列目: # marker sets
      • このクレードのcollacated marker genesの数。 # markers から隣りあう遺伝子同士をクラスタリングしたもの。
    • 6列目: 0
      • クエリゲノム中で観測されなかった遺伝子の数。
    • 7列目: 1
      • クエリゲノム中で1回だけ観測された遺伝子の数。
    • 8列目: 2
    • 9列目: 3
    • 10列目: 4
    • 11列目: 5+
    • 12列目: Completeness
      • 定義を参照
    • 13列目: Contamination
      • 定義を参照
    • 14列目: Strain heterogeneity
      • 定義を参照

定義

Completeness

ある系統内の生物がもつはずのマーカー遺伝子セットのうち、クエリゲノムにはどれだけ存在するか

$$\dfrac{\sum_{s \in M}\dfrac{|s \cap G_{M}|}{|s|}}{|M|}$$

  • $s$ : 近くの位置にあるマーカー遺伝子を1セットとしてまとめたもの。collacated marker genes
  • 「2つの遺伝子が並んでいる」の定義
    • 5kbp以内の距離に隣り合って存在する
  • 近くの位置に例: geneA-geneB, geneB-geneC がそれぞれ並んでいるときは、 $s_{1} = \mathrm{{geneA, geneB, geneC}}$ となる。
  • $M$ : すべてのcollocated marker set
  • $G_{M}$ : クエリのゲノムから予測された全マーカー遺伝子
  • $| \cdot |$ : 集合の要素の個数

範囲は 0 ≤ completeness ≤ 1 で、CheckMでは%表示で報告される

Contamination

ある系統内の生物がもつはずの1コピーしかないマーカー遺伝子が、クエリゲノムには平均してどれだけ余分に存在するか

$$\dfrac{\sum_{s \in M}\dfrac{\sum_{g \in s}C_{g}}{|s|}}{|M|}$$

  • $s$ : 近くの位置にあるマーカー遺伝子を1セットとしてまとめたもの。collacated marker genes
  • $M$ : すべてのcollocated marker set
  • $C_{g}$ : 遺伝子 $g$ が何回現れたか。 $N$ 回カウントされたら $C_{g}=\max(0, N-1)$ である。

範囲は 0 ≤ contamination で、CheckMでは%表示で報告される。上限は1とは限らない。

  • 例:contamination = 600% の場合
    • クエリゲノムの中にコンプリートゲノムが7株入っていたとすると、シングルコピーマーカー遺伝子が平均して7回カウントされる(つまり平均して6個だけ余分に存在する)ので、600% となる。
    • しかし、600%だからといって必ずしも7株分のゲノムが入っているとは限らないので注意。

Strain heterogeneity

同種の複数株や近縁種が1つのfastaファイルに入ってしまっていることによって生じるcontaminationを評価する指標。マルチコピーと判定された遺伝子のAverage Amino acid Identity(AAI) を用いて計算される。

手法の概要

  1. tree : リファレンスツリー上にビンを配置する
    • Prodigalで遺伝子予測
    • HMMERで43個のユニバーサル遺伝子を同定
    • すべての遺伝子をつなげる
    • つなげた遺伝子をpplacerでリファレンスツリー上に配置する
  2. lineage_set : 系統特異的なマーカーセットの同定
  • 1.の結果をもとに、各ビンを評価するのに適した系統マーカー遺伝子セットを同定する
  1. analyze : ビンごとに系統マーカー遺伝子の同定と各指標の計算
  • HMMERでビン中の3.で決めたマーカーセットを同定
  • Completeness, Contamination, Strain heterogeneity を計算する
  1. qa : output
  • いい感じのファイルにまとめる

TODO

  • 系統樹はどうやって書いているのか?
    • 論文の
  • ユニバーサル遺伝子、系統特異的マーカー遺伝子セットはどうやって決めているのか?

GTDB-Tk

手法の概要

  1. identify
    • prodigalで遺伝子予測、HMMERでマーカー遺伝子を同定
  2. align
    • アラインメント済みマーカー遺伝子を連結し、フィルタリングを行う
  3. classify
    • pplacerで連結マーカー遺伝子群をリファレンスツリー上に配置する

実際の系統は、リファレンスツリー上の配置、その相対的な進化的多様性、リファレンスゲノムへのANIを組み合わせてアサインされる。

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment