Snippy

快速单倍体变异检测与核心基因组比对

作者

Torsten Seemann

概述

Snippy 可在单倍体参考基因组与您的 NGS 测序读段之间寻找 SNPs。它能同时识别替换（snps）和插入/缺失（indels）。在单台计算机上，它可以利用您提供的所有 CPU 核心（已测试至 64 核）。其设计以速度为首要考量，并在单个文件夹中生成一套完整一致的输出文件。此外，它还能对使用相同参考基因组的多组 Snippy 结果进行处理，生成核心 SNP 比对结果（最终可用于构建系统发育树）。

快速开始

% snippy --cpus 16 --outdir mysnps --ref Listeria.gbk --R1 FDA_R1.fastq.gz --R2 FDA_R2.fastq.gz
<cut>
Walltime used: 3 min, 42 sec
Results folder: mysnps
Done.

% ls mysnps
snps.vcf snps.bed snps.gff snps.csv snps.tab snps.html 
snps.bam snps.txt reference/ ...

% head -5 mysnps/snps.tab
CHROM  POS     TYPE    REF   ALT    EVIDENCE        FTYPE STRAND NT_POS AA_POS LOCUS_TAG GENE PRODUCT EFFECT
chr      5958  snp     A     G      G:44 A:0        CDS   +      41/600 13/200 ECO_0001  dnaA replication protein DnaA missense_variant c.548A>C p.Lys183Thr
chr     35524  snp     G     T      T:73 G:1 C:1    tRNA  -   
chr     45722  ins     ATT   ATTT   ATTT:43 ATT:1   CDS   -                    ECO_0045  gyrA DNA gyrase
chr    100541  del     CAAA  CAA    CAA:38 CAAA:1   CDS   +                    ECO_0179      hypothetical protein
plas      619  complex GATC  AATA   GATC:28 AATA:0  
plas     3221  mnp     GA    CT     CT:39 CT:0      CDS   +                    ECO_p012  rep  hypothetical protein

% snippy-core --prefix core mysnps1 mysnps2 mysnps3 mysnps4 
Loaded 4 SNP tables.
Found 2814 core SNPs from 96615 SNPs.

% ls core.*
core.aln core.tab core.tab core.txt core.vcf

安装

Conda

安装 Bioconda，然后：

conda install -c conda-forge -c bioconda -c defaults snippy

Homebrew

安装 Homebrew（适用于 MacOS） 或 LinuxBrew（适用于 Linux），然后：

brew install brewsci/bio/snippy

来源

此命令将直接从 Github 安装最新版本。 您需要将 Snippy 的 bin 目录添加到您的 $PATH 中。

cd $HOME
git clone https://github.com/tseemann/snippy.git
$HOME/snippy/bin/snippy --help

检查安装情况

请确认已安装所需版本：

snippy --version

检查所有依赖项是否已安装并正常运行：

snippy --check

调用 SNPs

输入要求

• FASTA 或 GENBANK 格式的参考基因组（可包含多个 contig）
• FASTQ 或 FASTA 格式的序列读取文件（可为 .gz 压缩格式）
• 用于存放结果的文件夹

输出文件

扩展名	描述
.tab	所有变异的简单制表符分隔摘要
.csv	.tab 文件的逗号分隔版本
.html	.tab 文件的HTML版本
.vcf	VCF格式的最终注释变异
.bed	BED格式的变异
.gff	GFF3格式的变异
.bam	BAM格式的比对结果。包括未比对、多重比对的 reads。不包括重复序列。
.bam.bai	.bam 文件的索引
.log	包含运行的命令及其输出的日志文件
.aligned.fa	参考基因组的一个版本，但在深度=0 的位置用“-”表示，在 0 < 深度 < --mincov 的位置用“N”表示（不包含变异）
.consensus.fa	参考基因组的一个版本，其中所有变异均已实例化
.consensus.subs.fa	参考基因组的一个版本，其中仅替换变异已实例化
.raw.vcf	来自 Freebayes 的未过滤变异调用
.filt.vcf	来自 Freebayes 的已过滤变异调用
.vcf.gz	通过BGZIP压缩的 .vcf 文件
.vcf.gz.csi	通过 bcftools index 生成的 .vcf.gz 索引

