NBT群落物种级丰度鉴定新方法sylph

简介

Sylph是一个程序（https://github.com/bluenote-1577/sylph），用于对元基因组随机抽样样本执行超快速（1）平均核苷酸一致性（ANI）查询或（2）元基因组分析。

ANI查询：Sylph可以搜索您的样本中的基因组，例如大肠杆菌。如果Sylph输出97%的ANI估计值，则表示您的样本中包含一个与查询基因组具有97%ANI的大肠杆菌。
元基因组分析：Sylph可以确定您样本中的物种/分类群及其丰度，类似于Kraken或MetaPhlAn。能够在几秒钟内对1 Gbp的小鼠肠道读取进行针对85,205个基因组的分析。

为什么选择Sylph？

• 精确的物种级分析：我们的测试表明，Sylph的假阳性率低于Kraken，并且与标记基因方法（如MetaPhlAn, mOTUs）一样精确和敏感。
• 超快速、多线程、多样本处理：Sylph在处理多个样本时比其他方法快50倍以上。在对整个GTDB-R220数据库（110k个基因组）进行分析时，Sylph仅需约15GB的RAM。
• 准确的（包含）ANI信息：Sylph能够在覆盖度低至0.1x的情况下，提供参考基因组和您的元基因组样本之间准确的ANI估计值。
• 可定制的数据库和预建数据库：我们提供原核生物、病毒、真核生物的预建数据库。构建自定义数据库（例如，使用您自己的MAGs）也很简单。可以将分类学信息纳入下游的传统分析报告。
• 支持短读长和长读长：虽然Sylph主要针对短读长进行了基准测试，但在Oxford Nanopore的独立基准测试中，Sylph也是最准确的方法。

Sylph的工作原理

sylph 的新颖之处在于使用基于kmer统计技术对低覆盖率基因组的 ANI 进行校正，ANI大于95%的则归为同种，从而为低丰度基因组提供准确的结果。 有关 sylph 能做什么和不能做什么的更多信息，请参见此处。

• 图 1，使用 sylph sketch 选项将读数和参考基因组分解成 k-mers。k-mers 按 c 的分数（默认 = 200）进行下采样。
• 图 1，使用 sylph query 或 sylph profile，将每个参考基因组中的 k-mers 与读数中的 k-mers 进行核对。
• 图2，Sylph 使用统计数据估算每个参考基因组与元基因组之间的包含 ANI。
• 图 3，sylph 查询：报告上一步中所有 ANI 较高（默认值大于 90%）的基因组。 没有丰度。
• 图 3，sylph profile：如果 ANI > 95%，则使用 k-mer 重映射算法计算丰度并报告物种级的现有基因组。

几秒钟内对 85,205 个基因组分析 1 Gbp 的小鼠肠道读数

Install

(base) [yut@io02 GGG4622_add_jmk423_dRep0.95]$ mamba --version
mamba 1.5.8
conda 24.3.0[yut@n60 ~]$ mamba create -n sylph  -c bioconda sylph       

• database

# download GTDB-R220 pre-built database (~13 GB)
wget http://faust.compbio.cs.cmu.edu/sylph-stuff/gtdb-r220-c200-dbv1.syldb[yut@workstation Database]$ ll gtdb-r220-c200-dbv1.syldb
-rw-rw-r-- 1 yut yut 13G Oct 10 19:30 gtdb-r220-c200-dbv1.syldb# /datanode02/yut/Database/SylphDB//gtdb-r220-c200-dbv1.syldb # 兰州服务器# 安装解析成gtdb脚本
(sylph) [yut@io02 results]$ conda install -c bioconda sylph-tax
# 下载50M数据库
(sylph) [yut@io02 Database]$ sylph-tax download --download-to /datanode02/yut/Database/SylphDB/

使用

• 几个样本

