samtools操作步骤

SAM由头文件和map结果组成。头文件由一行行以@起始的注释构成。而map结果是类似下面的东西：

HWI-ST1001:137:C12FPACXX:7:1115:14131:66670 0 chr1 12805 1 42M4I5M * 0 0 TTGGATGCCCCTCCACACCCTCTTGATCTTCCCTGTGATGTCACCAATATG CCCFFFFFHHGHHJJJJJHJJJJJJJJJJJJJJJJIJJJJJJJJJJJJIJJ AS:i:-28 XN:i:0 XM:i:2 XO:i:1XG:i:4 NM:i:6 MD:Z:2C41C2 YT:Z:UU NH:i:3 CC:Z:chr15 CP:i:102518319 XS:A:+ HI:i:0   HWI-ST1001:137:C12FPACXX:7:2313:17391:30032 272 chr1 13494 1 51M * 0 0 ACTGCCTGGCGCTGTGCCCTTCCTTTGCTCTGCCCGCTGGAGACAGTGTTT CFFFFHHJJJJIJJJJIJJJJJJJJJJJJJJJJJJJJJHHHHFA+FFFC@B AS:i:-3 XN:i:0 XM:i:1 XO:i:0 XG:i:0NM:i:1 MD:Z:44G6 YT:Z:UU XS:A:+ NH:i:3 CC:Z:chr15 CP:i:102517626 HI:i:0   HWI-ST1001:137:C12FPACXX:7:1109:17518:53305 16 chr1 13528 1 51M * 0 0 CGCTGGAGCCGGTGTTTGTCATGGGCCTGGGCTGCAGGGATCCTGCTACAA #############AB=?:*B?;A?<2+233++;A+A2+<7==@7,A<A<=> AS:i:-5 XN:i:0 XM:i:2 XO:i:0 XG:i:0NM:i:2 MD:Z:8A21T20 YT:Z:UU XS:A:+ NH:i:4 CC:Z:chr15 CP:i:102517592 HI:i:0

看上去很类似fastq文件，它也有read名称，序列，质量等信息，但是又不完全一样。首先，每个read只占一行，只是它被tab分成了很多列，一共有12列，分别记录了：

1. read名称
2. SAM标记
3. chromosome
4. 5′端起始位置
5. MAPQ（mapping quality，描述比对的质量，数字越大，特异性越高）
6. CIGAR字串，记录插入，删除，错配以及splice junctions（后剪切拼接的接头）
7. mate名称，记录mate pair信息
8. mate的位置
9. 模板的长度
10. read序列
11. read质量
12. 程序用标记

显然，其中chromosome至CIGAR的信息都是非常重要的。但是这些对我们不重要，我们只需要了解SAM/BAM文件是什么，就可以了。重要的是如果进行下游的操作。

要操作SAM/BAM文件，首先需要安装samtools。它的安装过程和所有的linux/unix程序一样，都是经过make之后生成可执行程序，然后把它的路径告知系统，或者放在系统可以找到的位置就可以了。
比如：

tar zxvf samtools-0.1.18.tar.bz2 cd samtools-0.1.18/ make samtoolpath=`pwd` PATH=PATH:$samtoolpath

然后就可以按照samtools主页上介绍的工具进行各种操作了。我们最常见的几步操作比如
0. SAM，BAM转换

samtools view -h file.bam > file.sam samtools view -b -S file.sam > file.bam

1. sorting BAM文件。大多数下游程序都要求BAM文件是被排过序的。

samtools sort file.bam outputPrefix

2. 创建BAM index。这也是被大多数下游程序所要求。

samtools index sorted.bam

3. index模板基因组。这也是被大多数下游程序所要求。

samtools faidx Homo_sapiens_assembly19.fasta

在很多时候，我们还会看到一种扩展名为BED的mapping文件。其具体格式也是几经变化，但是现在以UCSC的描述为准。从BAM文件转换成BED文件，我们需要安装BEDtools。下载安装就不多说了。示例一个如何从BAM文件转换成BED文件的命令：

bamToBed -i reads.bam > reads.bed

samtools 可以应用于linux的命令管道里。以“-”作为标准输入或输出。

view   从bam/sam文件中提取/打印部分比对结果。默认为所有的区域，也可以染色体区域（1-based，须sort并index）。

例如：

samtools view -bt ref_list.txt -o aln.bam aln.sam.gz

samtools view aln.sorted.bam chr2:20,100,000-20,200,000

tview 

mpileup 列举每条reads比对的indel，SNP等信息。

“.”表示match;

“，”表示反向match;

“<”或“>”表示reference skip;

"ATCGN"表示正向mismatch;

"atcgn"表示反向mismatch;

‘+[0-9]+[ACGTNacgtn]+’ insertion;

‘-[0-9]+[ACGTNacgtn]+’ 表示deletion;

“^”标记reads起始;

“$”标记reads segment结尾;

samtools pileup -vcf ref.fasta aln.sorted.bam

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

samtools mpileup -C50 -gf ref.fasta -r chr3:1,000-2,000 in1.bam in2.bam

merge 合并bam文件

samtools merge out.bam in1.bam in2.bam in3.bam

tview 图形化的界面查看mapping的结果

? -help

g -到达指定的染色体区域

samtools tview aln.sorted.bam ref.fasta

sort 将比对结果排序

samtools sort aln.bam aln.sorted

reheader 将bam文件的头文件替换。

index 为sort后的bam文件建立所以，这是在 view 比对结果时指定染色体区域必须做的。

samtools index aln.sorted.bam

rmdup 去除可能的PCR重复。将比对到相同位置的reads去重，仅保留比对质量最好的一个。一般只对paired end reads（fr）有效，并须指定ISIZE。

phase Call and phase heterozygous SNPs