癌症分析【第3课】测序下机数据处理（下）

第三课：从测序下机数据开始处理（下）

本篇将会介绍，如何把质控过的reads比对到参考基因组上，标记重复序列（Duplicates），并对比对文件进行质量值校正（BQSR）。

CNGBdb已经把动手实验室课程所用的文件放在ftp上啦！本课程使用的文件可从ftp://ftp.cngb.org/pub/Hands-on-Lab/Lesson3 上下载~

从ftp.cngb.org/pub/Hands-on-Lab找到课程文件

Step1 准备参考基因组

首先我们需要下载一个参考基因组文件。第1课热身有介绍过，如何使用google搜索找到UCSC的hg19参考基因组的fasta文件，但这次我们不下载完整的hg19，只需下载chr6.fa.gz文件（文件大小52 M）。在reads比对之前，我们需要先对参考基因组文件进行索引（index）。

思考题：如何找到UCSC的chr6.fa.gz参考基因组文件文件？

这一次我们使用的软件是BWA，是一个非常常用的比对工具。

BWA页面 http://bio-bwa.sourceforge.net/

• 下载bwa docker镜像

  docker pull biocontainers/bwa: v0.7.15_cv4

• 进入容器

docker run -it -v D:\CancerPipeline\L2_Reads_Clean:/data biocontainers/bwa:v0.7.15_cv4

• 解压参考基因组文件

gzip -dc chr6.fa.gz > chr6.fasta

• 建立参考基因组的索引

bwa index chr6.fasta

在bwa容器中建立参考基因组索引

Step2 用BWA比对

BWA有三种算法，BWA-BACKTRACK, BWA-SW和BWA-MEM，其中BWA-MEM是最新的（目前还没有发表文章），适合用于70bp到1Mbp长度的reads比对，而且速度快、性能好、准确性高。本次我们使用BWA-MEM算法进行比对。

对参考基因组建好了索引（Index）以后，我们可以开始BWA-MEM比对。比对时间较长（测试使用约5h），可以用笔记本插上电源，让程序跑过夜：

• 比对产生sam文件

bwa mem chr6.fasta V100007471_L03_539_1.clean.fq.gz V100007471_L03_539_2.clean.fq.gz > aln-pe.sam

思考题：演示只用了BWA-MEM默认参数。若要加上read groups信息（后文BQSR会使用到），可使用可选参数。阅读BWA-MEM参数文档，选择合适的参数，加上read groups信息。

比对完以后，还需要做一个很重要的步骤，处理重复序列（Duplicates），这一步主要是为了减少数据产生过程中（如，PCR扩增）带入的偏差。此处使用GATK的SortSam (Picard)和MarkDuplicates (Picard)工具。这些工具的具体用法，都可以在GATK Tool Documentaion中查询到哦！

从GATK Tool Documentation找到工具帮助文档

重复序列标记工具-MarkDuplicates

• 注：另外，在一些商业工具如Edico Genome和MegaBOLT的流程中，标记/去除重复序列可在比对的操作中同步进行。所以，如果你使用的比对工具MarkDuplicates/RemoveDuplicates的参数可设置，则可省略这个步骤。

• Convert sam to bam

samtools view -Sb aln-pe.sam > aln-pe.chr6.sort.bam

• Sort bam

gatk SortSam --java-options "-Xmx2g -XX:ParallelGCThreads=1 -Dsamjdk.compression_level=5 -XX:-UseGCOverheadLimit -Djava.io.tmpdir=/data/tmp" -I aln-pe.chr6.bam -O aln-pe.chr6.sort.bam -SO coordinate> gatk MarkDuplicates -I aln-pe.chr6.bam -O aln-pe.chr6.marked_duplicates.bam -M marked_dup_metrics.txt

• Mark Duplicates

gatk MarkDuplicates --java-options "-Xmx2g -XX:ParallelGCThreads=1 -Dsamjdk.compression_level=5 -XX:-UseGCOverheadLimit -Djava.io.tmpdir=/data/tmp" -I aln-pe.chr6.sort.bam -O marked_duplicates.bam -M marked_dup_metrics.txt

• 注：如果你运行gatk程序的时候，有内存溢出（out of memory）的报错信息，可以加上--java-options参数，设置内存限制（可选）。

Step3 GATK BQSR质量值校正

Sam/bam文件格式，每一列的意义可参考bwa官网SAM ALIGNMENT FORMAT部分内容（视频里也有介绍）。Sam/bam是一种存储比对reads的格式，了解每一列的意义，可以尝试将bam文件重新转成fastq格式（有时我们确实需要这样的操作），参考程序（Perl代码）在github上的地址：

https://github.com/Jemimacat/CancerPipeline/blob/master/scripts/bam2fastq.pl

Sam/bam flag计算工具：https://www.samformat.info/sam-format-flag

变异检测高度依赖碱基质量值，而质量值受很多因素影响，包括测序仪以及周围环境。GATK的BQSR (Base Quality Score Recalibration，碱基质量值校正)工具就是为解决这个问题而生的。GATK BQSR主要分两步：

1. 根据已知的变异集合，计算校正表格（Generates recalibration table for Base Quality Score Recalibration ，BQSR）。根据已知变异位点，计算每一个位点的经验错误概率（empirical probability of error），并输出到表格中。

碱基质量分数校正表格计算工具 BaseRecalibrator

1. 准备参考基因组.fai和.dict文件

   samtools fadix chr6.fasta> gatk CreateSequenceDictionary -R chr6.fasta

2. 添加read groups信息（可选）

   gatk AddOrReplaceReadGroups -I marked_duplicates.bam -O marked_duplicates.rg.bam --RGLB lib1 --RGPL mgiseq --RGPU unit1 --RGSM test

3. 碱基质量值校正表格计算

   gatk BaseRecalibrator -I marked_duplicates.rg.bam -R chr6.fasta --known-sites dbsnp_138.hg19.chr6.vcf.gz -O recal_data.table

4. 注：若bwa mem比对这一步没有加入read groups信息，需要在这一步用AddOrReplaceReadGroups (Picard)工具，加入read groups信息。

需要提前下载已知变异位点集合（.vcf），可从GATK Resource Bundle上下载：

从GATK Download找到资源（Resource Bundle）
2. 根据上一步的产生的表格，对bam文件进行质量值校正，并输出校正后的bam。这一步骤将校正测序仪引入的质量分数系统偏差。

校正碱基质量分数工具 ApplyBQSR

gatk ApplyBQSR -R chr6.fasta -I marked_duplicates.rg.bam --bqsr-recal-file recal_data.table -O bqsr.bam

这一步输出的校正后的bam文件（bqsr.bam），就可以用于后面的变异检测了。这一次的操作，bwa和GATK BQSR都会耗费一些时间和计算资源，稍微一点耐性哈~

延伸阅读

除了bwa+GATK的方案以外，目前还有很多加速分析的商业工具，例如，Illumina的Edico Genome，Sentieon和MGI的MegaBOLT等。这些商业工具把测序数据处理、比对、变异检测等设计成完整的流程，并采用不同的加速方案，使你的分析更流畅。

本次小课堂就到这里，下次再见~