本篇将会介绍,如何把质控过的reads比对到参考基因组上,标记重复序列(Duplicates),并对比对文件进行质量值校正(BQSR)。 CNGBdb已经把动手实验室课程所用的文件放在ftp上啦!本课程使用的文件可从ftp://ftp.cngb.org/pub/Hands-on-Lab/Lesson3 上下载~ 从ftp.cngb.org/pub/Hands-on-Lab找到课程文件 首先我们需要下载一个参考基因组文件。第1课热身有介绍过,如何使用google搜索找到UCSC的hg19参考基因组的fasta文件,但这次我们不下载完整的hg19,只需下载chr6.fa.gz文件(文件大小52 M)。在reads比对之前,我们需要先对参考基因组文件进行索引(index)。 思考题 如何找到UCSC的chr6.fa.gz参考基因组文件文件? 这一次我们使用的软件是BWA,是一个非常常用的比对工具。 BWA页面 http://bio-bwa./ ## 下载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容器中建立参考基因组索引 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参数,设置内存限制(可选)。 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./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 *注:若bwa mem比对这一步没有加入read groups信息,需要在这一步用AddOrReplaceRead Groups (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都会耗费一些时间和计算资源,稍微一点耐性哈~ |
|