[yut@n101 Batch1_192Samples_20240801]$ time sylph profile /data/users/yut/Database/gtdb-r220-c200-dbv1.syldb -1 Trimgalore/ZXGL2_1.fq.gz Trimgalore/2208QL1_1.fq.gz -2 Trimgalore/ZXGL2_2.fq.gz  Trimgalore/2208QL1_2.fq.gz -t 10 -o profiling1.tsv 2>test_sylph.log&[yut@n101 Batch1_192Samples_20240801]$ jobs
[1]+  Running                 time sylph profile /data/users/yut/Database/gtdb-r220-c200-dbv1.syldb -1 Trimgalore/ZXGL2_1.fq.gz -2 Trimgalore/ZXGL2_2.fq.gz -t 10 > profiling.tsv &
2024-10-10T11:45:24.582Z INFO  [sylph::contain] Trimgalore/ZXGL2_1.fq.gz taxonomic profiling; reassigning k-mers for 287 genomes...
2024-10-10T11:45:28.775Z INFO  [sylph::contain] Trimgalore/ZXGL2_1.fq.gz has 194 genomes passing profiling threshold.
2024-10-10T11:45:28.783Z INFO  [sylph::contain] Finished paired sample Trimgalore/ZXGL2_1.fq.gz.
2024-10-10T11:45:28.783Z INFO  [sylph::contain] sylph finished.
real    7m38.125s
[1]+  Done                    time sylph profile /data/users/yut/Database/gtdb-r220-c200-dbv1.syldb -1 Trimgalore/ZXGL2_1.fq.gz -2 Trimgalore/ZXGL2_2.fq.gz -t 10 > profiling.tsv

• 多个，通配符