⚠️ ❌ Snippy 4.x 不会生成 Snippy 3.x 曾生成的以下文件

扩展名	描述
.vcf.gz.tbi	通过TABIX生成的 .vcf.gz 索引
.depth.gz	.bam 文件的 samtools depth -aa 输出
.depth.gz.tbi	.depth.gz 文件的索引

TAB/CSV/HTML 格式中的列

名称	描述
CHROM	发现变异的序列，例如 FASTA 参考序列中 > 后面的名称
POS	序列中的位置，从 1 开始计数
TYPE	变异类型：snp、mnp、ins、del、complex
REF	参考序列中的核苷酸
ALT	reads 支持的替代核苷酸
EVIDENCE	REF 和 ALT 的频率计数

如果您提供 Genbank 文件作为 --reference（而非 FASTA 文件），Snippy 将通过基因组注释来填充这些额外的列，告知您变异影响了哪些特征：

名称	描述
FTYPE	受影响的特征类别：CDS、tRNA、rRNA……
STRAND	特征所在的链：+、-、.
NT_POS	变异在特征内的核苷酸位置 / 核苷酸长度
AA_POS	氨基酸残基位置 / 氨基酸长度（仅当 FTYPE 为 CDS 时）
LOCUS_TAG	特征的 /locus_tag（如果存在）
GENE	特征的 /gene 标签（如果存在）
PRODUCT	特征的 /product 标签（如果存在）
EFFECT	该变异的 snpEff 注释结果（.vcf 中的 ANN 标签）

TXT 格式中的列

名称	描述
ID	参考序列 + 样本
LENGTH	参考序列的长度
ALIGNED	已比对的位点数量
UNALIGNED	未比对的位点数量
VARIANT	与参考序列不同的位点数量
HET	杂合位点数量或低质量基因型（用 n 表示，--minqual）
MASKED	参考序列中被屏蔽的位点数量（用 X 表示，--mask）
LOWCOV	本样本中低覆盖度的位点数量（用 N 表示，--mincov）

变异类型

类型	名称	示例
snp	单核苷酸多态性	A => T
mnp	多核苷酸多态性	GC => AT
ins	插入	ATT => AGTT
del	缺失	ACGG => ACG
complex	snp/mnp 的组合	ATTC => GTTA

变异调用器

变异调用由 Freebayes 完成。用户可控制的关键参数如下：

• --mincov - 考虑一个位点所需的最小覆盖 reads 数（默认值=10）
• --minfrac - 与参考序列不同的 reads 所占的最小比例
• --minqual - 最小 VCF 变异调用“质量值”（默认值=100）

使用snippy-vcf_report详细查看变异

如果您在运行Snippy时使用--report选项，程序将自动运行snippy-vcf_report并生成一个snps.report.txt文件，其中针对snps.vcf中的每个SNP都包含如下所示的部分：

~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
>LBB_contig000001:10332 snp A=>T DP=7 Q=66.3052 [7]

         10301     10311     10321     10331     10341     10351     10361
tcttctccgagaagggaatataatttaaaaaaattcttaaataattcccttccctcccgttataaaaattcttcgcttat
........................................T.......................................
,,,,,,  ,,,,,,,,,,,,,,,,,,,,,t,,,,,,,,,,t,,t,,,,,,,,,,,,,,,,g,,,,,,,g,,,,,,,,,t,
,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,, .......T..................A............A.......
.........................A........A.....T...........    .........C..............
.....A.....................C..C........CT.................TA.............
,a,,,,,a,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,t,t,,,g,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
,,,,,ga,,,,,,,c,,,,,,,t,,,,,,,,,,g,,,,,,t,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
                            ............T.C..............G...............G......
                                                    ,,,,,,,g,,,,,,,,g,,,,,,,,,,,
                                                           g,,,,,,,,,,,,,,,,,,,,

如果您希望在运行 Snippy 之后 生成此报告，可以直接运行它：

