新浪博客

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'表示到达第
10scaffold的第1000个碱基位点处。
使用H(左)J(上)K(下)L(右)移动显示界面。大写字母移动快,小写字母移动慢。
使用空格建向左快速移动(和 L 类似),使用Backspace键向左快速移动(和 H 类似)。
Ctrl+H 向左移动
1kb碱基距离; Ctrl+L 向右移动1kb碱基距离
可以用颜色标注比对质量,碱基质量,核苷酸等。3040的碱基质量或比对质量使用白色表示;
2030黄色;1020绿色;010蓝色。
使用点号'.'切换显示碱基和点号;使用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)
#同上一个,只是其中比对质量>=5reads的数量
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文件,适用于非sortedbam文件
$ 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等。这时需要将bamsam文件转换为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进行SNPIndel的分析。bcftoolssamtool中附带的软件,在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_12847后插入了9个长度的碱基acg

我的更多文章

下载客户端阅读体验更佳

APP专享