科普

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

2019-07-11 137动手实验室

第三课:从测序下机数据开始处理(下)

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

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


癌症分析【第3课】测序下机数据处理(下)_1

从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,是一个非常常用的比对工具。

癌症分析【第3课】测序下机数据处理(下)_2

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

癌症分析【第3课】测序下机数据处理(下)_3

在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中查询到哦!

癌症分析【第3课】测序下机数据处理(下)_4

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

癌症分析【第3课】测序下机数据处理(下)_5

重复序列标记工具-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上的地址:

癌症分析【第3课】测序下机数据处理(下)_6

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

癌症分析【第3课】测序下机数据处理(下)_7

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),并输出到表格中。

癌症分析【第3课】测序下机数据处理(下)_8

碱基质量分数校正表格计算工具 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上下载:

癌症分析【第3课】测序下机数据处理(下)_9

从GATK Download找到资源(Resource Bundle)

2. 根据上一步的产生的表格,对bam文件进行质量值校正,并输出校正后的bam。这一步骤将校正测序仪引入的质量分数系统偏差。
癌症分析【第3课】测序下机数据处理(下)_10

校正碱基质量分数工具 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等。这些商业工具把测序数据处理、比对、变异检测等设计成完整的流程,并采用不同的加速方案,使你的分析更流畅。

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

上一篇下一篇