更快的处理bam数据—Sambamba

生信菜鸟团

发布于 2023-12-28 13:40:13

3.5K00

代码可运行

文章被收录于专栏：生信菜鸟团生信菜鸟团

运行总次数：0

代码可运行

工欲善其事必先利其器

1sambamba

sambamba 主要是由Artem Tarasov开发的一款高效的生物信息学工具，主要用于处理大规模的测序数据，尤其是针对SAM/BAM格式的文件。这个软件的设计目的是为了提供比现有工具(samtools)更快的性能，特别是在多核处理器系统上，它利用多核处理并显著缩短处理时间。其具有以下特性：

多线程处理：Sambamba能够利用多核处理器的优势，通过并行处理来加快数据的读取、排序和索引等操作。
兼容性：Sambamba设计上与常用的生信工具如Samtools兼容，使得它可以容易地集成到现有的分析流程中。
主要编程语言：D语言
支持平台：支持Linux和MacOS，不支持Windows

2发表文章

题目：Sambamba: fast processing of NGS alignment formats 期刊：Bioinformatics 日期：2015/6/15 作者&单位：Artem Tarasov DOI：https://doi.org/10.1093/bioinformatics/btv098

性能测试对比

从图中可以看到，受限于服务IO的读取速度，sambamba虽然支持多线程，但也并不是调用线程越多越快，建议一般使用4个线程即可。

3简要用途

sambamba 子命令

4如何安装

conda安装

#codna create -n wes #先创建小环境，如果已经创建，可以忽略 
conda activate wes 
conda install -y sambamba

预编译的二进制版安装

sambamba同时提供预编译好的二进制版安装包，下载解压，可直接运行，也非常方便。

wget -c https://github.com/biod/sambamba/releases/download/v1.0.1/sambamba-1.0.1-linux-amd64-static.gz
##解压到指定文件夹
gunzip  -c sambamba-1.0.1-linux-amd64-static.gz > ~/software/bin/sambamba-1.0.1
##给文件添加可执行权限
chmod +x sambamba-1.0.1

## 检查是否安装成功
./sambamba-1.0.1 --help

5子命令用法及参数

sambamba 最经常用到的功能应该就是标记重复

markdup — 标记去重

识别并标记（默认）或移除在测序数据中出现的重复reads。重复reads通常是测序或样本准备过程中的 PCR 扩增产生的，它们可能会影响后续变异检测和其他生物信息学分析的准确性。在判断一个读取是否为重复时，采用的是与 Picard 工具相同的标准。这些标准通常包括比对的起始位置、方向和库ID等因素。如果两个或多个读取具有相同的起始位置和方向，并且来自同一个库，它们通常会被认为是重复的。

sambamba markdup -t 4 --tmpdir ~/test d0.bam  ~/test/d0_mkdup.bam

##其余参数
-r: #移除重复的reads，而不仅仅是标记它们。【直接从数据中清除被识别为重复的reads】
-t: #设定使用的线程数量
-l: #指定结果文件的压缩级别，范围从 0（无压缩）到 9（最大压缩）
-p: #在标准错误输出 (STDERR) 中显示进度条。这有助于监控长时间运行操作的进度
--tmpdir=TMPDIR: #指定临时文件的存储目录
--sort-buffer-size=SORT_BUFFER_SIZE: #设定用于排序过程的总内存量；默认值为2048M，增加它将减少创建的临时文件数量以及主线程中花费的时间
--io-buffer-size=BUFFER_SIZE: #在第二遍读取和写入 BAM 时，使用两个 BUFFER_SIZE 的缓冲区（默认是128M）。这会影响数据读写的效率和速度

sort — 排序

用于对 BAM 文件进行排序，这是许多生物信息学分析的关键步骤。默认会同时输出排序后的文件 .sorted.bam，以及排序后的索引文件 .sorted.bam.bai

sambamba sort -t 4 d0.bam --tmpdir ~/test 