cd snippydir
snippy-vcf_report --cpus 8 --auto > snps.report.txt

如果您需要可在网页浏览器中查看的 HTML 版本，请使用 --html 选项：

cd snippydir
snippy-vcf_report --html --cpus 16 --auto > snps.report.html

其工作原理是为每个变异位点运行 samtools tview，如果有数千个变异位点，这个过程可能会非常缓慢。建议尽可能将 --cpus 设置为最高值。

选项

• --rgid 将设置 BAM 和 VCF 文件中的读段组（RG）ID（ID）和样本（SM）。如果未提供，将使用 --outdir 文件夹名称同时作为 ID 和 SM。

• --mapqual 是变异检测中可接受的最低比对质量值。BWA MEM 使用 60 表示读段“唯一比对”。

• --basequal 是用于变异检测的核苷酸所需的最低质量值。我们使用 13，对应约 5% 的错误概率。这是 SAMtools 的传统值。

• --maxsoft 是允许比对中软剪辑的碱基数，超过此数值将丢弃该比对。这是为了优先选择全局比对而非局部比对，并会传递给 samclip 工具。

• --mincov 和 --minfrac 用于在现有统计量之外，对比对检测应用硬性阈值。最佳值取决于测序深度和污染率。常用值为 10 和 0.9。

• --targets 接受一个 BED 文件，并仅在那些区域内检测变异。通常不需要，除非你只对特定基因座（如 AMR 基因）的变异感兴趣，但仍在进行全基因组测序（WGS）而非扩增子测序。

• --contigs 允许你从 contig 而非读段中检测 SNPs。它会将 contig 拆分为合成读段，以便在多样本分析中，使这些调用与其他读段样本处于同等基础。

核心 SNP 系统发育

如果你从同一参考序列为多个分离株检测 SNPs，你可以生成“核心 SNPs”的比对序列，用于构建高分辨率的系统发育树（忽略可能的重组）。“核心位点”是在所有样本中都存在的基因组位置。核心位点可以在每个样本中具有相同的核苷酸（“单态性”），或者某些样本可能不同（“多态性”或“变异”）。如果我们忽略“ins”、“del”等变异类型的复杂性，仅使用变异位点，这些就是“核心 SNP 基因组”。

输入要求

• 一组使用相同 --ref 序列的 Snippy 文件夹。

使用 snippy-multi

为简化对同一参考序列运行一组分离株序列（reads 或 contigs）的流程，您可以使用 snippy-multi 脚本。该脚本需要一个如下所示的制表符分隔的输入文件，并且能够处理双端 reads、单端 reads 和组装好的 contigs。

# input.tab = ID R1 [R2]
Isolate1	/path/to/R1.fq.gz	/path/to/R2.fq.gz
Isolate1b	/path/to/R1.fastq.gz	/path/to/R2.fastq.gz
Isolate1c	/path/to/R1.fa		/path/to/R2.fa
# single end reads supported too
Isolate2	/path/to/SE.fq.gz
Isolate2b	/path/to/iontorrent.fastq
# or already assembled contigs if you don't have reads
Isolate3	/path/to/contigs.fa
Isolate3b	/path/to/reference.fna.gz

然后运行此命令以生成输出脚本。 第一个参数应为 input.tab 文件。 其余参数应为任何剩余的共享 snippy 参数。ID 将用于每个分离株的 --outdir。

% snippy-multi input.tab --ref Reference.gbk --cpus 16 > runme.sh

% less runme.sh   # check the script makes sense

% sh ./runme.sh   # leave it running over lunch

最后，它还会运行snippy-core以生成核心基因组SNP比对文件core.*。

输出文件

扩展名	描述
.aln	--aformat格式的核心SNP比对文件（默认FASTA格式）
.full.aln	全基因组SNP比对文件（包含不变位点）
.tab	核心SNP位点的制表符分隔列表，包含等位基因，但无注释
.vcf	多样本VCF文件，包含所有已发现等位基因的基因型GT标签
.txt	比对/核心大小统计信息的制表符分隔列表
.ref.fa	--ref的FASTA版本/副本
.self_mask.bed	如果使用--mask auto，则生成此BED文件。

