samtools 学习记录

SAMtools不仅仅用来call snp。从samtools的软件名就能看出，是对SAM格式文件进行操作的工作，比如讲sam转成bam格式，index，rmdup等等。samtools结合linux命令比如grep，awk和SAM格式描述的flag，tag，亦是非常非常非常强大，比如根据flag过滤duplicate的reads，根据XA tag过滤multiple hit的reads。本文在此只介绍一下samtools的统计命令，能快速对bam文件进行各种统计。

samtools flagstat

给出BAM文件的比对结果

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

全部reads数

0 + 0 secondary
146415 + 0 supplementary

5889334 + 0 duplicates

重复

97455019 + 0 mapped (99.92% : N/A)

总体上reads的匹配率

97391220 + 0 paired in sequencing

有多少reads是属于paired reads

48695610 + 0 read1

reads1中的reads数

48695610 + 0 read2

reads2中的reads数

95449424 + 0 properly paired (98.01% : N/A)

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

97249016 + 0 with itself and mate mapped

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

59588 + 0 singletons (0.06% : N/A)

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

1552228 + 0 with mate mapped to a different chr

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

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

同上一个，只是其中比对质量>=5的reads的数量

samtools stats <bam file>

输出的信息比较多，部分如下Summary Numbersraw total sequences，filtered sequences, reads mapped, reads mapped and paired,reads properly paired等信息

输出的信息比较多，部分如下Summary Numbersraw total sequences，filtered sequences, reads mapped, reads mapped and paired,reads properly paired等信息

SN raw total sequences:97391220
SN filtered sequences:0
SN sequences: 97391220
SN is sorted: 0
SN 1st fragments:48695610
SN last fragments:48695610
SN reads mapped:97308604
SN reads mapped and paired:97249016# paired-end technology bit set + both mates mapped
SN reads unmapped:82616
SN reads properly paired:95449424# proper-pair bit set
SN reads paired:97391220# paired-end technology bit set
SN reads duplicated:0# PCR or optical duplicate bit set
SN reads MQ0: 2180859 # mapped and MQ=0
SN reads QC failed:0
SN non-primary alignments:0
SN total length:11989013517# ignores clipping
SN bases mapped:11981756912# ignores clipping
SN bases mapped (cigar):11959098931# more accurate
SN bases trimmed:0
SN bases duplicated:0
SN mismatches:41386410# from NM fields
SN error rate:3.460663e-03# mismatches / bases mapped (cigar)
SN average length:123
SN maximum length:125
SN average quality:35.2
SN insert size average:213.4
SN insert size standard deviation:61.5
SN inward oriented pairs:45465946
SN outward oriented pairs:40738
SN pairs with other orientation:26409
SN pairs on different chromosomes:776114

samtools  tview  <.bam>  <ref.fasta>

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等
还有很多其它的使用说明，具体按 ？ 键来查看。

samtools faidx

 从参考序列获取一段序列

http://www.htslib.org/doc/samtools.html