macでインフォマティクス

macでインフォマティクス

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

メタゲノムデータから株レベルの多様性を検出する metaSNV

2021 Bioinformatics metagenome SNV genetic distance sub-populations bacteria


 微生物群集の分類学的解析は、種および株レベルで十分に支持されている。しかし、種内には顕著な表現型の多様性が存在し、株は世界的な集団間で広く共有されることは稀である。種と株の間の多様性を層別化することで、「>亜種」と呼ばれる有用な中間層を特定できる。亜種のハイスループットな同定とプロファイリングは、微生物叢分野ではまだサポートされていない。本研究では、種内の一塩基多型（SNV）パターンに基づく操作上の定義を用いて>、メタゲノム中の亜種を同定・プロファイリングし、それらの特徴的なSNVおよび遺伝子を明らかにする。この手法をmetaSNV v2に組み込み、既存のSNV検出ソフトウェアを拡張して、集団遺伝学のためのさらなるSNV解釈を可能に>した。これらの新機能は、SNVプロファイルと宿主の表現型や環境、ニッチ特異性を結びつける微生物叢解析を支援する。海洋および糞便メタゲノムにおける亜種同定を実証し、後者では7524人の成人および乳児の70種を解析し、>人間の腸内微生物叢における一般的な亜種集団構造を支持するとともに、亜種同定のいくつかの限界を示した。ソースコード、ドキュメント、チュートリアル、テストデータは https://github.com/metasnv-tool/metaSNV および https://metasnv.embl.de で入手できる｡

 

インストール

Github

mamba create --name metaSNV -c bioconda -c conda-forge 'metasnv>=2.0.3' -y
conda activate metaSNV
mamba install -c bioconda -c conda-forge 'metasnv>=2.0.3' -y

> metaSNV.py --help

usage: metaSNV.py [-h] [--db_ann DB_ANN_FILE] [--print-commands] [--threads INT] [--n_splits INT] [--use_prev_cov] [--min_pos_cov INT] [--min_pos_snvs INT] DIR FILE REF_DB_FILE

 

Compute SNV profiles

 

positional arguments:

  DIR                   The output directory that metaSNV will create e.g. "outputs". Can be a path.

  FILE                  File with an input list of bam files, one file per line

  REF_DB_FILE           reference multi-sequence FASTA file used for the alignments.

 

options:

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

  --db_ann DB_ANN_FILE  Database gene annotation.

  --print-commands      Instead of executing the commands, simply print them out

  --threads INT         Number of jobs to run simmultaneously. Will create same number of splits, unless n_splits set differently.

  --n_splits INT        Number of bins to split ref into

  --use_prev_cov        Use "cov/" and "outputs.all_cov.tab" and "outputs.all_perc.tab" data produced by previous metaSNV run

  --min_pos_cov INT     minimum coverage (mapped reads) per position for snpCall.

  --min_pos_snvs INT    minimum number of non-reference nucleotides per position for snpCall.

 

> metaSNV_DistDiv.py --help

usage: metaSNV_DistDiv.py [-h] --filt : Filtered frequency files [--dist] [--div] [--divNS] [--matched] [--n_threads : Number of Processes]

 

metaSNV distance and diversity computation

 

