转录组数据分析笔记（3）——samtools用法小结

本篇笔记主要记录samtools的用法。

1. sam文件转化bam

bam文件是二进制文件，占用磁盘空间小，运算速度快，samtools操作是针对bam文件的，所以我们要进行数据转化。samtools sort指令可以将bam文件进行排序，这个指令同时也可以将sam文件转化成bam文件：

#!/bin/bash
ls *.sam | while read id
do
    samtools sort -l 0 -@ 5 -o $(basename $id ".sam").bam $id     # 指定输出文件，改后缀.bam
done

运行脚本，将当前目录的sam文件全转换成bam文件并排序（这里的排序不是按名称，用HTseq还要再按照read名称排序，使用参数-n）。

samtools sort参数

samtools sort # 对bam文件进行排序（sam文件排序不会变）
-l # 设置输出文件压缩等级，0-9，0是不压缩,9是最高等级压缩
-@ # 设置线程数
-o # 设置排序后输出的文件名
最后接输入的bam或者sam格式文件

2. 构建索引文件

2.1 构建bam文件索引

在bam文件目录下，排序后的bam文件可以建立索引：

$ ls *.bam | xargs -i samtools index {}

注意下xargs -i的用法，和管道不一样，是传递参数给后一个命令的花括号中，后一个命令中不存在歧义的时候可省略参数-i和花括号。

如图生成的bai文件就是索引文件。其实到了这一步，前面的sam文件就可以删除（节省电脑空间），只留下bam文件就行，bam文件无法直接查看，可以通过samtools view命令查看bam文件。

2.2 构建参考基因组fa文件索引

在参考基因组文件目录下，对参考基因组的fa文件建立索引：

$ samtools faidx Arabidopsis_thaliana.dna.genome.fa

参考基因组文件名注意改成自己的，生成的索引文件是.fai结尾的

3. bam文件qc质控

samtools转化生成的bam文件需要进行质控，看看比对情况如何。在bam文件目录下，我们创建一个samtools自带qc质控指令samtools flagstat运行脚本：

#!/bin/bash
ls *.bam | while read id
do
    samtools flagstat -@ 4 $id > $(basename $id ".bam").flagstat     # 自定义输出文件
done

$ samtools flagstat bam文件 > 输出文件 # 这种格式，其他参数都一样

运行脚本文件可以获得16个.flagstat质控文件，和fastqc一样，我们还可以做完后用multiqc命令集合成一个html格式的总的qc报告网页。和fastqc不同之处是，fastqc是做下机数据质控，samtools是做比对参考基因组的质控。如下图所示，可以比较直观地看出大部分reads都是map上的。

生成的每一个flagstat文件我们也可以直接点开。

每一行统计数据都是以通过QC的reads数量和未通过QC的reads数量组成，以我点开的这个文件为例，主要信息有以下几个：

13992629个reads都是合格的
12328290个reads只比对到参考基因组一个位置上
13988737个reads比对到参考基因组（99.97%）
12332182个reads是成对的
12201338个reads可以正确配对（98.94%）
2846条reads成对但只有一条能比对上参考基因组
12398个配对的reads可以比对到别的染色体上

可以自己将所有的flagstat运行结束后的文件放在一个目录下，运用paste命令全部按列粘贴在一起，用cut或者awk提取所需的列数据自己做比对情况表格，这里不再赘述。