## 其余参数
-m: 指定所有线程的总内存限制，默认为2GB。这个参数可以控制 `sambamba sort` 在排序过程中使用的内存量，以避免耗尽系统资源
--tmpdir=TMPDIR: 指定临时文件的存储目录；默认是系统的临时文件目录
-o: 指定输出文件名（可直接定义输出文件的位置和名称）；如果未提供，则结果写入一个以 `.sorted.bam` 为扩展名的文件
-n: 按read名而不是坐标排序（字典顺序）。这种排序对于某些特定分析可能更有用，尤其是当read名中的信息对于后续处理很重要时
--sort-picard: 像 Picard 工具一样按 query name 排序。这可以确保与使用 Picard 工具时的兼容性和一致性
-N: 按read name 而不是坐标进行所谓的“natural”排序（如 samtools 中的排序）。这与 `-n` 类似，但排序方法更接近人类直觉地理解数字和字符的组合 (不建议使用，因为此种排序方式可能会与GATK流程不兼容，见：https://cloud.tencent.com/developer/article/2009777)
-M, --match-mates: 在按read name排序时，将同一比对的配对read放在一起。通常用于需要分析或处理配对末端read的情况
-l: 设置排序后的 BAM 文件的压缩级别，从0（无压缩）到9（最大压缩）
-u: 将排序后的 BAM不压缩输出（默认是以压缩级别1写入），在某些情况下这可能更快，但会使用更多的磁盘空间
-p: 在 STDERR 中显示进度条
-t, --nthreads=NTHREADS: 使用指定数量的线程
-F: 仅保留满足 FILTER 条件的read。在排序过程中进行read过滤，仅保留对后续分析有用的数据

index — 建索引

用于为按坐标排序的 BAM 文件创建索引。在运行 sambamba index 之前，BAM 文件必须已经按照参考序列的坐标进行了排序。默认输出与 BAM 文件同名的文件，但扩展名为 .bai

sambamba index -t 4  d0.bam 

##其余参数
-t: #设定使用的线程数量
-l: #指定结果文件的压缩级别，范围从 0（无压缩）到 9（最大压缩）
-p: #在标准错误输出 (STDERR) 中显示进度条。这有助于监控长时间运行操作的进度
-c: #检查bins（一种存储数据位置信息的结构）是否设置正确；这是一种完整性检查，确保索引的准确性和有效性
-F, --fasta-input: #指定输入文件为 FASTA 格式。默认情况下，`sambamba index` 输入是 BAM 文件。如果你需要为FASTA 文件创建索引（例如，基因组参考序列），则需要使用此选项

view — 查看、过滤

主要用于高效地过滤 BAM 文件以及访问 SAM 头部信息和参考序列信息。默认情况下，sambamba view 期望输入文件为 BAM 格式。要使用 SAM 格式的文件，你需要显示指定 -S 或 --sam-input 参数，因为sambamba view 不会尝试从文件扩展名猜测文件格式。如果遇到不符合 SAM/BAM 规范的标签或字段，sambamba view 会跳过这些标签并将无效字段设置为默认值。这意味着即使源文件有些小错误或不规范的地方，工具也能继续运行，但可能会忽略或修改某些数据。

## 查看bam文件，不显示heads 信息
sambamba view d0.bam|less -SN
## 查看bam文件，显示heads 信息
sambamba view  -h d0.bam|less -SN
## 仅输出heads信息
sambamba view  -H  d0.bam >head_info.txt
## 计算出映射质量（mapping quality）不低于50的read数量
sambamba view -c -F "mapping_quality >= 50" d0.bam
## 统计在参考基因组第一条染色体的100到200区域内，所有正确配对并有重叠部分的read的数量
sambamba view -c -F  "proper_pair" d0.bam chr1:100-200


-F: #对比对结果进行自定义过滤。你可以根据需要指定各种过滤条件，如特定的比对质量、标记或其它特征   
-f: #指定输出的格式，默认为 SAM。也可以选择 BAM、JSON 或解压缩的 BAM（unpack）
-h: #在reads之前打印头部信息（对于 BAM 输出总是这样做）。这对于保持文件的上下文信息很有用
-H: #仅将头部信息输出到标准输出（如果格式为 BAM，则头部信息以 SAM 格式输出）。这对于获取文件的元数据很有用
-I: #以 JSON 格式输出参考序列的名称和长度到标准输出。这有助于快速检索关于参考序列的信息
-L：#输出与 BED 文件中的某些区域重叠的读取。这对于关注特定基因区域或目标序列非常有用
-c : #输出匹配记录的计数到标准输出(忽略 hHI参数)。这可以快速提供过滤后的比对数量
-v: #仅输出有效的比对。如果你需要对比对进行完整性验证，可以使用选项。这将更严格地检查数据的有效性，确保所有比对都符合预期的质量和格式标准
-S: #指定输入格式为 SAM
-T: #指定写入时使用的参考文件（默认为 空）。可以确保比对正确对应于参考序列
-p : #在 STDERR 中显示进度条（仅当没有指定区域且文件为 BAM 时有效）。这有助于监控长时间运行的操作的进度
-l : #指定压缩级别（从0到9，仅对 BAM 输出有效）
-o : #指定输出文件名，可以直接定义输出文件的位置和名称
-t : #设置使用的最大线程数。建议4个线程即可
-s: #对读取（读取对）进行抽样。这是减少数据量以进行快速分析或测试的一种方法
--subsampling-seed=SEED : #设置抽样的种子。这可以确保了抽样的可重复性

merge —合并

主要用途是将多个排序过的 BAM 文件合并成一个单一的 BAM 文件。所有输入文件必须具有相同的排序顺序（例如，都是按坐标或按read name 排序）。就像 Picard 等合并工具一样，SAM 文件的 headers（包含关于参考序列、程序参数等的元数据）会自动合并。这意味着来自所有输入文件的重要信息都会被保留并整合到最终合并的文件中，确保了文件的完整性和可用性

##合并2个bam
sambamba merge -t 4 out_merge.bam d0.sorted.bam chr3_200k-300k.sorted.bam 

-t: #设定使用的线程数量
-l: #指定结果文件的压缩级别，范围从 0（无压缩）到 9（最大压缩）
-p: #在标准错误输出 (STDERR) 中显示进度条。这有助于监控长时间运行操作的进度
-H: #将合并后的 head 信息以SAM格式输出到标准输出（stdout），其他选项将被忽略；主要用于调试，用户可以查看和验证合并后的头部信息，确保所有必要的信息都被正确地合并  
-F, --filter=FILTER: #仅保留满足 FILTER 条件的read；在合并过程中对read进行过滤，仅保留对后续分析有用的数据

slice — 切片

用于从BAM 或 FASTA 文件中提取指定区域的reads 。虽然 sambamba view 也可以用来提取指定区域的read，但 sambamba slice 在这个任务上通常会更快

区域以标准形式给出，即 ref:beg-end，其中 ref 是参考序列的名字，beg 和 end 是区域的开始和结束位置。这允许精确指定想要提取的序列区域。
对于那些没有参考序列的读取，可以使用特殊的区域 '*' 来指定。

## 提取d0.sorted.bam 中 chr3 的 200k 到 300k 区域内的所有reads
sambamba slice -o chr3_200k-300k.bam d0.sorted.bam chr3:200000-300000

-o: #指定输出的 BAM 或 FASTA 文件名。如果不指定，输出默认是到标准输出（STDOUT）
-L, --regions=FILENAME: #仅输出与 BED 文件中的某些区域重叠的读取。BED 文件是一种常用的格式，用于指定一系列的基因组区域。该参数允许用户基于复杂的区域列表进行操作，而不用手动指定每个区域
-F, --fasta-input: #显示指定输入文件为 FASTA 格式

flagstat — 统计

从read flags 中提取和输出统计信息

统计信息第一行是过质量控制（QC-passed）和未通过质量控制（QC-failed）的 read 数量，然后分别对通过和未通过的read进行统计

 sambamba flagstat -b d0.bam > d0_stat.csv

-l: #指定结果文件的压缩级别，范围从 0（无压缩）到 9（最大压缩）
-p: #在标准错误输出 (STDERR) 中显示进度条
-b: #以 CSV 格式输出结果

统计信息

depth — 覆盖度统计

用于计算 BAM 文件中指定区域覆盖深度

其有三种模式：base、region和window，每种模式都有其特定的应用场景和参数

共同参数

-F, --filter=FILTER: 设置对比对的自定义过滤条件。默认值是 'mapping_quality > 0 and not duplicate and not failed_quality_control'，这意味着只计算那些映射质量大于0、非重复、质量控制通过的read
-o: 指定输出文件名，默认输出到标准输出
-t: 设定线程
-c, --min-coverage=MINCOVERAGE: 设置输出的最小平均覆盖深度，默认为0（region/window模式）或1（base模式）。只有平均覆盖度达到这个阈值的区域才会被报告
-C, --max-coverage=MAXCOVERAGE: 设置输出的最大平均覆盖深度。这有助于识别和排除异常高覆盖的区域
-q, --min-base-quality=QUAL: 不计算低于此质量值的碱基。这有助于提高覆盖深度计算的准确性
--combined: 输出所有样本的组合统计。通常用于比较多个样本的覆盖深度
-a, --annotate: 添加额外的列来标记是否满足给定的标准，而不是跳过不满足条件的记录
-m, --fix-mate-overlaps: 检测配对读取的重叠部分，并在每个碱基的基础上处理它们；这有助于更准确地计算覆盖度

base模式特定选项

-L, --regions=FILENAME|REGION: 指定感兴趣区域的列表或单个区域的形式（例如 chr:beg-end）。通常用于分析特定基因或区域

region模式特定选项

-L, --regions=FILENAME|REGION: （必需）指定感兴趣区域的列表或单个区域的形式
-T, --cov-threshold=COVTHRESHOLD: 提供一个或多个覆盖度阈值，对于每个阈值，会添加一个额外的列，显示区域中覆盖度超过此值的碱基的百分比

window模式特定选项

-w, --window-size=WINDOWSIZE: 窗口的宽度，以碱基对（bp）为单位（必需）；这定义了计算覆盖度的窗口大小
--overlap=OVERLAP: 连续窗口之间的重叠，以碱基对（bp）为单位（默认是0）；这可以帮助平滑覆盖度的变化
-T, --cov-threshold=COVTHRESHOLD: 与 'region' 子命令中的含义相同，表示覆盖度的阈值

subsample —抽样

用于对 BAM 文件进行抽样

注：目前测试发现该功能柜并没有设置随机数种子的参数，且非常耗费内存

## 随机抽取大约10%的reads 生成一个新BAM 文件
sambamba-1.0.1 subsample --type=fasthash --max-cov=10  -o subsampled.bam d0.bam 

--type [fasthash]: #选择用于抽样的算法。`fasthash` 是一个选项，它提供了一种快速的抽样方法；默认情况下，不使用特定的算法
--max-cov [depth]: #设置所需的最大覆盖深度（approx）。这个参数允许你控制输出样本的覆盖深度，以便在保持足够数据的同时减少数据量
-o: #设置输出文件名；默认情况下，输出是到标准输出（STDOUT）
-r: #从输出中移除过度采样的read；通过移除那些超过指定深度的read，可以帮助确保抽样结果更接近所需的覆盖度

参考