01 背景
1.1 SAM、BAM文件及SAMtools
随着生物信息学、分子生物学等学科研究的深入,以及人类基因计划的完成,越来越多的人类基因和其他模式生命体的基因被测序。序列比对是处理测序结果的方法,可以发现生物序列之间存在的结构、功能和进化的关系,是生物信息学的基础。 随着这些测序项目的展开,每天都有海量的DNA序列数据产生,DNA序列数据经过序列比对处理,比对结果数据也随之出现。虽然存储设备的快速发展已经在一定程度上缓解了相关数据量急剧膨胀的问题。然而随着比对研究的深入,单纯依靠增加硬件设备已经无法满足DNA比对结果数据量快速增长的需求,存储和使用这些数据的成本也终将增加至无法承担的规模。因此,不论是从存储方面,还是应用方面考虑,序列比对结果的压缩在DNA数据的存储、管理和传输中起到了重要作用。DNA序列数据的压缩目前已经引起了国内外学术界的广泛关注,samtools作为管道pipeline,实现了上述目的!
SAM(Sequence Alignment/Map)和BAM(Binary Alignment/Map)都是存储DNA序列比对信息的文件格式。SAM格式是人类可读的文本格式,而BAM格式是二进制格式,通常用于存储大量的比对数据,因为它在空间上更有效率。
SAM格式文件包含了比对序列的原始信息,例如序列名、比对的质量评分、比对位置、CIGAR字符串(表示比对序列中的插入、删除、错配等信息)、序列和质量信息等。
BAM格式文件是SAM格式文件的二进制表示形式,通常比SAM格式文件更加紧凑,因此在存储和处理大量数据时更为高效。它们使用了索引来加快访问速度,允许在大型基因组上进行快速的搜索和检索。
1.2 SNP的研究需要
在高通量技术和SNP分子标记技术日渐普及的今天,越来越多的研究需要从高通量测序技术产生的reads中获得有效的SNP遗传信息,以进行进一步研究。随着分子实验技术的进一步发展,更多的有效、快速并且高通量的SNP检测方法不停地被研究者开发出来。但目前对于各种SNP calling的方法比较还比较少,究竟什么样的测序策略配合什么种类的SNP calling方法能获得最有效的结果依然不明朗。在大多数研究工作者将该步骤交给测序公司处理分析的背景下,这一点尤为突出。
方法:通过比较目前较为常用的六种SNP calling软件(Varscan, Altas-snp2, GATK, Freebayes, SOAPsnp2和SAMtools)在两种数据集,模拟数据集和真实数据集中calling SNP的表现,对这一问题进行了一定程度上的解读,为研究者对于目前的SNP calling情况进行了解提供了便利。
结果:研究结果显示,这六种SNP calling软件产生的分析结果差异较大:在使用真实数据进行SNP检测的结果中,SOAPsnp能够检测出最多的SNP,Freebays和Atlas-snp2检测出的SNP数量较少;是否去除低质量碱基对这款软件具有一定影响,去除低质量碱基能够检测出更多的SNP;这六种软件calling出的SNP一致性较差,db SNP与Non-db SNP相比通常具有更高的一致性,这表明低覆盖测序数据中被鉴定的Non-db SNP的质量较低。在使用模拟数据进行SNP检测的结果中,测序深度对SNP检测具有显著影响,而测序的reads长度则对SNP检测没有明显影响;如果仅从实验结果看,250bp的测序策略结果最好;在5×下,除了Varscan和Altas-snp2以外,都能得到较好的结果,其余4个软件的SNP检出率都在85%以上;如果测序深度达到7×,基本上除了Varscan,其余软件都能得到很好的结果,检出率都在90%以上。并且所有的SNP检测软件的SNP检测正确率都在97%以上,完全满足科研工作中的需求。
结论:总的来说,本次研究推荐的测序策略为7X测序深度,250bp测序reads长度,在获得数据后的SNP calling过程中,本次研究推荐使用多种SNP calling软件组合的方法,以便获得最优的结果,方便后续的对基因组的分析。
02 参考
https://github.com/samtools/samtools/ #参看官网
03 安装
#下载
wget -c https://github.com/samtools/samtools/releases/download/1.17/samtools-1.17.tar.bz2
tar jxvf samtools-1.17.tar.bz2 #特殊解压格式
mkdir samtoolslocation #创建安装目录
cd samtools-1.17
./configure --prefix=./samtoolslocation #安装在本目录,并创建新文件夹,这里也可自定义绝对路径新文件夹
make
make install
#因为是胶水,所以怕环境冲突,最后安装
conda install -c bioconda samtools
04 使用
./samtools
程序:samtools(用于SAM格式的比对工具)
版本:1.17(使用htslib 1.17)
用法:samtools <命令> [选项]
命令:
-- 索引
dict 创建序列字典文件
faidx 索引/提取FASTA
fqidx 索引/提取FASTQ
index 索引比对结果
-- 编辑
calmd 重新计算MD/NM标签和'='基
fixmate 修复mate信息
reheader 替换BAM头文件
targetcut 切割fosmid区域(仅用于fosmid池)
addreplacerg 添加或替换RG标签
markdup 标记重复
ampliconclip 从读取的末端剪切寡核苷酸
-- 文件操作
collate 按名称对比对结果进行洗牌和分组
cat 连接BAM文件
consensus 生成一致的Pileup/FASTA/FASTQ
merge 合并排序后的比对结果
mpileup 多路比对结果
sort 对比对结果文件进行排序
split 按read group分割文件
quickcheck 快速检查SAM/BAM/CRAM文件是否完整
fastq 将BAM转换为FASTQ
fasta 将BAM转换为FASTA
import 将FASTA或FASTQ文件转换为SAM/BAM/CRAM
reference 从比对数据生成参考序列
reset 恢复读取的比对变化
-- 统计
bedcov BED区域的读取深度
coverage 比对深度和百分比覆盖率
depth 计算深度
flagstat 简单统计
idxstats BAM索引统计
cram-size 列出CRAM内容ID和数据系列大小
phase 对杂合子进行相位分析
stats 生成统计信息(之前的bamcheck)
ampliconstats 生成引物特异性统计信息
-- 查看
flags 解释BAM标志
head 头文件查看器
tview 文本比对查看器
view SAM<->BAM<->CRAM转换
depad 将带填充的BAM转换为不带填充的BAM
samples 列出一组SAM/BAM/CRAM文件中的样本
-- 其他
help [cmd] 显示此帮助消息或[cmd]的帮助 钥匙在这里!!!
version 详细的版本信息
05 常用命令行
一般为许多pipeline的依赖或者内置程序。一般多用view查看sam,bam文件或实现两者格式互换。
#将sam转换为bam
./samtools view -S in.sam -b > out.bam
#将bam转换为sam
./samtools view -h in.bam -o out.sam
用法:samtools view [选项] <in.bam>|<in.sam>|<in.cram> [区域...]
输出选项:
-b, --bam 输出BAM格式
-C, --cram 输出CRAM格式(需要 -T)
-1, --fast 使用快速BAM压缩(并默认为 --bam)
-u, --uncompressed 未压缩的BAM输出(并默认为 --bam)
-h, --with-header 在SAM输出中包含头文件
-H, --header-only 仅打印SAM头文件(无比对记录)
--no-header 仅打印SAM比对记录 [默认]
-c, --count 仅打印匹配记录的计数
-o, --output FILE 将输出写入文件 FILE [标准输出]
-U, --unoutput FILE, --output-unselected FILE
将未通过过滤器选择的读取输出到文件 FILE
-p, --unmap 对未被选择的读取设置为UNMAP标志,然后写入输出文件。
-P, --fetch-pairs 即使在区域外也检索完整的配对
输入选项:
-t, --fai-reference FILE 列出参考名称和长度的文件
-M, --use-index 使用索引和多区域迭代器进行区域选择
--region[s]-file FILE 使用索引仅包含与文件 FILE 重叠的读取
-X, --customized-index 在 <in.bam> 后期望额外的索引文件参数
过滤选项(仅包含输出中的读取...):
-L, --target[s]-file FILE ...与文件 FILE 中的区域(BED)重叠
-r, --read-group STR ...位于读取组 STR 中
-R, --read-group-file FILE ...位于文件 FILE 中列出的读取组中
-N, --qname-file FILE ...其读取名称在文件 FILE 中列出
-d, --tag STR1[:STR2] ...具有标签 STR1(及其关联值 STR2)
-D, --tag-file STR:FILE ...具有值列在 FILE 中的标签 STR
-q, --min-MQ INT ...具有 >= INT 的比对质量
-l, --library STR ...位于库 STR 中
-m, --min-qlen INT ...覆盖 >= INT 查询碱基数(通过CIGAR测量)
-e, --expr STR ...与过滤表达式 STR 匹配
-f, --require-flags FLAG ...具有所有 FLAG 标志
-F, --excl[ude]-flags FLAG ...不具有任何 FLAG 标志
--rf, --incl-flags, --include-flags FLAG
...具有某些 FLAG 标志
-G FLAG 排除具有所有 FLAG 标志的读取
--subsample FLOAT 仅保留 FLOAT 比例的模板/读取对
--subsample-seed INT 影响哪些读取在子采样中被保留 [0]
-s INT.FRAC 等同于 --subsample 0.FRAC --subsample-seed INT
处理选项:
--add-flags FLAG 将 FLAG 添加到读取
--remove-flags FLAG 从读取中删除 FLAG
-x, --remove-tag STR
要删除的逗号分隔的读取标签(可重复)[空]
--keep-tag STR
要保留的逗号分隔的读取标签(可重复)[空]。
等同于 "-x ^STR"
-B, --remove-B 折叠向后的CIGAR操作
-z, --sanitize FLAGS 在记录上执行检查和修复。
FLAGS 用逗号分隔(参见手册)。[关闭]
常规选项:
-?, --help 打印长帮助,包括区域指定的注释
-S 被忽略(输入格式会自动检测)
--no-PG 不添加PG行
--input-fmt-option OPT[=VAL]
以 OPTION 或 OPTION=VALUE 形式指定单个输入文件格式选项
-O, --output-fmt FORMAT[,OPT[=VAL]]...
指定输出格式(SAM、BAM、CRAM)
--output-fmt-option OPT[=VAL]
以 OPTION 或 OPTION=VALUE 形式指定单个输出文件格式选项
-T, --reference FILE
参考序列FASTA文件 [空]
-@, --threads INT
使用的附加线程数 [0]
--write-index
自动为输出文件建立索引 [关闭]
--verbosity INT
设置详细程度
The original samtools package has been split into three separate but tightly coordinated projects:
- htslib: C-library for handling high-throughput sequencing data
- samtools: mpileup and other tools for handling SAM, BAM, CRAM
- bcftools: calling and other tools for handling VCF, BCF
See also samtools · GitHub
他也是trinity的依赖项之一,详见Trinity安装与使用-Trinity-v2.15.1(bioinfomatics tools-006)
06 参考文献
Twelve years of SAMtools and BCFtools
Petr Danecek, James K Bonfield, Jennifer Liddle, John Marshall, Valeriu Ohan, Martin O Pollard, Andrew Whitwham, Thomas Keane, Shane A McCarthy, Robert M Davies, Heng Li
GigaScience, Volume 10, Issue 2, February 2021, giab008, Twelve years of SAMtools and BCFtools | GigaScience | Oxford Academic
陈灿. DNA序列比对结果的存储与压缩[D].复旦大学,2015.
何熠. 基于二代测序技术对SNP检测软件的比较研究[D].石河子大学,2020.