メタゲノムをデータベースと照合してプロファイリングすることで、アセンブルが不可能な低存在量でも微生物の検出と定量が可能になる。本著者らは、ゼロインフレートポアソンk-mer統計量を用いてゲノム間平均ヌクレオチド同一性(ANI)を推定し、ANIに基づく分類群の検出を可能にする、種レベルのメタゲノムプロファイラであるsylphを紹介する。Critical Assessment of Metagenome Interpretation II (CAMI2) Marineデータセットでは、sylphはテストした7つの手法の中で最も正確なプロファイリング手法であった。マルチサンプルのプロファイリングにおいて、sylphはKraken2と比較して中央演算処理装置の処理時間を10倍以上短縮し、メモリ使用量を30倍削減した。SylphのANI推定値はアバンダンスとorthogonalなシグナルを提供し、289,232のゲノムに対するANIベースのメタゲノム全体のパーキンソン病(PD)関連研究を可能にし、同時に既知の酪酸-PD関連を菌株レベルで確認した。Sylphは、85,205の原核生物ゲノムと2,917,516のウイルスゲノムに対してメタゲノムをプロファイリングするのに1分未満と16GBのランダムアクセスメモリを要し、RefSeqと比較してヒト腸内のウイルス配列を30倍多く検出した。Sylphは、低カバレッジゲノムでも正確な封じ込めANI推定を行い、正確で効率的なプロファイリングを提供する。
特徴(レポジトリより)
-
sylphは、メタゲノムショットガンサンプルに対して、超高速で(1)ANIクエリサーチ、または(2)メタゲノムプロファイリングを行う。
1)ANIクエリ:sylphはサンプルに対して大腸菌などのゲノムを検索することができる。
2)メタゲノミックプロファイリング: KrakenやMetaPhlAnと同様に、サンプル中の種/菌種とその存在量を決定できる。 -
正確な種レベルのプロファイリング: sylphはKrakenよりも偽陽性が少なく、マーカー遺伝子メソッド(MetaPhlAn、mOTUs)と同程度の精度と感度を持つ。
-
超高速、マルチスレッド、マルチサンプル: sylphは他の方法より50倍以上速い。GTDB-R220データベース全体(110kゲノム)に対するプロファイリングに必要なRAM容量はわずか15GB。
-
正確なANI情報: sylphはリファレンスゲノムとメタゲノムサンプルとの間のANIを0.1xカバレッジまで正確に推定可能
-
カスタマイズ可能なデータベースと構築済みデータベース: 原核生物、ウイルス、真核生物の構築済みデータベースを提供している。カスタムデータベース(独自のMAGを使用するなど)も簡単に構築可能
-
ショートリードでもロングリードでも、Sylphは、Oxford Nanoporeが独自に実施したベンチマークでも最も精度の高い手法だった(link)
Introduction: what is sylph and how does it work?
https://github.com/bluenote-1577/sylph/wiki/Introduction:-what-is-sylph-and-how-does-it-work%3F
Pre-Build databases
https://github.com/bluenote-1577/sylph/wiki/Pre%E2%80%90built-databases
cookbook
https://github.com/bluenote-1577/sylph/wiki/sylph-cookbook
Tutorial
https://github.com/bluenote-1577/sylph/wiki/5%E2%80%90minute-sylph-tutorial
インストール
#本体 (conda)
mamba create -n sylph -y
conda activate sylph
mamba install -c bioconda sylph -y
#sylph-tax (tax assignに使うサポートコマンド)
mamba install -c bioconda sylph-tax -y
#本体 x86-64 binary
wget https://github.com/bluenote-1577/sylph/releases/download/latest/sylph
chmod +x sylph
./sylph -h
> sylph -h
sylph 0.8.1
Ultrafast genome ANI queries and taxonomic profiling for metagenomic shotgun samples.
--- Preparing inputs by sketching (indexing)
## fastq (reads) and fasta (genomes all at once)
## *.sylsp found in -d; *.syldb given by -o
sylph sketch -t 5 sample1.fq sample2.fq genome1.fa genome2.fa -o genome1+genome2 -d sample_dir
## paired-end reads
sylph sketch -1 a_1.fq b_1.fq -2 b_2.fq b_2.fq -d paired_sketches
--- Nearest neighbour containment ANI
sylph query *.syldb *.sylsp > all-to-all-query.tsv
--- Taxonomic profiling with relative abundances and ANI
sylph profile *.syldb *.sylsp > all-to-all-profile.tsv
USAGE:
sylph <SUBCOMMAND>
OPTIONS:
-h, --help Print help information
-V, --version Print version information
SUBCOMMANDS:
sketch Sketch sequences into samples (reads) and databases (genomes). Each sample.fq ->
sample.sylsp. All *.fa -> *.syldb
profile Species-level taxonomic profiling with abundances and ANIs
query Coverage-adjusted ANI querying between databases and samples
inspect Inspect sketched .syldb and .sylsp files
> sylph sketch -h
sylph-sketch
Sketch sequences into samples (reads) and databases (genomes). Each sample.fq -> sample.sylsp. All
*.fa -> *.syldb
USAGE:
sylph sketch [OPTIONS] [--] [FILES]...
OPTIONS:
--debug Debug output
-h, --help Print help information
-t <THREADS> Number of threads [default: 3]
--trace Trace output (caution: very verbose)
INPUT:
-l, --list <LIST_SEQUENCE>
Newline delimited file with inputs; fastas -> database, fastq -> sample
--lS <LIST_SAMPLE_NAMES>
Newline delimited file; read sketches are renamed to given sample names
-S, --sample-names <SAMPLE_NAMES>...
Read sketches are renamed to given sample names
<FILES>...
fasta/fastq files; gzip optional. Default: fastq file produces a sample sketch (*.sylsp)
while fasta files are combined into a database (*.syldb).
OUTPUT:
-d, --sample-output-directory <SAMPLE_OUTPUT_DIR>
Output directory for sample sketches [default: ./]
-o, --out-name-db <DB_OUT_NAME>
Output name for database sketch (with .syldb appended) [default: database]
GENOME INPUT:
-g, --genomes <GENOMES>... Genomes in fasta format
--gl <LIST_GENOMES> Newline delimited file; inputs assumed genomes
-i, --individual-records Use individual records (contigs) for database construction
SINGLE-END INPUT:
-r, --reads <READS>... Single-end fasta/fastq reads
PAIRED-END INPUT:
-1, --first-pairs <FIRST_PAIR>...
First pairs for paired end reads
-2, --second-pairs <SECOND_PAIR>...
Second pairs for paired end reads
--l1 <LIST_FIRST_PAIR>
Newline delimited file; inputs are first pair of PE reads
--l2 <LIST_SECOND_PAIR>
Newline delimited file; inputs are second pair of PE reads
ALGORITHM:
-c <C>
Subsampling rate [default: 200]
--fpr <FPR>
False positive rate for read deduplicate hashing; valid values in [0,1). [default:
0.0001]
-k <K>
Value of k. Only k = 21, 31 are currently supported [default: 31]
--min-spacing <MIN_SPACING_KMER>
Minimum spacing between selected k-mers on the genomes [default: 30]
--no-dedup
Disable read deduplication procedure. Reduces memory; not recommended for illumina data
> sylph profile -h
sylph-profile
Species-level taxonomic profiling with abundances and ANIs
USAGE:
sylph profile [OPTIONS] [--] [FILES]...
ARGS:
<FILES>... Pre-sketched *.syldb/*.sylsp files. Raw single-end fastq/fasta are allowed and
will be automatically sketched to .sylsp/.syldb
OPTIONS:
--debug
Debug output
-h, --help
Print help information
--log-reassignments
Output information for how k-mers for genomes are reassigned during `profile`. Caution:
can be verbose and slows down computation.
-s, --sample-threads <SAMPLE_THREADS>
Number of samples to be processed concurrently. Default: (# of total threads / 3) + 1
for profile, 1 for query
-t <THREADS>
Number of threads [default: 3]
--trace
Trace output (caution: very verbose)
INPUT/OUTPUT:
-l, --list <FILE_LIST> Newline delimited file of file inputs
-o, --output-file <OUT_FILE_NAME> Output to this file (TSV format). [default: stdout]
ALGORITHM:
-I, --read-seq-id <SEQ_ID>
Sequence identity (%) of reads. Only used in -u option and overrides automatic
detection.
-m, --minimum-ani <MINIMUM_ANI>
Minimum adjusted ANI to consider (0-100). Default is 90 for query and 95 for profile.
Smaller than 95 for profile will give inaccurate results.
-M, --min-number-kmers <MIN_NUMBER_KMERS>
Exclude genomes with less than this number of sampled k-mers [default: 50]
--min-count-correct <MIN_COUNT_CORRECT>
Minimum k-mer multiplicity needed for coverage correction. Higher values gives more
precision but lower sensitivity [default: 3]
-u, --estimate-unknown
Estimate true coverage and scale sequence abundance in `profile` by estimated unknown
sequence percentage
SKETCHING:
-r, --reads <READS>...
Single-end raw reads (fastx/gzip)
-1, --first-pairs <FIRST_PAIR>...
First pairs for raw paired-end reads (fastx/gzip)
-2, --second-pairs <SECOND_PAIR>...
Second pairs for raw paired-end reads (fastx/gzip)
-c <C>
Subsampling rate. Does nothing for pre-sketched files [default: 200]
-i, --individual-records
Use individual records (e.g. contigs) for database construction instead. Does nothing
for pre-sketched files
-k <K>
Value of k. Only k = 21, 31 are currently supported. Does nothing for pre-sketched files
[default: 31]
--min-spacing <MIN_SPACING_KMER>
Minimum spacing between selected k-mers on the database genomes. Does nothing for
pre-sketched files [default: 30]
> sylph profile -h
sylph-profile
Species-level taxonomic profiling with abundances and ANIs
USAGE:
sylph profile [OPTIONS] [--] [FILES]...
ARGS:
<FILES>... Pre-sketched *.syldb/*.sylsp files. Raw single-end fastq/fasta are allowed and
will be automatically sketched to .sylsp/.syldb
OPTIONS:
--debug
Debug output
-h, --help
Print help information
--log-reassignments
Output information for how k-mers for genomes are reassigned during `profile`. Caution:
can be verbose and slows down computation.
-s, --sample-threads <SAMPLE_THREADS>
Number of samples to be processed concurrently. Default: (# of total threads / 3) + 1
for profile, 1 for query
-t <THREADS>
Number of threads [default: 3]
--trace
Trace output (caution: very verbose)
INPUT/OUTPUT:
-l, --list <FILE_LIST> Newline delimited file of file inputs
-o, --output-file <OUT_FILE_NAME> Output to this file (TSV format). [default: stdout]
ALGORITHM:
-I, --read-seq-id <SEQ_ID>
Sequence identity (%) of reads. Only used in -u option and overrides automatic
detection.
-m, --minimum-ani <MINIMUM_ANI>
Minimum adjusted ANI to consider (0-100). Default is 90 for query and 95 for profile.
Smaller than 95 for profile will give inaccurate results.
-M, --min-number-kmers <MIN_NUMBER_KMERS>
Exclude genomes with less than this number of sampled k-mers [default: 50]
--min-count-correct <MIN_COUNT_CORRECT>
Minimum k-mer multiplicity needed for coverage correction. Higher values gives more
precision but lower sensitivity [default: 3]
-u, --estimate-unknown
Estimate true coverage and scale sequence abundance in `profile` by estimated unknown
sequence percentage
SKETCHING:
-r, --reads <READS>...
Single-end raw reads (fastx/gzip)
-1, --first-pairs <FIRST_PAIR>...
First pairs for raw paired-end reads (fastx/gzip)
-2, --second-pairs <SECOND_PAIR>...
Second pairs for raw paired-end reads (fastx/gzip)
-c <C>
Subsampling rate. Does nothing for pre-sketched files [default: 200]
-i, --individual-records
Use individual records (e.g. contigs) for database construction instead. Does nothing
for pre-sketched files
-k <K>
Value of k. Only k = 21, 31 are currently supported. Does nothing for pre-sketched files
[default: 31]
--min-spacing <MIN_SPACING_KMER>
Minimum spacing between selected k-mers on the database genomes. Does nothing for
pre-sketched files [default: 30]
> sylph inspect -h
sylph-inspect
Inspect sketched .syldb and .sylsp files
USAGE:
sylph inspect [OPTIONS] [FILES]...
ARGS:
<FILES>... Pre-sketched *.syldb/*.sylsp files.
OPTIONS:
-h, --help Print help information
-o, --output-file <OUT_FILE_NAME> Output to this file (YAML format). [default: stdout]
$ sylph-tax -h
usage: sylph-tax [-h] [--version] {taxprof,merge,download} ...
sylph-tax - A taxonomy management tool for the sylph metagenome profiler
positional arguments:
{taxprof,merge,download}
taxprof Generate a taxonomic profile by integrating taxonomic information into sylph output files
merge Merge taxonomic profiles from sylph-tax taxprof
download Download sylph-compatible taxonomy data for a collection of genomic databases.
optional arguments:
-h, --help show this help message and exit
--version show program's version number and exit
> sylph-tax taxprof -h
usage: sylph-tax taxprof [-h] [-o STRING] -t FILE [FILE ...] [-a] [-f] SYLPH-FILE [SYLPH-FILE ...]
positional arguments:
SYLPH-FILE sylph result files (TSV)
optional arguments:
-h, --help show this help message and exit
-o STRING, --output-prefix STRING
Append this prefix to the outputs. Output files will be 'prefix + Sample_file_column + .sylphmpa'
-t FILE [FILE ...], --taxonomy-metadata FILE [FILE ...]
Taxonomy metadata inputs. If multiple are provided, they will be merged. Provided taxonomies: [FungiRefSeq-2024-07-25, GTDB_r214, GTDB_r220, IMGVR_4.1, OceanDNA, SoilSMAG, TaraEukaryoticSMAG]. Custom metadata files can be used as well; see online manual.
-a, --annotate-virus-hosts
Add additional column(s) by integrating viral-host information available (currently available for IMGVR4.1)
-f, --add-folder-information
Include directory/folder path information in the output files. This is needed if your samples have the same read name but different directory structures.
データベース
https://github.com/bluenote-1577/sylph/wiki/Pre%E2%80%90built-databases
ここではGTDB r220 database (113,104 species representative genomes)を使用する。
#GTDB-R220 pre-built database
wget http://faust.compbio.cs.cmu.edu/sylph-stuff/gtdb-r220-c200-dbv1.syldb
GTDB r220 databaseのファイルサイズは14GB程度。
実行方法
1,複数サンプルのプロファイリング。20スレッド指定。
#paired-end
sylph profile gtdb-r220-c200-dbv1.syldb -1 *_1.fastq.gz -2 *_2.fastq.gz -t 20 > profiling.tsv
#simgle-end
sylph profile gtdb-r220-c200-dbv1.syldb *.fastq -t 20 > profiling.tsv
- -t Number of threads [default: 3]
上記のデータベースをSATA3 HDDに配置し、400MBx2のgzip圧縮ペアエンドfastqをプロファイリングしたところ、22秒で終了した(5995WX、20スレッド指定時)。
出力例
profiling.tsv(右側にさらに列あり)
Taxonomic_abundanceが正規化された分類学的存在量。このサンプルは単離されたバクテリアのシークエンシングリードだが、2行に2つの分類がプリントされている。しかしSequence_abundance:(正規化された配列の存在量。Kraken abundanceと同じ)が1行目で99%以上となっており、1行目の分類が主要であることが分かった。他にもカバレッジ調整後ANI値、Naive_ANI(カバレッジ調整なしのANI)、信頼区間などの統計量がプリントされている。Contig_nameはゲノム中の最初のコンティグ名。Containment_indはcontainment(封じ込め)指数で、サンプルで見つかったk-merの数を全k-merで割ったもの、を示す。メタゲノムが単一ゲノムの場合、封じ込めANIは標準ANIに近似する。メタゲノムがゲノムの集合体である場合、封じ込めANIは「最近傍ANI」として解釈できる(マニュアルより)。
"出力について"
https://github.com/bluenote-1577/sylph/wiki/Output-format
krakenスタイルの分類群情報をアサインするには、sylph profileコマンドの後にsylph-taxコマンドを実行する。
2,初回はDB(metadata.tsv)をダウンロードする。
mkdir sylph-tax_DB
sylph-tax download --download-to ./sylph-tax_DB
sylph-tax_DB
3,sylphのTSVからtaxonomic profilesを発生させる。1で指定したDBを-tで指定する。
sylph-tax taxprof sylph_results/*.tsv -t GTDB_r220
- -o Append this prefix to the outputs. Output files will be 'prefix + Sample_file_column + .sylphmpa'
- -t Taxonomy metadata inputs. If multiple are provided, they will be merged. Provided taxonomies: [FungiRefSeq-2024-07-25, GTDB_r214, GTDB_r220, IMGVR_4.1, OceanDNA, SoilSMAG, TaraEukaryoticSMAG]. Custom metadata files can be used as well; see online manual.
出力例
.sylphmpaファイルがTSVごとに出力される。”-o <prefix>”なしで実行すると、prefixには元のfastqファイル名が使用される。
sample1_R1.fq.gz.sylphmpa
5つの重要なカラムがある:
1,clade_name: クレードを表す。つまり、d__Bacteria|p__Actinomycetota|c__Acidimicrobiia|o__Acidimicrobiales|f__Ilumatobacteraceaeのような文字列
2,relative_abundance: クレードの分類学的相対存在量
3,sequence_abundance:クレードの配列存在度、すなわち割り当てられたリードの割合
4,ANI:株レベル(t__strain)以外はNA。それ以外はsylphのANI推定値
5,カバレッジ: sylphの出力のEff_covまたはTrue_cov列
"出力について"
他のコマンドについてはマニュアルを確認してください。マニュアルやチュートリアルが充実していて、とてもユーザーフレンドリーな作りになっています。
引用
Rapid species-level metagenome profiling and containment estimation with sylph
Jim Shaw, Yun William Yu
Nat Biotechnol. 2024 Oct 8.