166-为什么使用samtools检查bam质量时报错？

Feb 14, 2020 3 min read cnposts

刘小泽写于2020.2.14 在做数据比对后，会得到bam文件，那么这个比对的结果如何呢？需要检查一下，检查的方法一般是用samtools自带的命令，另外还有用qualimap。我们今天来看看容易报错的samtools 命令

首先需要得到一个统计结果

一般在进行统计这一步之前，比对好的bam会进行sort和index，下面这个脚本是为了应付多样本循环操作的，只需要制定好fq1、fq2的路径，以及对应的name，就能批量操作多样本比对，并完成排序和构建索引

# INDEX就是bwa构建好的索引文件
bwa mem -M -t 4 -R "@RG\tID:${name}\tSM:${name}\tLB:WXS\tPL:Illumina" ${INDEX} ${fq1} ${fq2} | samtools sort -@ 4 -o ${name}.sort.bam - && samtools index ${name}.sort.bam

有了bam文件后，先使用samtools stats根据bam得到统计结果

REF=~/reference/Homo_sapiens_assembly38.fasta
samtools stats -@ 4 --reference $REF ${name}.sort.bam > ${name}.stat

当然这个结果我们不需要肉眼去看，可以再利用multiqc 去得到一个总体比对的表格

multiqc *stat

把multiqc_data和multiqc_report.html导出来，直接打开html查看

可能报错的地方在于plot-bamstats

samtools对bam的检验方便就方便在于，它能直接根据上一步的${name}.stat 进行画图

这一点，在官网介绍就提到了：https://www.htslib.org/doc/samtools.html【其中包含了samtools各个功能的介绍】

其实这个samtooll内置的命令很简单，一行就搞定，只需要指定一个文件前缀，然后提供.stat文件

plot-bamstats -p ${name} ${name}.stat

但是，很有可能使用conda安装的samtools，会因为版本导致这个功能不能用

第一个错误：缺少 gnuplot

因为plot-bamstats 需要调用gnuplot去画图

conda install -c bioconda gnuplot -y

第二个错误：缺少某个库

比如报错：

gnuplot: error while loading shared libraries: libwebp.so.6: cannot open shared object file: No such file or directory

搜索便知，我们同样可以用conda完成安装

conda install libwebp=0.5.2 -y

第三个错误：软件版本

可能琢磨好半天，也不如降级来的快一些

对于samtools，使用一个更稳定的版本或许真的是个好主意。

conda默认安装的最新版本是1.9（官网提供了1.10源代码安装），如果别人能成功作图，而自己却不能，就要考虑，是不是版本不一致导致的。并不是版本越新越好！

那就试试1.8版本

conda create -n old samtools=1.8 gnuplot
conda activate old

判断标准

其实运行完，应该能看到生成html，如果没有那么就要从上面找找原因：

另外，除了samtools，还可以用qualimap去质控

其实网上搜索，会出来很多的结果，这里我也是把默认的代码用法放一下，方便查阅

qualimap bamqc --java-mem-size=10G -gff $REF -nr 1000 -nw 400 -nt 10 -bam ${name}.sort.bam -outdir ${name}.quanlimap.stat

	# -nr: Number of reads analyzed in a chunk (default 1000)
	# -nw: Number of windows (default 400)
	# -nt: threads (default 48)

值得注意的是，这个-gff 指定的参考文件，取决于测序数据类型

外显子：需要自己根据CCDS文件制作exon.bed
转录组：使用gencode的GTF文件即可

具体怎么操作，看这篇就好啦：https://cloud.tencent.com/developer/article/1168326

samtools