我用到的Samtools介绍

Posted 2023-05-09

参考技术A 记录一下我用到的samtools的用法。

samtools的说明文档： http://samtools.sourceforge.net/samtools.shtml

bam文件优点：bam文件为二进制文件，占用的磁盘空间比sam文本文件小；利用bam二进制文件的运算速度快。

首先需要意识到的是samtools是一个非常强大的工具，想要熟练的使用它，还需要不断的摸索。

samtools的用法

（1）View

samtools view -bS abc.sam > abc.bam    #将sam文件转换为bam文件

 参数：

-b bam 输出bam

-S sam 输入sam

-@ 线程

在比对完成的sam文件中，包含着mapped reads 和unmapped reads

$ samtools view -bF 4  abc.bam > abc.F.bam       #提取没有比对到参考序列上的比对结果，步包含标签

$ samtools view -bF 12 abc.bam > abc.F12.bam   #提取paired reads中两条reads都比对到参考序列上的比对结果，只需要把两个4+8的值12作为过滤参数即可

$ samtools view -bf 4 abc.bam > abc.f.bam    #提取没有比对到参考序列上的比对结果，包含标签

$ samtools view abc.bam scaffold1 > scaffold1.sam     #提取bam文件中比对到caffold1上的比对结果，并保存到sam文件格式

$ samtools view abc.bam scaffold1:30000-100000 $gt; scaffold1_30k-100k.sam    #提取scaffold1上能比对到30k到100k区域的比对结果

$ samtools view -T genome.fasta -h scaffold1.sam > scaffold1.h.sam    #根据fasta文件，将 header 加入到 sam 或 bam 文件中

samtools的view不就可以进行格式转换，还可以进行数据的提取

例：提取1号染色体上1234~123456区域的以对read

samtools view SRR3589957_sorted.bam chr1:1234-123456| head

 samtools view SRR3589957_sorted.bam chr1:1234-123456 > sub.bam 

使FLAG更具可读性

samtools view -X sample.sorted.bam | head -n 5

计算总的比对数量

samtools view sample.sorted.bam | wc -l

显示标题，-H选项

samtools view -H sample.sorted.bam

将bam文件转换为sam文件

samtools view -h abc.bam > abc.sam

（2）Sort

samtools sort对bam文件进行排序，不能对sam文件进行排序。

以leftmost coordinates的方式对比对结果进行排序，或者使用-n参数以read名称进行排序。将会添加适当的@HD-SO排序顺序标头标签或者如果有必要的话，将会更新现存的一个排序顺序标头标签。sort命令的输出默认是标准输出写入，或者使用-o参数时，指定bam文件输出名。sort命令还会在内存不足时创建临时文件tmpprefix.%d.bam。

也就是说：samtools的排序方式有两种（常用）

默认方式，按照染色体的位置进行排序

samtools sort test.bam default

参数-n则是根据read名进行排序。

samtools sort -n test.bam sort_left

usage: samtools sort [-l level] [-m maxMem] [-o out.bam] [-O format] [-n] [-T tmpprefix] [-@ threads] [in.sam|in.bam|in.cram]

例如：samtools sort abc.bam abc.sort

samtools sort -O bam -@ 2 SRR1909070.bam -o SRR1909070.sorted.bam

RNA-seq 的数据比对结果 BAM 文件使用 samtools 进行 sort 之后文件压缩比例变化会比DNA-seq 更甚。另外，samtools 对 BAM 文件进行排序之后那些没有比对上的 reads 会被放在文件的末尾。

 参数：

-l INT 设置输出文件压缩等级。0-9，0是不压缩，9是压缩等级最高。不设置此参数时，使用默认压缩等级；

-m INT 设置每个线程运行时的内存大小，可以使用K，M和G表示内存大小。

-n 设定排序方式按short reads的ID排序。默认下是按序列在fasta文件中的顺序（即header）和序列从左往右的位点排序。

-o FILE 设置最终排序后的输出文件名；

-O FORMAT 设置最终输出的文件格式，可以是bam，sam或者cram，默认为bam；

-T PREFIX 设置临时文件的前缀；

-@ INT 设置排序和压缩是的线程数量，默认是单线程。

（3）index

samtools index 建立索引，在建立索引之前应该先对bam文件进行排序。必须对bam文件进行 默认情况下的排序后 ，才能进行index。否则会报错。

建立索引后将产生后缀为.bai的文件，用于快速的随机处理。很多情况下需要有bai文件的存在，特别是显示序列比对情况下。比如samtool的tview命令就需要；gbrowse2显示reads的比对图形的时候也需要。

samtools index abc.sort.bam

如果想要建立索引的，具体可以看看比对的内部的算法，链接具体是怎么建立索引的

建立索引的目的应该是为了提高比对的效率

以下两种命令结果一样

$ samtools index abc.sort.bam

$ samtools index abc.sort.bam abc.sort.bam.bai

（4）flagstat

samtools flagstat  给出BAM文件的比对结果

samtools flagstat [options] <in.bam>

-@ 线程

-O FORMAT 设置最终输出的文件格式，可以是txt，json或者tsv，默认为json，tsv；

samtools flagstat输出结果解释:

11945742 + 0 in total (QC-passed reads + QC-failed reads)

#总共的reads数

0 + 0 duplicates

7536364 + 0 mapped (63.09%:-nan%)

#总体上reads的匹配率

11945742 + 0 paired in sequencing

#有多少reads是属于paired reads

5972871 + 0 read1

#reads1中的reads数

5972871 + 0 read2

#reads2中的reads数

