macでインフォマティクス

macでインフォマティクス

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

高速かつ低メモリ使用量でメタゲノムプロファイリングやANIサーチを行う sylph

 

 メタゲノムをデータベースと照合してプロファイリングすることで、アセンブルが不可能な低存在量でも微生物の検出と定量が可能になる。本著者らは、ゼロインフレートポアソン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

 

インストール

Github

#本体 (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

 

"出力について"

https://github.com/bluenote-1577/sylph/wiki/Incorporating-taxonomic-information-into-sylph-with-sylph%E2%80%90tax

 

他のコマンドについてはマニュアルを確認してください。マニュアルやチュートリアルが充実していて、とてもユーザーフレンドリーな作りになっています。

引用

Rapid species-level metagenome profiling and containment estimation with sylph

Jim Shaw, Yun William Yu

Nat Biotechnol. 2024 Oct 8.