time sylph profile /data/users/yut/Database/gtdb-r220-c200-dbv1.syldb -1 Trimgalore/*_1.fq.gz  -2 Trimgalore/*_2.fq.gz  -t 40 -o 192Samples_sylph_abundance/192Samples_sylph_abundance.profiling 2>192Samples_sylph.log

解析成gtdb格式

• 1.解析成gtdb格式的一个样本一个表

(sylph) [yut@node02 Others]$ mkdir  sylph_tables; time sylph-tax taxprof 66Samples_sylph_abundance.profiling  -o sylph_tables/ -t  /datanode02/yut/Database/SylphDB/gtdb_r220_metadata.tsv.gz

• 2.合并表格

(sylph) [yut@node02 sylph_tables]$ cd sylph_tables
(sylph) [yut@node02 sylph_tables]$ sylph-tax merge *sylphmpa --column relative_abundance -o 66Samples_sylph_output_table.tsv```# 输出
Sylph 输出一个 TSV（以制表符分隔的值）文件。 每一行是元基因组样本中检测到的一个基因组。```bash
Sample_file   Genome_file   Taxonomic_abundance   Sequence_abundance   Adjusted_ANI   Eff_cov   ANI_5-95_percentile   Eff_lambda   Lambda_5-95_percentile   Median_cov   Mean_cov_geq1   Containment_ind   Naive_ANI   Contig_name
reads.fq   genome.fa   78.1242   81.8234   97.53   264.000   NA-NA   HIGH   NA-NA   264   264.143   10281/22299   97.53   NC_016901.1 Shewanella baltica OS678, complete genome

• Sample_file：读数/样本的文件名
• Genome_file：检测到的基因组的文件名
• Taxonomic_abundance：归一化分类丰度百分比。 覆盖率归一化–与 MetaPhlAn 丰度相同 不用于 sylph 查询
• Sequence_abundance：归一化序列丰度百分比。 分配给每个基因组的 “读数百分比”–与 Kraken 丰度相同 不适用于 sylph 查询
• Adjusted_ANI：调整后的覆盖率 ANI 估计值。 如果可以调整覆盖率（cov < 3x cov）：返回调整后的覆盖率 ANI 如果覆盖率过低/过高：返回 Naive_ANI（见下文 Eff_cov/True_cov：有效覆盖率估计值，如果指定 -u，则返回真实覆盖率估计值。 总是小数。
• ANI_5-95_percentile：[5%,95%] 置信区间。 如果可以调整覆盖率：浮动-浮动，例如 98.52-99.55 如果覆盖率过低/过高：
• Eff_lambda：有效覆盖率参数的估计值。 如果可以调整覆盖率：给出 lambda 估计值 如果覆盖率过低/过高，则输出 LOW 或 HIGH
• Lambda_5-95_percentile: lambda 的 [5%, 95%] 置信区间。 与 ANI_5-95_percentile 的格式规则相同。
• Median_cov：多重性 >= 1 的 k-mer多重性中位数。
• Mean_cov_geq1：多重性 >= 1 的 k-mer的平均 k-mer多重性。
• Containment_ind：显示包含指数（样本中发现的 k-mer数除以 k-mer总数）的 int/int，例如 959/1053。例如 959/1053。
• Naive_ANI：不进行覆盖率调整的包含 ANI。 kmers_reassigned：从基因组中重新分配的 k-mers 数量。 sylph 查询中不存
• Contig_name：基因组中第一个contig的名称（或 -i 选项中的等位组名称）。

sylph 能做什么？

• 定量元基因组：sylph 可以利用参考数据库计算样本中基因组的丰度，这与 Kraken 或 MetaPhlAn 的输出类型相同。
• 根据元基因组搜索基因组：sylph 可以检查您的样本中是否含有某个基因组（例如，我的样本中是否含有这个大肠杆菌基因组？
• ANI 查询：sylph 可以估算参考基因组与您样本中基因组的平均核苷酸同一性（ANI）。
• 使用自定义参考数据库： 真核生物、病毒和任何 fasta 文件集都可以使用。
• 可使用长读数：sylph 可以高精度地利用纳米孔或 PacBio 读数。 最近的一项研究发现，在他们的数据中，sylph 是最精确的定量方法
• 计算覆盖率：sylph 可以估算数据库中基因组的覆盖率（而不仅仅是丰度）。
• 计算数据库中检测到的物种级读数百分比：sylph 可以检查数据库 "捕获 "了多少元基因组。

sylph 不能做什么？

• 比对reads: 与 Kraken 不同，sylph 不会对每个读数进行分类。
• 查找超低丰度基因组。 Sylph 要求细菌基因组的覆盖率至少 > 0.01-0.05 倍。 所有细菌基因组都需要至少几百个短读数。
• 能可靠地找到属或属以上级别的基因组（如果没有种级别的基因组）。 如果您的样本没有被数据库很好地描述，sylph 可能会很吃力。 注意：这也适用于大多数定量软件。
• Compare genomes to genomes / metagenomes to metagenomes / contigs to genomes
• 处理 16S / ITS 数据

ANI定义

在 sylph 中，ANI 暗指基因组与元基因组之间的包含 ANI。 包含 ANI 是根据元基因组中包含的参考基因组中的 k-mers 数量计算得出的。 如果 "元基因组 "是单个基因组，则包含 ANI 近似于标准 ANI。 如果 "元基因组 "是基因组集合，则包含 ANI 可解释为 “近邻 ANI”（profiling）： 如果 "元基因组 "是从基因组集合中读取的，则 sylph 会用统计模型估算出调整后（包含）的 ANI，该 ANI 与情况 2 相同。 注意：包含 ANI 会略微高估真实 ANI。 参见我们论文中的补充图。

如何使用 sylph-utils 生成包含分类信息的配置文件

默认情况下，sylph 的 TSV 输出不包含分类信息。 不过，sylph 支持整合几个主要的基因组数据库，以获得分类概况。 在 sylph-utils 资源库中，sylph_to_taxprof.py 脚本的使用方法如下
在 sylph-utils 软件仓库中，sylph_to_taxprof.py 脚本的使用方法如下

git clone https://github.com/bluenote-1577/sylph-utils
sylph profile  gtdb-r220-c200-dbv1.syldb sample.sylsp -o result.tsv
python sylph-utils/sylph_to_taxprof.py -s result.tsv -m sylph-utils/prokaryote/gtdb_r220_metadata.tsv.gz
ls *.sylphmpa# 兰大超算
mkdir results && cd results # 先创建输出目录，切换到对应目录再分析
[yut@n101 Batch1_192Samples_20240801]$ time sylph_to_taxprof.py -s profiling.tsv -m ~/Software/sylph-utils/prokaryote/gtdb_r220_metadata.tsv.gz
Reading metadata files: ['/home/yut/Software/sylph-utils/prokaryote/gtdb_r220_metadata.tsv.gz'] ...
Processing sylph output file: profiling.tsv ...
Writing output to: ZXGL2_1.fq.gz.sylphmpa ...real    0m51.314syut@n101 192Samples_sylph_abundance]$ head ../192Samples_sylph.log
2024-10-10T23:25:21.267Z INFO  [sylph::contain] Obtaining sketches...
2024-10-10T23:25:56.852Z INFO  [sylph::contain] Finished obtaining genome sketches.
2024-10-10T23:31:23.396Z INFO  [sylph::contain] Trimgalore/2208QL3_1.fq.gz taxonomic profiling; reassigning k-mers for 183 genomes...
2024-10-10T23:31:49.080Z INFO  [sylph::contain] Trimgalore/2208QL3_1.fq.gz has 110 genomes passing profiling threshold.
2024-10-10T23:31:49.139Z INFO  [sylph::contain] Finished paired sample Trimgalore/2208QL3_1.fq.gz.
2024-10-10T23:32:16.937Z INFO  [sylph::contain] Trimgalore/2305QL2_1.fq.gz taxonomic profiling; reassigning k-mers for 195 genomes...
2024-10-10T23:32:16.937Z INFO  [sylph::contain] Trimgalore/2304QL1-MIX_1.fq.gz taxonomic profiling; reassigning k-mers for 287 genomes...
2024-10-10T23:32:16.937Z INFO  [sylph::contain] Trimgalore/2305QL2-16m_1.fq.gz taxonomic profiling; reassigning k-mers for 107 genomes...
2024-10-10T23:32:46.438Z INFO  [sylph::contain] Trimgalore/2305QL2-16m_1.fq.gz has 48 genomes passing profiling threshold.
2024-10-10T23:32:46.441Z INFO  [sylph::contain] Finished paired sample Trimgalore/2305QL2-16m_1.fq.gz.
024-10-11T00:27:22.802Z INFO  [sylph::contain] Trimgalore/LGCL3-S_1.fq.gz has 66 genomes passing profiling threshold.
2024-10-11T00:27:22.803Z INFO  [sylph::contain] Finished paired sample Trimgalore/LGCL3-S_1.fq.gz.
2024-10-11T00:27:22.803Z INFO  [sylph::contain] sylph finished.

• 192个宏基因组样本40核心1h跑完
  gtdb_r220_metadata.tsv.gz 文件将 GTDB 分类法与每个基因组联系起来。 它将所有结果按样本名称分组，并输出一个名为 sample.sylphmpa 的文件，格式如下所示。
  MetaPhlAn-like or CAMI-like outputs
  *.sylphmpa files look like this:

#SampleID       biofilm_reads/SRR24442552_1.fastq.gz
clade_name     relative_abundance      sequence_abundance      ANI (if strain-level)
d__Bacteria     99.99999999999997       3.4775000000000014      NA
d__Bacteria|p__Actinomycetota   1.7266  0.0414  NA
d__Bacteria|p__Actinomycetota|c__Acidimicrobiia 1.3727  0.0313  NA
d__Bacteria|p__Actinomycetota|c__Acidimicrobiia|o__Acidimicrobiales     1.3727  0.0313  NA
d__Bacteria|p__Actinomycetota|c__Acidimicrobiia|o__Acidimicrobiales|f__Ilumatobacteraceae       1.3727  0.0313  NA
d__Bacteria|p__Actinomycetota|c__Acidimicrobiia|o__Acidimicrobiales|f__Ilumatobacteraceae|g__Casp-actino5       1.0267  0.0172  NA
d__Bacteria|p__Actinomycetota|c__Acidimicrobiia|o__Acidimicrobiales|f__Ilumatobacteraceae|g__Casp-actino5|s__Casp-actino5       1.0267  0.0172  NA
d__Bacteria|p__Actinomycetota|c__Acidimicrobiia|o__Acidimicrobiales|f__Ilumatobacteraceae|g__Casp-actino5|s__Casp-actino5|t__GCA_017859785.1    1.0267  0.0172  98.54
d__Bacteria|p__Actinomycetota|c__Acidimicrobiia|o__Acidimicrobiales|f__Ilumatobacteraceae|g__Ilumatobacter      0.346   0.0141  NA
....

提示 这是一个有效的 TSV 文件，但前缀为 # 的行是注释。 您可以用 python 读取 .sylphmpa 文件，如 pd.read_csv(‘output.sylphmpa’,sep=‘\t’,comment=‘#’)。 有五个重要列：clade_name：描述支系的字符串，如 d__Bacteria|p__Actinomycetota|c__Acidimicrobiia|o_Acidimicrobiales|f__Ilumatobacteraceae。

• relative_abundance（相对丰度）：该支系在分类学上的相对丰度
• sequence_abundance（序列丰度）：该支系的序列丰度，即分配 ANI 的读数百分比：除菌株水平（t__strain）外，该值为 NA。
• Coverage：（sylph_to_taxprof.py 0.2 版新增）这是 sylph 输出中的 Eff_cov 或 True_cov 列。
  在 sylph_to_taxprof.py v0.2 中，病毒宿主信息可用于 IMG/VR 4.1。 -a选项在.sylphmpa文件中增加了一列将病毒基因组与其宿主相关联的新内容。 例如，一行可以是这样的：r__Duplodnaviria|k__Heunggongvirae|p__Uroviricota|c__Caudoviricetes|||||t__IMGVR_UViG_2503982007_000001 … d__Bacteria;p__Firmicutes;c__Bacilli;o__Staphylococcales;f__Staphylococcaceae;g__Staphylococcus;s__Staphylococcus epidermidis 其中 IMGVR_UVIG_2503982007 的宿主是表皮葡萄球菌。

sylph_to_taxprof.py 文件非常简单。 简而言之：-m 选项指定的 gtdb_r220_metadata.tsv 文件（可以压缩，也可以不压缩）是一个有两列的文件。 第一列是你的 fasta 文件的名称，第二列是一个分类字符串，看起来像 d__Archaea;p__Methanobacteriota_B;c__Thermococcci;o__Thermococcales;f__Thermococcaceae;g__Thermococcus_A;s__Thermococcus_A alcaliphilus。 假设你想在 GTDB 分类中添加 genome1.fa 和 genome2.fa，两个新的 MAG。 您可以通过添加两行新内容，将 gtdb_r220_metadata.tsv 文件改成这样：
…
…
GCA_945889495.1 d__Bacteria;p__Desulfobacterota_B;c__Binatia;o__UBA9968;f__UBA9968;g__DP-20;s__DP-20 sp945889495
GCA_934500585.1 d__Bacteria;p__Bacillota_A;c__Clostridia;o__UBA1381;f__UBA1381;g__RQCD01;s__RQCD01 sp008668455
genome1.fa d__Archaea;(…);s__My new species name genome2.fa d__Bacteria;(...);g__My_genus_name;s__My species name2

并将其用于 python 脚本中的 -m 选项。 不要添加 t__STRAIN 行，脚本会自动添加。 警告 对于 Genbank/RefSeq 汇编，必须小心处理文件名。 如果基因组文件名中含有 _genomic 或 _ASM，例如 GCF_002863645.1_ASM286364v1_genomic.fna.gz，请使用 “GCF_002863645.1”，而不是第一列中的整个文件名。

耗时：66个样本耗时1h

转成easymicroplot可用数据

# 只要样本和相对丰度列
[yut@workstation raw_sylphmpa]$ for i in *.sylphmpa;do awk '{print $1"\t"$2}' $i >../processed_for_easymicroplot/${i};done# 合并表格，并且替换|为;并去掉样本后缀名称
[yut@workstation 192Samples_sylph_abundance]$ Merged_abundance_table_by_tab.sh processed_for_easymicroplot/ |sed "s/|/;/g"|sed "s/_1.fq.gz.sylphmpa//g" >192Samples_sylph_mpa4.tsv
# 去掉多余行
[yut@workstation 192Samples_sylph_abundance]$ Merged_abundance_table_by_tab.sh processed_for_easymicroplot/ |sed "s/|/;/g"|sed "s/_1.fq.gz.sylphmpa//g"|awk 'NR!=2 && NR!=3'|sed "s/gene/SampleID/"  >192Samples_sylph_mpa4.tsv# 样本名称统一添加字母s_
[yut@workstation 192Samples_sylph_abundance]$ awk 'NR==1 {OFS="\t";for(i=2; i<=NF; i++) $i="s_" $i; print} NR>1 {print}' 192Samples_sylph_mpa4.tsv >192Samples_sylph_mpa4_for_easymicroplot.tsv