options:

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

  --filt : Filtered frequency files

                        Folder containing /*.filtered.freq (default: None)

  --dist                Compute distances (default: False)

  --div                 Compute Diversity and FST (default: False)

  --divNS               Computing piN and piS (default: False)

  --matched             Computing on matched positions only (default: False)

  --n_threads : Number of Processes

                        Number of jobs to run simmultaneously. (default: 1)

 

> metaSNV_Filtering.py --help

usage: metaSNV_filtering.py [-h] [-b FLOAT] [-d FLOAT] [-m INT] [-c FLOAT] [-p FLOAT] [--ind] [--n_threads : Number of Processes] Proj

 

metaSNV filtering step

 

positional arguments:

  Proj                  project name

 

options:

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

  -b FLOAT              Coverage breadth: minimal horizontal genome coverage percentage per sample per species (default: 40.0)

  -d FLOAT              Coverage depth: minimal average vertical genome coverage per sample per species (default: 5.0)

  -m INT                Minimum number of samples per species (default: 2)

  -c FLOAT              FILTERING STEP II:minimum coverage per position per sample per species (default: 5.0)

  -p FLOAT              FILTERING STEP II:required proportion of informative samples (coverage non-zero) per position (default: 0.5)

  --ind                 Compute individual SNVs (default: False)

  --n_threads : Number of Processes

                        Number of jobs to run simultaneously. (default: 1)

 

Note:

 

 

テストラン1

1､テストデータのダウンロードと解凍

bamファイル群､リファレンス､出力ファイルなどが含まれている｡

wget http://swifter.embl.de/~ralves/metaSNV_test_data/testdata.tar.xz
tar xvf testdata.tar.xz && rm -f testdata.tar.xz

解凍したディレクトリ

all_samplesはbamのパスを1行ずつ記載したリストファイル

> head testdata/

testdata/maps/m0.insilicoRefs.unique.sorted.bam

testdata/maps/m1.insilicoRefs.unique.sorted.bam

testdata/maps/m2.insilicoRefs.unique.sorted.bam

testdata/maps/m3.insilicoRefs.unique.sorted.bam

testdata/maps/m4.insilicoRefs.unique.sorted.bam

testdata/maps/m5.insilicoRefs.unique.sorted.bam

testdata/maps/m6.insilicoRefs.unique.sorted.bam

testdata/maps/m7.insilicoRefs.unique.sorted.bam

testdata/maps/m8.insilicoRefs.unique.sorted.bam

testdata/maps/m9.insilicoRefs.unique.sorted.bam

 

以下のようなbamファイルが使われている｡map/に配置されている｡

> ls -l testdata/maps/ |head

total 339360

-rw-r--r-- 1 kazu kazu 3145441  2月 13  2020 m0.insilicoRefs.unique.sorted.bam

-rw-r--r-- 1 kazu kazu 3121649  2月 13  2020 m1.insilicoRefs.unique.sorted.bam

-rw-r--r-- 1 kazu kazu 3143212  2月 13  2020 m2.insilicoRefs.unique.sorted.bam

-rw-r--r-- 1 kazu kazu 3150718  2月 13  2020 m3.insilicoRefs.unique.sorted.bam

-rw-r--r-- 1 kazu kazu 3139221  2月 13  2020 m4.insilicoRefs.unique.sorted.bam

-rw-r--r-- 1 kazu kazu 3109410  2月 13  2020 m5.insilicoRefs.unique.sorted.bam

-rw-r--r-- 1 kazu kazu 3130979  2月 13  2020 m6.insilicoRefs.unique.sorted.bam

-rw-r--r-- 1 kazu kazu 3121314  2月 13  2020 m7.insilicoRefs.unique.sorted.bam

-rw-r--r-- 1 kazu kazu 3141209  2月 13  2020 m8.insilicoRefs.unique.sorted.bam

リファレンスも含まれている｡ref/に配置されている｡

>  ls -l testdata/ref/ |head

total 316

-rw-r--r-- 1 kazu kazu 305049  2月 13  2020 allReferenceGenomes.fasta

-rw-r--r-- 1 kazu kazu  15534  9月 18  2021 allReferenceGenomes_geneAnno.txt

 

2､metaSNVのラン

option､出力ディレクトリ､bamkリスト､リファレンスゲノムの順番で指定する｡計算には1分ほどかかる｡

metaSNV.py --threads 8 output testdata/all_samples testdata/ref/allReferenceGenomes.fasta
  • --threads    Number of jobs to run simmultaneously. Will create same number of splits, unless n_splits set differently.

出力

SNVは現在、output/snpCaller/に書き出される｡SNVがコドンの変化を引き起こすかどうかも検出するには、--db_ann testdata/ref/allReferenceGenomes_geneAnno.txt を追加する(マニュアルより)。all_cov.tabはcontigごとの 平均カバレッジ

 

3､フィルタリング

アリル頻度(変異型の頻度)の計算を行い､低カバレッジ領域のフィルタリングを行う｡

metaSNV_Filtering.py output_test/ | tee filtering.log

フィルタリング条件は標準出力されるのでlogに保存しておく｡

  • breadth ≥ 40%
  • depth ≥ 5x
  • Min. samples per taxid ≥ 2
  • Min. coverage per position ≥ 5x
  • Min. proportion of covered samples ≥ 0.5

(デフォルト設定)

出力

出力ディレクトリにfiltered/が追加される｡

filtered/pop/

 

4､株間距離行列ファイルの作成

このコマンドは、フィルタリングされたSNVのアリル頻度に基づいてサンプル間のペアワイズ非類似度を計算する。

metaSNV_Filtering.py output_test/

サンプル間の遺伝的距離行列や多様性統計量が出力される。

フィルタリングされたSNVアリル頻度はoutput/distancesディレクトリに保存される。

distances/

サンプル間の遺伝的距離行列や多様性統計量が出力されている。リファレンスの種ごとに､160個のbamサンプル分（ヘッダー行を含め161行）の行列ファイルが作成されている。

1つ開いた｡160x160となっている｡

 

5､種内の亜集団クラスターの検出

metaSNV_subpopr.R --procs 3 -i output -g testdata/abunds/geneAbundances.tsv -a testdata/abunds/speciesAbundances.tsv

結果はresults/params.hr10.hs80.ps80.gs80/output/resultsSummary.htmlとresults/params.hr10.hs80.ps80.gs80/output/summary_allResults.csvに記載される。R version 4.4.3で試したときはエラーが発生した｡本来はHTMLレポートが出力され､PCAのグラフなどで亜集団クラスタが可視化される｡

 

テストラン2

レポジトリではProGenomes2のspecies-level representative genomeをリファレンスとした実践的なラン手順も説明されている｡256,346個の配列(コンティグも含む)が含まれている｡

 

1､ゲノムとアノテーションのダウンロード

wget http://swifter.embl.de/~ralves/metaSNV_reference_data/progenomes2_speciesReps_genomes.fna
wget http://swifter.embl.de/~ralves/metaSNV_reference_data/progenomes2_speciesReps_annotations.txt

25GBある｡こちらをリファレンスに使用する｡

 

コメント

metaSNVは 複数サンプル間での株間多様性（サブポピュレーション）を解析する目的で設計されているので､,fastqが1個だけしかないときは得られる情報はかなり限定的になります｡ご注意ください(注; snpCaller/called_SNPsにはSNPs頻度のテキストが配置されているので､これを使ってsampleの株間多様性が高いかどうか評価することはできると思いますが､それならsamtoolsとbcftoolsのほうがシンプルでしょうか)｡

引用

metaSNV v2: detection of SNVs and subspecies in prokaryotic metagenomes Open Access

Thea Van Rossum , Paul I Costea , Lucas Paoli , Renato Alves , Roman Thielemann , Shinichi Sunagawa , Peer Bork

Bioinformatics, Volume 38, Issue 4, February 2022, Pages 1162–1164

 

関連