6412042 + 0 properly paired (53.68%:-nan%)

#完美匹配的reads数：比对到同一条参考序列，并且两条reads之间的距离符合设置的阈值

6899708 + 0 with itself and mate mapped

#paired reads中两条都比对到参考序列上的reads数

636656 + 0 singletons (5.33%:-nan%)

#单独一条匹配到参考序列上的reads数，和上一个相加，则是总的匹配上的reads数。

469868 + 0 with mate mapped to a different chr

#paired reads中两条分别比对到两条不同的参考序列的reads数

243047 + 0 with mate mapped to a different chr (mapQ>=5)

#paired reads中两条分别比对到两条不同的参考序列的reads数,并且其中比对质量>=5的reads的数量

（5）depth

得到每个碱基位点的测序深度,并输出到标准输出。

usage: samtools depth [options] in.bam [in.bam ...]

注意 ：做depth之前必须做samtools index；

示例：

samtools depth in.bam  >  out.depth.txt

注意： in.bam 必须经过了排序。

（6）samtools rmdup

NGS上机测序前需要进行PCR一步，使一个模板扩增出一簇，从而在上机测序的时候表现出为1个点，即一个reads。若一个模板扩增出了多簇，结果得到了多个reads，这些reads的坐标(coordinates)是相近的。在进行了reads比对后需要将这些由PCRduplicates获得的reads去掉，并只保留最高比对质量的read。使用rmdup命令即可完成.

Usage:

samtools rmdup[-sS]

-s对single-end reads。默认情况下，只对paired-endreads

-S将Paired-endreads作为single-endreads处理。

$samtools rmdup input.sorted.bam output.bam

（7）mpileup

samtools还有个非常重要的命令mpileup，以前为pileup。该命令用于生成bcf文件，再使用bcftools进行SNP和Indel的分析。bcftools是samtool中附带的软件，在samtools的安装文件夹中可以找到。

最常用的参数有2个：

-f来输入有索引文件的fasta参考序列；

-g输出到bcf格式。用法和最简单的例子如下

Usage:samtoolsmpileup[-EBug][-CcapQcoef][-rreg][-fin.fa][-llist][-McapMapQ][-QminBaseQ][-qminMapQ]in.bam[in2.bam[...]]

$samtoolsmpileup-fgenome.fastaabc.bam>abc.txt

$samtoolsmpileup-gSDfgenome.fastaabc.bam>abc.bcf

$samtoolsmpileup-guSDfgenome.fastaabc.bam|\bcftoolsview-cvNg->abc.vcf

mpileup不使用-u或-g参数时，则不生成二进制的bcf文件，而生成一个文本文件(输出到标准输出)。该文本文件统计了参考序列中每个碱基位点的比对情况；该文件每一行代表了参考序列中某一个碱基位点的比对结果。比如：

（8）faidx

对fasta文件建立索引，比如基因组的文件，生成的索引文件以.fai后缀结尾。该命令也能依据索引文件快速提取fasta文件中的某一条（子）序列

Usage: samtools faidx <in.bam> [ [...]]

对基因组文件建立索引

$ samtools faidx genome.fasta

生成了索引文件genome.fasta.fai,是一个文本文件，分成了5列。

第一列是子序列的名称；

第二列是子序列的长度；

第三列是序列所在的位置，因为该数字从上往下逐渐变大，最后的数字是genome.fasta文件的大小；

第4和5列不知是啥意思。于是通过此文件，可以定

位子序列在fasta文件在磁盘上的存放位置，直接快速调出子序列。

由于有索引文件，可以使用以下命令很快从基因组中提取到fasta格式的子序列

$ samtools faidx genome.fasta scffold_10 > scaffold_10.fasta

拓展：bcftools软件

bcftools和samtools类似，用于处理vcf(variant call format)文件和bcf(binary call format)文件。前者为文本文件，后者为其二进制文件。

bcftools使用简单，最主要的命令是view命令，其次还有index和cat等命令。index和cat命令和samtools中类似。此处主讲使用view命令来进行SNP和Indel calling。该命令的使用方法和例子为：

$ bcftools view -cvNg abc.bcf > snp_indel.vcf

生成的结果文件为vcf格式，有10列，分别是：1 参考序列名；2 varianti所在的left-most位置；3 variant的ID（默认未设置，用’.'表示)；4 参考序列的allele；5 variant的allele(有多个alleles，则用’,'分隔);6 variant/reference QUALity;7 FILTers applied;8 variant的信息，使用分号隔开；9 FORMAT of the genotype fields, separated by colon (optional)； 10 SAMPLE genotypes and per-sample information (optional)。

参考链接：

原文链接：https://blog.csdn.net/u013553061/article/details/53179945

https://www.cnblogs.com/emanlee/p/4316581.html

http://events.jianshu.io/p/794d82bccf6c

http://blog.sina.com.cn/s/blog_13de3725c0102v7rd.html

https://www.cnblogs.com/shuaihe/articles/6802246.html

无root权限下利用miniconda2成功安装samtools1.5

刚开始学生信，首先碰到的问题就是在没有root权限的前提下服务器上安装各种生信软件，一开始就碰到一个钉子。搜索了网络上各大论坛关于samtools安装教程，从中学到了很多，在此表示对各位战友的详细解说表示感谢。在此介绍一个比较简单实用的samtools安装方法。首先，安装miniconda2，安装过程可参考博文“BioConda --生信工作者的福音”，网址：http://blog.sciencenet.cn/blog-2970729-1074395.html；其次，运行命令conda install samtools=1.5,即可。提示：不要安装最新版本，往往会有问题。