为什么core.full.aln像一锅“字母汤”？

core.full.aln文件是FASTA格式的多序列比对文件。它包含一条参考序列，以及参与核心基因组计算的每个样本各一条序列。每条序列的长度与参考序列相同。

字符	含义
ATGC	与参考序列相同
atgc	与参考序列不同
-	该样本中覆盖度为零 或 相对于参考序列的缺失
N	该样本中覆盖度低（基于--mincov）
X	参考序列的屏蔽区域（来自--mask）
n	杂合或低质量基因型（在snps.raw.vcf中GT=0/1或QUAL < --minqual）

您可以使用附带的snippy-clean_full_aln工具移除所有“奇怪”字符并将其替换为N。当您需要将文件传递给建树工具或重组去除工具时，这非常有用：

% snippy-clean_full_aln core.full.aln > clean.full.aln
% run_gubbins.py -p gubbins clean.full.aln
% snp-sites -c gubbins.filtered_polymorphic_sites.fasta > clean.core.aln
% FastTree -gtr -nt clean.core.aln > clean.core.tree

选项

• 如果您想对基因组的特定区域进行屏蔽，可以使用 --mask 参数提供一个 BED 文件。这些区域内的任何 SNP 都将被排除。这在处理类似 结核分枝杆菌（M.tuberculosis）的基因组时非常常见，因为其烦人的 PE/PPE/PGRS 重复基因可能导致假阳性结果，或者用于屏蔽噬菌体区域。Snippy 中提供了一个用于 M.tb 的 --mask BED 文件，位于 etc/Mtb_NC_000962.3_mask.bed 文件夹中。该文件来源于 https://gph.niid.go.jp/tgs-tb/ 网站的 XLSX 文件。
• 如果您使用 snippy --cleanup 选项，参考文件将会被删除。这意味着 snippy-core 无法“自动查找”参考序列。在这种情况下，您只需使用 snippy-core --reference REF 来提供 FASTA 格式的参考序列即可。

高级用法

读取过多时提高速度

有时，您的测序深度可能远超 SNP calling 所需。一个常见的问题是，对单个细菌分离株使用整个 MiSeq 测序芯片，产生的 2500 万条 reads 可能导致基因组测序深度高达 2000x。这会使 Snippy 的运行速度远低于必要水平，因为大多数 SNP 在 50-100x 的深度下就能被检测到。如果您确定手头的数据量是实际需求的 10 倍，Snippy 可以对您的 FASTQ 数据进行随机抽样：

# have 1000x depth, only need 100x so sample at 10%
snippy --subsample 0.1  ...
<snip>
Sub-sampling reads at rate 0.1
<snip>

仅在特定区域检测SNP

如果您要寻找特定的SNP（例如参考基因组中特定基因内与AMR相关的SNP），仅在这些区域检测变异可以节省大量时间。只需将感兴趣的区域放入BED文件即可：

snippy --targets sites.bed ...

在contigs之间寻找SNP

有时，你的某个样本仅能获得contigs，而没有相应的FASTQ reads。你仍然可以将这些contigs与Snippy结合使用，以找到相对于参考序列的变异。其原理是将contigs切割成250 bp的单端reads，覆盖度为2 × --mincov的均匀覆盖度。

要使用此功能，无需提供--R1和--R2，而是使用--ctgs选项并指定contigs文件：

% ls
ref.gbk mutant.fasta

% snippy --outdir mut1 --ref ref.gbk --ctgs mut1.fasta
Shredding mut1.fasta into pseudo-reads.
Identified 257 variants.

% snippy --outdir mut2 --ref ref.gbk --ctgs mut2.fasta
Shredding mut2.fasta into pseudo-reads.
Identified 413 variants.

% snippy-core mut1 mut2 
Found 129 core SNPs from 541 variant sites.

% ls
core.aln core.full.aln ...

此输出文件夹与 snippy-core 完全兼容，因此您可以混合基于 FASTQ 和 contig 的 snippy 输出文件夹来生成比对结果。

校正组装错误

