STEP5——比对数据处理:Samtools(2)
2016-09-15 12:27阅读:
原文转载自:http://www.cnblogs.com/emanlee/p/4316581.html
6. tview
tview能直观的显示出reads比对基因组的情况,和基因组浏览器有点类似。
Usage: samtools tview [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
$ 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] [...]
-r
后面跟染色体号(region)
注意:做depth之前必须做samtools
index;
9.reheader
替换bam文件的头
$ samtools reheader
cat
连接多个bam文件,适用于非sorted的bam文件
$ samtools cat [-h header.sam] [-o out.bam] [ ...
]
idxstats
统计一个表格,4列,分别为”序列名,序列长度,比对上的reads数,unmapped
reads
number”。第4列应该是paired
reads中有一端能匹配到该scaffold上,而另外一端不匹配到任何scaffolds上的reads数。
$ samtools idxstats
10.
将bam文件转换为fastq文件
有时候,我们需要提取出比对到一段参考序列的reads,进行小范围的分析,以利于debug等。这时需要将bam或sam文件转换为fastq格式。
该网站提供了一个bam转换为fastq的程序:
http://www.hudsonalpha.org/gsl/information/software/bam2fastq
$ wget
http://www.hudsonalpha.org/gsl/static/software/bam2fastq-1.1.0.tgz
$ tar zxf bam2fastq-1.1.0.tgz
$ cd bam2fastq-1.1.0
$ make
$ ./bam2fastq
11. mpileup
samtools还有个非常重要的命令mpileup,以前为pileup。该命令用于生成bcf文件,再使用bcftools进行SNP和Indel的分析。bcftools是samtool中附带的软件,在samtools的安装文件夹中可以找到。
最常用的参数有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
$ 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文件,而生成一个文本文件(输出到标准输出)。该文本文件统计了参考序列中每个碱基位点的比对情况;该文件每一行代表了参考序列中某一个碱基位点的比对结果。比如:
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?AAAA
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
%.+....-..)
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个长度的碱基acg