肿瘤外显子数据处理系列教程（四）比对结果的质控
大家好，肿瘤外显子数据处理系列教程将继续更新啦~
前面的系列教程，请看通知：肿瘤外显子数据处理系列教程阅读量太低-停更整改
回顾
前面我们已经走了比对了，拿到了bam文件，然而，在仔细检查bam文件（本文图表均以样本case1_biorep_A_techrep作为演示）的时候，发现其头文件有这样的一条记录：
@PG ID:bwa PN:bwa VN:0.7.17-r1198-dirty CL:bwa mem -M -t 11 -R @RG\tID:case1_biorep_A_techrep\tSM:case1_biorep_A_techrep\tLB:WXS\tPL:Illumina /home/llwu/reference/fasta/human/GRCh38 case1_biorep_A_techrep_1_val_1.fq.gz case1_biorep_A_techrep_1_val_1.fq.gz
这里是有问题的，比对时用到的两个fq文件是同一个，即case1_biorep_A_techrep_1_val_1.fq.gz
这是一个致命的错误，不过好在我们把每一步都记录下来的，报错我们可以从笔记中发现问题，前面的推文：肿瘤外显子数据处理系列教程（三）比对 记录了这个错误。所以大家实战时一定要做好记录！！！
GRCh38，在这里虽然可行，但是我们分析的是肿瘤全外显子测序的数据，后面要走GATK流程，而GATK要求的参考基因组比较特殊，需要从GATK的数据库下载，如果你不知道怎么下载，请参考这个教程：GATK4操作文档修改代码
重新修改代码后，提交任务，脚本如下：
# bwa.sh
cat ../case |while read id
do
fq1=${id}_1_val_1.fq.gz
fq2=${id}_2_val_2.fq.gz
echo $fq1
echo $fq2
bwa mem -M -t 11 -R "@RG\tID:$id\tSM:$id\tLB:WXS\tPL:Illumina" /public/biosoft/GATK/resources/bundle/hg38/bwa_index/gatk_hg38 $fq1 $fq2 | samtools sort -@ 11 -m 1G -o ../4.clean/${id}.bam -
done
最后，把任务挂到后台运行
nohup bash bwa.sh &
检查比对结果
同样的，我们也要对结果进行检查，代码在：肿瘤外显子数据处理系列教程（三）比对 ，和之前错误的结果进行一下比较，比如说bam文件的header中同一条记录，输入的fq文件应该是pair-end的fq1和fq2，即应该是这样的：
@PG ID:bwa PN:bwa VN:0.7.17-r1198-dirty CL:bwa mem -M -t 11 -R @RG\tID:case1_biorep_A_techrep\tSM:case1_biorep_A_techrep\tLB:WXS\tPL:Illumina /public/biosoft/GATK/resources/bundle/hg38/bwa_index/gatk_hg38 case1_biorep_A_techrep_1_val_1.fq.gz case1_biorep_A_techrep_2_val_2.fq.gz
然后就是对bam文件的一些统计信息，如用qualimap对case1_biorep_A_techrep的比对结果拿到bam文件进行统计，可以看到reads覆盖度条形图：
这和作者在文章中提到的基本吻合，在文章中关于测序深度的描述是：
还有基因组的覆盖信息，可以看到超过85%的位点的测序深度大于50x，不过需要明确一点，这里只是外显子区域的覆盖信息，因为测序策略就是WES，检查比对结果的时候也指定了外显子坐标文件/home/llwu/reference/index/human/CCDS/GRCh38.bed
然后就是测序reads的length：
也和文献中描述的差不多：
更让人无语的是，指出了这个学徒的错误，居然第二次仍然错！！！不可救药
下一期我们将通过igv来对这样的比对结果进行可视化以进一步了解bam文件，顺便学习一下igv这个强大的可视化工具，希望你不要错过，我们下周一再见。
实际上这个qualimap结果，如果样本比较多的话，那么报告就比较多，可以考虑multiqc软件可视化一下！
比如下面这个例子，就是数据覆盖度不够！
全国巡讲约你
第1-11站北上广深杭，西安，郑州， 吉林，武汉，成都，港珠澳（全部结束）
一年一度的生信技能树单细胞线下培训班（已结束）