肿瘤外显子数据处理系列教程（四）比对结果的质控

大家好，肿瘤外显子数据处理系列教程将继续更新啦~

前面的系列教程，请看通知：肿瘤外显子数据处理系列教程阅读量太低-停更整改

回顾

前面我们已经走了比对了，拿到了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
这是一个致命的错误，不过好在我们把每一步都记录下来的，报错我们可以从笔记中发现问题，前面的推文：肿瘤外显子数据处理系列教程（三）比对 记录了这个错误。所以大家实战时一定要做好记录！！！

1566615907119

这并不是我们想要的结果，继续往下分析都是没有意义的。因此需要重新修改代码，而就在检查这个代码的同时，又发现了一个细节上的错误，参考基因组（索引）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覆盖度条形图：

genome_coverage_histogram

这和作者在文章中提到的基本吻合，在文章中关于测序深度的描述是：

1561381825372

还有基因组的覆盖信息，可以看到超过85%的位点的测序深度大于50x，不过需要明确一点，这里只是外显子区域的覆盖信息，因为测序策略就是WES，检查比对结果的时候也指定了外显子坐标文件/home/llwu/reference/index/human/CCDS/GRCh38.bed

genome_coverage_quotes

然后就是测序reads的length：

genome_insert_size_histogram

也和文献中描述的差不多：

1561382264619

当然了，在qualimap结果也可以看到学徒的这个错误：

更让人无语的是，指出了这个学徒的错误，居然第二次仍然错！！！不可救药

下一期我们将通过igv来对这样的比对结果进行可视化以进一步了解bam文件，顺便学习一下igv这个强大的可视化工具，希望你不要错过，我们下周一再见。

写在后面

实际上这个qualimap结果，如果样本比较多的话，那么报告就比较多，可以考虑multiqc软件可视化一下！

比如下面这个例子，就是数据覆盖度不够！

全国巡讲约你

一年一度的生信技能树单细胞线下培训班（已结束）