本文目录一览：
- 1、sam文件转换为bam文件——SAMtools
- 2、samtools常用命令详解
- 3、利用samtools将sam格式的文件与bam格式的文件进行相互转换
- 4、fasta、sam、bam相互转换
- 5、我用到的Samtools介绍
- 6、samtools使用大全
sam文件转换为bam文件——SAMtools
在转录组数据与参考基因组进行比对后，得到sam文件，后续分析需要将sam转换为bam，这里用到的工具是SAMtools。
序列比对 —— Hisat2 - (jianshu.com)
SAMtools的主要功能是读取、写出、编辑、查看SAM/BAM/CRAM格式的数据文件。
SAMtools官方网站：
链接：
安装：
这里看到基础命令主要分为六个模块，分别为索引、编辑、文件操作、统计、视图以及其他。这里想要SAM格式转为BAM格式，主要用到的是Viewing模块中的view。
说明书：
在老版本1.9的samtools中，需要用-s 指定sam文件，1.14中不需要指定。
samtools view命令完成sam转为bam。
-@8：8个线程
-b：输出格式bam文件 ：输出文件名
samtools常用命令详解
samtools是一个用于操作sam和bam文件的工具集合。
view命令的主要功能是： 将输入文件转换成输出文件，通常是将比对后的sam文件转换为bam文件 ，然后对bam文件进行各种操作，比如数据的排序(和提取(这些操作是对bam文件进行的，因而当输入为sam文件的时候，不能进行该操作)。-bam文件转换成是sam文件格式
view命令中，对sam文件头部的输入( -t 或 -T ）和输出( -h )是单独的一些参数来控制的。
默认情况下不加 region，则是输出所有的 region.
Options:
-b output BAM
默认下输出是 SAM 格式文件，该参数设置输出 BAM 格式
-h print header for the SAM output
默认下输出的 sam 格式文件不带 header，该参数设定输出sam文件时带 header 信息
-H print header only (no alignments)
-S input is SAM
默认下输入是 BAM 文件，若是输入是 SAM 文件，则最好加该参数，否则有时候会报错。
-u uncompressed BAM output (force -b)
该参数的使用需要有-b参数，能节约时间，但是需要更多磁盘空间。
-c Instead of printing the alignments, only count them and print the total number. All filter options, such as ‘-f’, ‘-F’ and ‘-q’ , are taken into account.-bam文件转换成是sam文件格式
-1 fast compression (force -b)
-x output FLAG in HEX (samtools-C specific)
-X output FLAG in string (samtools-C specific)
-c print only the count of matching records
-L FILE output alignments overlapping the input BED FILE [null]
-t FILE list of reference names and lengths (force -S) [null]
使用一个list文件来作为header的输入
-T FILE reference sequence file (force -S) [null]
使用序列fasta文件作为header的输入
-o FILE output file name [stdout]
-R FILE list of read groups to be outputted [null]
-f INT required flag, 0 for unset [0]
-F INT filtering flag, 0 for unset [0]
Skip alignments with bits present in INT [0]
数字4代表该序列没有比对到参考序列上
数字8代表该序列的mate序列没有比对到参考序列上
-q INT minimum mapping quality [0]
-l STR only output reads in library STR [null]
-r STR only output reads in read group STR [null]
-s FLOAT fraction of templates to subsample; integer part as seed [-1]
-? longer help/pre
如果没有指定option和region选项，则会在屏幕中显示sam格式的文件，如下所示。
例子：
（1）将sam文件和bam文件相互转换。注意： 没有header的sam文件不能转换成bam文件。
（2） 提取比对到参考序列上的比对结果 F
（3）提取paired reads中两条reads都比对到参考序列上的比对结果，只需要把两个4+8的值12作为过滤参数即可
（4） 提取没有比对到参考序列上的比对结果 f
（5）提取paired reads中比对到参考基因组上的数据
（6）提取bam文件中比对到caffold1上的比对结果，并保存到sam文件格式
（7）提取scaffold1上能比对到30k到100k区域的比对结果
（8）根据fasta文件，将 header 加入到 sam 或 bam 文件中
（9）sam 和 bam 格式转换
对bam文件进行排序，不能对sam文件进行排序。以leftmost coordinates的方式对比对结果进行排序，或者使用-n参数以read名称进行排序。将会添加适当的@HD-SO排序顺序标头标签或者如果有必要的话，将会更新现存的一个排序顺序标头标签。sort命令的输出默认是标准输出写入，或者使用-o参数时，指定bam文件输出名。sort命令还会在内存不足时创建临时文件tmpprefix.%d.bam-bam文件转换成是sam文件格式
Options:
-m 设置每个线程运行时的内存大小，可以使用K，M和G表示内存大小。默认下是 500,000,000 即500M。对于处理大数据时，如果内存够用，则设置大点的值，以节约时间。
-n 设置按照read名称进行排序。默认下是按序列在fasta文件中的顺序（即header）和序列从左往右的位点排序。
-l INT 设置输出文件压缩等级。0-9，0是不压缩，9是压缩等级最高。不设置此参数时，使用默认压缩等级；
-o FILE 设置最终排序后的输出文件名；
-O FORMAT 设置最终输出的文件格式，可以是bam，sam或者cram，默认为bam；
-T PREFIX 设置临时文件的前缀；
-@ INT 设置排序和压缩是的线程数量，默认是单线程。
例子：
当有多个样本的bam文件时，可以使用samtools的merge命令将这些bam文件进行合并为一个bam文件。。Merge命令将多个已经排序后的bam文件合并成为一个排序的且保持所有输入记录并保持现有排序顺序的bam文件。-bam文件转换成是sam文件格式
Options:
-1 指定压缩等级；
-b FILE 输入文件列表，一个文件一行；
-f overwrite the output BAM if exist 强制覆盖同名输出文件；
-h FILE copy the header in FILE to out.bam [in1.bam]。 指定FILE内的’@’头复制到输出bam文件中并替换输出文件的文件头。否则，输出文件的文件头将从第一个输入文件复制过来；-bam文件转换成是sam文件格式
-n sort by read names。设定输入比对文件是以read名进行排序的而不是以染色体坐标排序的；
-R STR merge file in the specified region STR [all]。合并输入文件的指定区域；
-r attach RG tag (inferred from file names)。使RG标签添加到每一个比对文件上，标签值来自文件名；
-u uncompressed BAM output。输出的bam文件不压缩；
-c 当多个输入文件包含相同的@RG头ID时，只保留第一个到合并后输出的文件。当合并多个相同样本的不同文件时，非常有用。
-p 与-c参数类似，对于要合并的每一个文件中的@PG ID只保留第一个文件中的@PG。
例子
Note: Samtools' merge does not reconstruct the @RG dictionary in the header. Users must provide the correct header with -h, or uses Picard which properly maintains the header dictionary in merging.-bam文件转换成是sam文件格式
为了能够快速访问bam文件，可以为已经基于坐标排序后bam或者cram的文件创建索引，生成以.bai或者.crai为后缀的索引文件。必须使用排序后的文件，否则可能会报错。另外，不能对sam文件使用此命令。如果想对sam文件建立索引，那么可以使用tabix命令创建。-bam文件转换成是sam文件格式
Options:
-b 创建bai索引文件，未指定输出格式时，此参数为默认参数；
-c 创建csi索引文件，默认情况下，索引的最小间隔值为2^14，与bai格式一致；
-m INT 创建csi索引文件，最小间隔值2^INT；
例子：
以下两种命令结果一样
对fasta文件建立索引，生成的索引文件以.fai后缀结尾。该命令也能依据索引文件快速提取fasta文件中的某一条（子）序列．如果没有指定区域，faidx命令就创建文件索引并生成后缀为.fai的索引文件。如果指定区域，那么就是生产并显示fasta格式的子序列。输入文件可以使BGZF压缩格式的文件。-bam文件转换成是sam文件格式
输入文件中的序列要有不同的名称。如果不是这样，即存在相同名称的序列，在建立索引的过程中将发出有关重复序列的警告而且生产的同名子序列的信息都要被第一个同名子序列的信息覆盖。
例：对基因组文件建立索引
上图中先显示了待处理的test.fasta文件的内容，由4个长度不一的序列组成，其中前两个重名。然后使用faidx命令进行处理。最后在显示生产索引文件test.fasta.fai的内容。
生成了索引文件fasta.fai,是一个文本文件，分成了5列。
于是通过此文件，可以定位子序列在fasta文件在磁盘上的存放位置，直接快速调出子序列。
由于有索引文件，可以使用以下命令很快从基因组中提取到fasta格式的子序列
提取序列如下图：
tview能直观的显示出reads比对基因组的情况，和基因组浏览器有点类似。
当给出参考基因组的时候，会在第一排显示参考基因组的序列，否则，第一排全用N表示。
按下 g ，则提示输入要到达基因组的某一个位点。例子“scaffold_10:1000"表示到达第10号scaffold的第1000个碱基位点处。
使用H(左）J（上）K（下）L（右）移动显示界面。大写字母移动快，小写字母移动慢。
使用空格建向左快速移动（和 L 类似），使用Backspace键向左快速移动（和 H 类似）。
Ctrl+H 向左移动1kb碱基距离； Ctrl+L 向右移动1kb碱基距离
可以用颜色标注比对质量，碱基质量，核苷酸等。30～40的碱基质量或比对质量使用白色表示；20～30黄色；10～20绿色；0～10蓝色。
使用点号'.'切换显示碱基和点号；使用r切换显示read name等
还有很多其它的使用说明，具体按 ？ 键来查看。
统计输入文件的相关数据并将这些数据输出至屏幕显示。每一项统计数据都由两部分组成，分别是QC pass和QC failed，表示通过QC的reads数据量和未通过QC的reads数量。以“PASS + FAILED”格式显示。-bam文件转换成是sam文件格式
例如
得到每个碱基位点或者区域的测序深度,并输出到标准输出。depth命令计算每一个位点的测序深度并在标准显示设备中显示。 注意 ：使用此命令之前必须先samtools index。
Options:
-a 输出所有位点，包括零深度的位点；-a –a，--aa 完全输出所有位点，包括未使用到的参考序列；
-b FILE 计算BED文件中指定位置或区域的深度；
-f FiLE 使用在FILE中的指定bam文件；
-I INT 忽略掉小于此INT值的reads；
-q INT 只计算碱基质量值大于此值的reads；
-Q INT 只计算比对质量值大于此值的reads；
-r CHR:FROM –TO 只计算指定区域的reads，后面跟染色体号（region）
示例：
注意： in.bam 必须经过了排序。
depth命令运行如下图所示：
一共得到3列以指标分隔符分隔的数据，第一列为染色体名称，第二列为位点，第三列为覆盖深度。
该命令用于生成bcf文件，再使用bcftools进行SNP和Indel的分析。bcftools是samtool中附带的软件，在samtools的安装文件夹中可以找到。
Options:
-f –来输入有索引文件的fasta参考序列；
-F –gap-frac FLOAT 含有间隔reads的最小片段。
-l –positions FILE BED文件或者包含区域位点的位置列表文件。位置文件包含两列，染色体和位置，从1开始计数。BED文件至少包含3列，染色体、开始位置和结束位置，开始端从0开始计数。-bam文件转换成是sam文件格式
-r –region STR 只在指定区域产生pileup，需要已建立索引的bam文件。通常和-l参数一起使用。
-o –output FILE 生成pileup格式文件或者VCF、BCF文件而不是默认的标准输出。
-g –BCF 计算基因型的似然值和输出文件格式为BCF。
-v –VCF 计算基因型的似然值和输出文件格式为VCF。
-C --adjust-MQ INT 用于降低比对质量的系数，如果reads中含有过多的错配。不能设置为零。BWA推荐值为50。
-A --count-orphans 在检测变异中，不忽略异常的reads对。
-D –输出每个样本的reads深度。
-V –输出每个样本未比对到参考基因组的reads数量。
-t –output-tags LIST设置FORMAT和INFO的列表内容，以逗号分割。
-u –uncompressed 生成未压缩的VCF和BCF文件。
-I –skip-indel 不检测INDEL。
-m –min-ireads INT 候选INDEL的最小间隔的reads。
下面是一个使用-r参数和-l参数生成vcf文件的实例：
mpileup不使用-u或-g参数时，则不生成二进制的bcf文件，而生成一个文本文件(输出到标准输出)。该文本文件统计了参考序列中每个碱基位点的比对情况；该文件每一行代表了参考序列中某一个碱基位点的比对结果。比如：-bam文件转换成是sam文件格式
mpileup生成的结果包含6行：参考序列名；位置；参考碱基；比对上的reads数；比对情况；比对上的碱基的质量。
其中第5列比较复杂,解释如下：
NGS上机测序前需要进行PCR一步，使一个模板扩增出一簇，从而在上机测序的时候表现出为1个点，即一个reads。若一个模板扩增出了多簇，结果得到了多个reads，这些reads的坐标(coordinates)是相近的。在进行了reads比对后需要将这些由PCR duplicates获得的reads去掉，并只保留最高比对质量的read。使用rmdup命令即可完成.-bam文件转换成是sam文件格式
Options:
-s 对single-end reads。默认情况下，只对paired-end reads
-S 将Paired-end reads作为single-end reads处理。
检索和打印已建立索引的bam文件的统计数据，包括参考序列名称、序列长度、比对上的read数量和未比对上的read数量。输出结果显示在屏幕上，以制表符分割。
如下图所示:
idxstats 统计一个表格，4列，分别为”序列名，序列长度，比对上的reads数，unmapped reads number”。第4列应该是paired reads中有一端能匹配到该scaffold上，而另外一端不匹配到任何scaffolds上的reads数。-bam文件转换成是sam文件格式
计算由BED文件指定的基因组区域内的总碱基数量。
如下图所示：
reheader 替换bam文件的头
cat 连接多个bam文件，适用于非sorted的bam文件
有时候，我们需要提取出比对到一段参考序列的reads，进行小范围的分析，以利于debug等。这时需要将bam或sam文件转换为fastq格式。
该网站提供了一个bam转换为fastq的程序：
1：从bam中过滤掉没有比对上的信息，并将比对上的部分保存到sam中
2：从bam中过滤没有比对上的信息，并保存到bam中：(要加-h,表示输出header）
注意：虽然有些没有比对上的结果 的flag值是141,77等，但是都可以用4过滤掉。
3：提取为fastq：主要参数为-f,-F,-G三个。
提取没有比对上的信息到fastq用：
提取比对上的reads到fastq用
REF:
利用samtools将sam格式的文件与bam格式的文件进行相互转换
bowtie2是当今流行的序列比对软件，其输出结果为sam后缀名的文件
sam格式是一种通用的比对格式，用来存储reads到参考序列的比对信息SAM是一种序列比对格式标准， 由sanger制定，是以TAB为分割符的文本格式。
主要应用于测序序列mapping到基因组上的结果表示，当然也可以表示任意的多
重比对结果
而bam格式文件可以理解为时sam格式文件的二进制保存
在进行下一步的转录本组装时要用到cufflinks软件，而cufflinks只接受bam格式的文件作为输入，所以我们要把sam格式的文件转换为bam格式的文件以便进行下一步操作 samtools可以有效地帮我们解决这个问题-bam文件转换成是sam文件格式
samtools view [-bhuHS] [-t in.reList] [-o output] [-f repFlag] [-F skipFlag] [-q minMapQ] [-l library]-bam文件转换成是sam文件格式
[-r read]
-b 以BAM格式输出，可以用于samtools的后续分析
-u 以未压缩的BAM格式输出，可以节约时间，一般在管道执行时使用
-h 在结果中包含头header
-H 只输出头 -S 输入文件为SAM格式，如果确实@SQ头，则需要-t选项
sam转化为bam
samtools view -bS aln.sam aln.bam
bam转化为sam
samtools view -h -o aln.sam aln.bam
另外在利用cufflinks对转录本进行拼接时，cufflinks还需要我们把转换后的bam格式文件进行排序
samtools sort aln.bam aln.sorted_bam
建议使用tophat2+cufflinks的软件组合进行转录组的比对和分析
具体教程会在后面更新
fasta、sam、bam相互转换
1. fasta = sam
2. fasta = sam
1. sam = bam
2. sam = bam
1. fasta = bam
2. fasta = bam
我用到的Samtools介绍
记录一下我用到的samtools的用法。
samtools的说明文档：
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作为过滤参数即可-bam文件转换成是sam文件格式
$ 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区域的比对结果-bam文件转换成是sam文件格式
$ samtools view -T genome.fasta -h scaffold1.sam scaffold1.h.sam #根据fasta文件，将 header 加入到 sam 或 bam 文件中-bam文件转换成是sam文件格式
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。-bam文件转换成是sam文件格式
也就是说：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]-bam文件转换成是sam文件格式
例如：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 会被放在文件的末尾。 -bam文件转换成是sam文件格式
参数：
-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的比对图形的时候也需要。-bam文件转换成是sam文件格式
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命令即可完成.-bam文件转换成是sam文件格式
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的安装文件夹中可以找到。 -bam文件转换成是sam文件格式
最常用的参数有2个：
-f来输入有索引文件的fasta参考序列；
-g输出到bcf格式。用法和最简单的例子如下
Usage:samtoolsmpileup[-EBug][-CcapQcoef][-rreg][-fin.fa][-llist][-McapMapQ][-QminBaseQ][-qminMapQ]in.bam[in2.bam[...]]-bam文件转换成是sam文件格式
$samtoolsmpileup-fgenome.fastaabc.bamabc.txt
$samtoolsmpileup-gSDfgenome.fastaabc.bamabc.bcf
$samtoolsmpileup-guSDfgenome.fastaabc.bam|\bcftoolsview-cvNg-abc.vcf
mpileup不使用-u或-g参数时，则不生成二进制的bcf文件，而生成一个文本文件(输出到标准输出)。该文本文件统计了参考序列中每个碱基位点的比对情况；该文件每一行代表了参考序列中某一个碱基位点的比对结果。比如：-bam文件转换成是sam文件格式
（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。该命令的使用方法和例子为：-bam文件转换成是sam文件格式
$ 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)。-bam文件转换成是sam文件格式
参考链接：
原文链接：
samtools使用大全
samtools是一个用于操作sam和bam文件（通常是短序列比对工具如bwa,bowtie2,hisat2,tophat2等等产生的，具体格式可以在消息框输入“SAM”查看）的工具合集，包含有许多命令。以下是常用命令的介绍。-bam文件转换成是sam文件格式
1.View
view命令的主要功能是：将sam文件与bam文件互换；然后对bam文件进行各种操作，比如数据的排序(sort)和提取(这些操作 是对bam文件进行的，因而当输入为sam文件的时候，不能进行该操作)；最后将排序或提取得到的数据输出为bam或sam（默认的）格式。-bam文件转换成是sam文件格式
bam文件优点：bam文件为二进制文件，占用的磁盘空间比sam文本文件小；利用bam二进制文件的运算速度快。
view命令中，对sam文件头部（序列ID）的输入(-t或-T）和输出(-h)是单独的一些参数来控制的。
Usage: samtools view [options] in.bam|in.sam [region1 [...]]
默认情况下不加 region，则是输出所有的 region.
options:
-b output BAM
默认下输出是 SAM 格式文件，该参数设置输出 BAM 格式
-h print header for the SAM output
默认下输出的 sam 格式文件不带 header，该参数设定输出sam文件时带 header 信息
-H print header only (no alignments)
仅仅输出文件的头
-S input is SAM
默认下输入是 BAM 文件，若是输入是 SAM 文件，则最好加该参数，否则有时候会报错。
-u uncompressed BAM output (force -b)
该参数的使用需要有-b参数，能节约时间，但是需要更多磁盘空间。
-c Instead of printing the alignments, only count them and print the
total number. All filter options, such as ‘-f’, ‘-F’ and ‘-q’ , are taken into account.
过滤功统计功能
-c print only the count of matching records
-L FILE output alignments overlapping the input BED FILE [null]
-t FILE list of reference names and lengths (force -S) [null]
使用一个list文件来作为header的输入
-T FILE reference sequence file (force -S) [null]
使用序列fasta文件作为header的输入
-o FILE output file name [stdout]
-F INT filtering flag, 0 for unset [0]
Skip alignments with bits present in INT [0]
数字4代表该序列没有比对到参考序列上
数字8代表该序列的mate序列没有比对到参考序列上
过滤功能。如F12过滤只有双端map的
-q INT minimum mapping quality [0]
比对的最低质量值，一般认为20就为unique比对了，可以结合上述-bF参数使用使用提取特定的比对结果
例子：
将sam文件转换成bam文件
samtools view -bS abc.sam abc.bam
BAM转换为SAM
samtools view -h -o out.sam out.bam
提取比对到参考序列上的比对结果
samtools view -bF 4 abc.bam abc.F.bam
提取paired reads中两条reads都比对到参考序列上的比对结果，只需要把两个4+8的值12作为过滤参数即可
samtools view -bF 12 abc.bam abc.F12.bam
提取没有比对到参考序列上的比对结果
samtools view -bf 4 abc.bam abc.f.bam
提取bam文件中比对到caffold1上的比对结果，并保存到sam文件格式
samtools view abc.bam scaffold1 scaffold1.sam
提取scaffold1上能比对到30k到100k区域的比对结果
samtools view abc.bam scaffold1:30000-100000 $gt; scaffold1_30k-100k.sam
根据fasta文件，将 header 加入到 sam 或 bam 文件中
samtools view -T genome.fasta -h scaffold1.sam scaffold1.h.sam
2.Sort
sort对bam文件进行排序。一些软件需要sort的bam或者sam文件，如stringtie，所以必须要sort使用；求depth时，也必须要sort；
Usage: samtools sort [-n] [-m maxMem] in.bam out.prefix
-m 内存参数默认下是 500,000,000 即500M（不支持K，M，G等缩写）。对于处理大数据时，如果内存够用，则设置大点的值，以节约时间。
-n 设定排序方式按short reads的ID排序。默认下是按序列在fasta文件中的顺序（即header）和序列从左往右的位点排序。
例子：
samtools sort accept.bam accept.sort最终产生accept.sort.bam
3.merge
将2个或2个以上的已经sort了的bam文件融合成一个bam文件。融合后的文件已经sort过了的。
Usage: samtools merge [-nr] [-h inh.sam] out.bam in1.bam in2.bam[...]
Options: -n sort by read names
-r attach RG tag (inferred from file names)
-u uncompressed BAM output
-f overwrite the output BAM if exist
-1 compress level 1
-R STR merge file in the specified region STR [all]
-h FILE copy the header in FILE to out.bam [in1.bam]
例子：
4.index
必须对bam文件进行默认情况下的排序后，才能进行index。否则会报错。
建立索引后将产生后缀为.bai的文件，用于快速的随机处理。很多情况下需要有bai文件的存在，特别是显示序列比对情况下。比如samtool的tview命令就需要；gbrowse2显示reads的比对图形的时候也需要。IGV显示比对情况也需要。-bam文件转换成是sam文件格式
Usage: samtools index in.bam [out.index]
例子：
以下两种命令结果一样
$ samtools index abc.sort.bam
$ samtools index abc.sort.bam abc.sort.bam.bai
5.faidx
对fasta文件建立索引,生成的索引文件以.fai后缀结尾。该命令也能依据索引文件快速提取fasta文件中的某一条（子）序列
Usage: samtools faidx in.bam [ [...]]
对基因组文件建立索引，方便提取序列。
例子：$ samtools faidx genome.fasta
由于有索引文件，可以使用以下命令很快从基因组中提取到fasta格式的子序列
$ samtools faidx genome.fasta scffold_10 scaffold_10.fasta
6.tview
tview能直观的显示出reads比对基因组的情况，和基因组浏览器有点类似。
需要事先利用利用上面讲的sort和建index命令执行完后，用下述命令。
Usage: samtools tview aln.bam [ref.fasta]
出参考基因组的时候，会在第一排显示参考基因组的序列，否则，第一排全用N表示。
按下 g ，则提示输入要到达基因组的某一个位点。例子“scaffold_10:1000"表示到达第
10号scaffold的第1000个碱基位点处。
使用H(左）J（上）K（下）L（右）移动显示界面。大写字母移动快，小写字母移动慢。
使用空格建向左快速移动（和 L 类似），使用Backspace键向左快速移动（和 H 类似）。
Ctrl+H 向左移动1kb碱基距离； Ctrl+L 向右移动1kb碱基距离
可以用颜色标注比对质量，碱基质量，核苷酸等。30～40的碱基质量或比对质量使用白色表示；
20～30黄色；10～20绿色；0～10蓝色。
使用点号'.'切换显示碱基和点号；使用r切换显示read name等
还有很多其它的使用说明，具体按 ？ 键来查看。
7.flagstat
给出BAM文件的比对结果
Usage: samtools flagstat in.bam
$ samtools flagstat example.bam
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)
#同上一个，只是其中比对质量=5的reads的数量
8.depth
得到每个碱基位点的测序深度,并输出到标准输出，所以要用大于号追加到一个文件。
Usage: bam2depth [-r reg] [-q baseQthres] [-Q mapQthres] [-b in.bed] in1.bam [...]
-r 后面跟染色体号（region）
-q ：计算深度时要求测序碱基质量最低质量值
-Q ：计算深度时要求比对的最低质量值
注意：做depth之前必须做samtools index；
例子
samtools depth accept.bam depth
9.其他命令
reheader：替换bam文件的头
$ samtools reheader in.header.sam in.bam
idxstats ：统计一个表格，4列，分别为”序列名，序列长度，比对上的reads数，unmapped reads number。第4列应该是paired reads中有一端能匹配到该scaffold上，而另外一端不匹配到任何scaffolds上的reads数。-bam文件转换成是sam文件格式
$ samtools idxstats aln.bam
rmdup：将由PCR duplicates获得的reads去掉，并只保留最高比对质量的read。
Usage: samtools rmdup [-sS]
-s 对single-end reads。默认情况下，只对paired-end reads
-S 将Paired-end reads作为single-end reads处理。
10. 将bam文件转换为fastq文件
有时候，我们需要提取出比对到一段参考序列的reads，进行小范围的分析，以利于debug等。这时需要将bam或sam文件转换为fastq格式。
该网站提供了一个bam转换为fastq的程序：
$ wget
$ tar zxf bam2fastq-1.1.0.tgz
$ cd bam2fastq-1.1.0
$ make
$ ./bam2fastq in.bam
11. mpileup
samtools还有个非常重要的命令mpileup，以前为pileup。该命令用于生成bcf文件，再使用bcftools进行SNP和Indel的分析。bcftools是samtool中附带的软件，在samtools的安装文件夹中可以找到。-bam文件转换成是sam文件格式
最常用的参数有2：
-f 来输入有索引文件的fasta参考序列； -g 输出到bcf格式。用法和最简单的例子如下
Usage: samtools mpileup [-EBug] [-C capQcoef] [-r reg] [-f in.fa] [-l list] [-M capMapQ] [-Q minBaseQ] [-q minMapQ] in.bam [in2.bam [...]]$ samtools mpileup -f genome.fasta abc.bam abc.txt-bam文件转换成是sam文件格式
$ samtools mpileup -gSDf genome.fasta abc.bam abc.bcf
$ samtools mpileup -guSDf genome.fasta abc.bam | \
bcftools view -cvNg - abc.vcf
mpileup不使用-u或-g参数时，则不生成二进制的bcf文件，而生成一个文本文件(输出到标准输出)。该文本文件统计了参考序列中每个碱基位点的比对情况；该文件每一行代表了参考序列中某一个碱基位点的比对结果。比如：-bam文件转换成是sam文件格式
scaffold_1 2841 A 11 ,,,...,.... BHIGDGIJ?FF
scaffold_1 2842 C 12 ,$,,...,....^I. CFGEGEGGCFF+
scaffold_1 2843 G 11 ,,...,..... FDDDDCD?DD+
scaffold_1 2844 G 11 ,,...,..... FA?AAAAAA+
scaffold_1 2845 G 11 ,,...,..... F656666166*
scaffold_1 2846 A 11 ,,...,..... (1.1111)11*
scaffold_1 2847 A 11 ,,+9acggtgaag.+9ACGGTGAAT.+9ACGGTGAAG.+9ACGGTGAAG,+9acggtgaag.+9ACGGTGAAG.+9ACGGTGAAG.+9ACGGTGAAG.+9ACGGTGAAG.+9ACGGTGAAG %.+....-..)-bam文件转换成是sam文件格式
scaffold_1 2848 N 11 agGGGgGGGGG !!$!!!!!!!!
scaffold_1 2849 A 11 c$,...,..... !0000000000
scaffold_1 2850 A 10 ,...,..... 353333333
mpileup生成的结果包含6行：参考序列名；位置；参考碱基；比对上的reads数；比对情况；比对上的碱基的质量。其中第5列比较复杂,解释如下：
1 ‘.’代表与参考序列正链匹配。
2 ‘,’代表与参考序列负链匹配。
3 ‘ATCGN’代表在正链上的不匹配。
4 ‘atcgn’代表在负链上的不匹配。
5 ‘*’代表模糊碱基
6 ‘^’代表匹配的碱基是一个read的开始；’^'后面紧跟的ascii码减去33代表比对质量；这两个符号修饰的是后面的碱基，其后紧跟的碱基(.,ATCGatcgNn)代表该read的第一个碱基。
7 ‘$’代表一个read的结束，该符号修饰的是其前面的碱基。
8 正则式’\+[0-9]+[ACGTNacgtn]+’代表在该位点后插入的碱基；比如上例中在scaffold_1的2847后插入了9个长度的碱基acggtgaag。表明此处极可能是indel。
9 正则式’-[0-9]+[ACGTNacgtn]+’代表在该位点后缺失的碱基；
12. 使用bcftools
bcftools和samtools类似，用于处理vcf(variant call format)文件和bcf(binary call format)文件。前者为文本文件，后者为其二进制文件。
bcftools使用简单，最主要的命令是view命令，其次还有index和cat等命令。index和cat命令和samtools中类似。此处主讲使用view命令来进行SNP和Indel calling。该命令的使用方法和例子为：-bam文件转换成是sam文件格式
$ bcftools view [-AbFGNQSucgv] [-D seqDict] [-l listLoci] [-s listSample]
[-i gapSNPratio] [-t mutRate] [-p varThres] [-P prior]
[-1 nGroup1] [-d minFrac] [-U nPerm] [-X permThres]
[-T trioType] in.bcf [region]
$ 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)。-bam文件转换成是sam文件格式
例如：
scaffold_1 2847 . A AACGGTGAAG 194 . INDEL;DP=11;VDB=0.0401;AF1=1;AC1=2;DP4=0,0,8,3;MQ=35;FQ=-67.5 GT:PL:GQ 1/1:235,33,0:63-bam文件转换成是sam文件格式
scaffold_1 3908 . G A 111 . DP=13;VDB=0.0085;AF1=1;AC1=2;DP4=0,0,5,7;MQ=42;FQ=-63 GT:PL:GQ 1/1:144,36,0:69-bam文件转换成是sam文件格式
scaffold_1 4500 . A G 31.5 . DP=8;VDB=0.0034;AF1=1;AC1=2;DP4=0,0,1,3;MQ=42;FQ=-39 GT:PL:GQ 1/1:64,12,0:21-bam文件转换成是sam文件格式
scaffold_1 4581 . TGGNGG TGG 145 . INDEL;DP=8;VDB=0.0308;AF1=1;AC1=2;DP4=0,0,0,8;MQ=42;FQ=-58.5 GT:PL:GQ 1/1:186,24,0:45-bam文件转换成是sam文件格式
scaffold_1 4644 . G A 195 . DP=21;VDB=0.0198;AF1=1;AC1=2;DP4=0,0,10,10;MQ=42;FQ=-87 GT:PL:GQ 1/1:228,60,0:99-bam文件转换成是sam文件格式
scaffold_1 4827 . NACAAAGA NA 4.42 . INDEL;DP=1;AF1=1;AC1=2;DP4=0,0,1,0;MQ=40;FQ=-37.5 GT:PL:GQ 0/1:40,3,0:3-bam文件转换成是sam文件格式
scaffold_1 4854 . A G 48 . DP=6;VDB=0.0085;AF1=1;AC1=2;DP4=0,0,2,1;MQ=41;FQ=-36 GT:PL:GQ 1/1:80,9,0:16-bam文件转换成是sam文件格式
scaffold_1 5120 . A G 85 . DP=8;VDB=0.0355;AF1=1;AC1=2;DP4=0,0,5,3;MQ=42;FQ=-51 GT:PL:GQ 1/1:118,24,0:45-bam文件转换成是sam文件格式
第8列中显示了对variants的信息描述，比较重要，其中的 Tag 的描述如下：
Tag Format Description
AF1 double Max-likelihood estimate of the site allele frequency (AF) of the first ALT allele
DP int Raw read depth (without quality filtering)
DP4 int[4] # high-quality reference forward bases, ref reverse, alternate for and alt rev bases
FQ int Consensus quality. Positive: sample genotypes different; negative: otherwise
MQ int Root-Mean-Square mapping quality of covering reads
PC2 int[2] Phred probability of AF in group1 samples being larger (,smaller) than in group2
PCHI2 double Posterior weighted chi^2 P-value between group1 and group2 samples
PV4 double[4] P-value for strand bias, baseQ bias, mapQ bias and tail distance bias
QCHI2 int Phred-scaled PCHI2
RP int # permutations yielding a smaller PCHI2
CLR int Phred log ratio of genotype likelihoods with and without the trio/pair constraint
UGT string Most probable genotype configuration without the trio constraint
CGT string Most probable configuration with the trio constraint
使用bcftools得到variant calling结果后。需要对结果再次进行过滤。主要依据比对结果中第8列信息。其中的 DP4 一行尤为重要，提供了4个数据：1 比对结果和正链一致的reads数、2 比对结果和负链一致的reads数、3 比对结果在正链的variant上的reads数、4 比对结果在负链的variant上的reads数。可以设定 （value3 + value4）大于某一阈值，才算是variant。比如：-bam文件转换成是sam文件格式
$ perl -ne 'print $_ if /DP4=(\d+),(\d+),(\d+),(\d+)/ ($3+$4)=10 ($3+$4)/($1+$2+$3+$4)=0.8' snp_indel.vcf snp_indel.final.vcf-bam文件转换成是sam文件格式
