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 secondary146415 + 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