从头组装过程会尝试将测序 reads 重构为它们所源自的原始 DNA 序列。这些重构的序列称为 contig（重叠群）或 scaffold（ scaffolds）。由于各种原因，组装得到的 contig 中可能会引入一些小错误，而这些错误并未得到组装过程中所用原始 reads 的支持。

一种常用的策略是将 reads 重新比对回 contig，以检查是否存在差异。这些错误会表现为变异（SNP 和 indel）。如果我们能够“逆转”这些变异，就能“校正”contig，使其与原始 reads 提供的证据相符。显然，如果在 read 比对的执行方式以及变异的筛选上不够谨慎，这种策略可能会出错。

Snippy 能够协助进行这种 contig 校正过程。实际上，它会生成一个 snps.consensus.fa FASTA 文件，该文件是在提供的 ref.fa 输入文件的基础上，应用了 snps.vcf 中发现的变异！

然而，Snippy 并非完美无缺，有时也会发现一些可疑的变异。通常，您可以复制一份 snps.vcf（我们称之为 corrections.vcf），然后删除那些对应于我们不信任的变异的行。例如，在校正 Roche 454 和 PacBio SMRT 测序数据组装的 contig 时，我们主要预期会发现均聚物错误，因此预期会看到更多 ins（插入）类型的变异，而非 snp（单核苷酸多态性）类型的变异。

在这种情况下，您需要使用以下步骤手动运行校正过程：

% cd snippy-outdir
% cp snps.vcf corrections.vcf
% $EDITOR corrections.vcf
% bgzip -c corrections.vcf > corrections.vcf.gz
% tabix -p vcf corrections.vcf.gz
% vcf-consensus corrections.vcf.gz < ref.fa > corrected.fa

您可能希望通过将 corrected.fa 用作 Snippy 重复运行时新的 --ref 来迭代此过程。有时纠正一个错误后，BWA 就能比对之前无法处理的序列，从而发现新的错误。

Snippy 可能并非校正组装结果的最佳工具——您应考虑使用专门的工具，例如 PILON 或 iCorn2，或调整 Quiver 参数（针对 Pacbio 数据）。

未比对的reads

有时您可能会关注那些未与参考基因组比对上的 reads。这些 reads 代表了您样本中潜在的新 DNA 序列，可能具有研究价值。一种标准策略是对未比对的 reads 进行从头组装，以发现这些新的 DNA 元件，它们通常包括质粒等可移动遗传元件。

默认情况下，Snippy不会保留未比对的 reads，即使在 BAM 文件中也不会。如果您希望保留它们，请使用 --unmapped 选项，未比对的 reads 将保存到一个压缩的 FASTQ 文件中：

% snippy --outdir out --unmapped ....

% ls out/
snps.unmapped.fastq.gz ....

信息

名称由来

Snippy 这个名字是以下几个元素的组合： SNP（发音为 "snip"，即单核苷酸多态性）、snappy（意为“迅速的”）以及 Skippy the Bush Kangaroo（用以代表其澳大利亚起源）。

许可证

Snippy 是自由软件，根据 GPL（版本 2） 协议发布。

问题反馈

请将建议和错误报告提交至 Issue Tracker。

运行要求

• perl >= 5.18
• bioperl >= 1.7
• bwa mem >= 0.7.12
• minimap2 >= 2.0
• samtools >= 1.7
• bcftools >= 1.7
• bedtools >= 2.0
• GNU parallel >= 2013xxxx
• freebayes >= 1.1（包括 freebayes、freebayes-parallel、fasta_generate_regions.py）
• vcflib >= 1.0（包括 vcfstreamsort、vcfuniq、vcffirstheader）
• vt >= 0.5
• snpEff >= 4.3
• samclip >= 0.2
• seqtk >= 1.2
• snp-sites >= 2.0
• any2fasta >= 0.4
• wgsim >= 1.8（仅用于测试 - wgsim 命令）

捆绑的二进制文件

对于 Linux（在 Ubuntu 16.04 LTS 上编译）和 macOS（在 High Sierra Brew 上编译），已包含部分二进制文件、JAR 包和脚本。