GATK4最佳实践-数据预处理篇

生信修炼手册

发布于 2020-05-11 10:42:11

1.7K0

发布于 2020-05-11 10:42:11

GATK4 官方针对不同的变异类型，给出了好几套用于参考的pipeline。所有的pipeline有一个共同点，就是数据预处理部分。数据预处理的目的，是将原始的fastq或者ubam 文件，经过一系列处理，得到用于变异识别的bam文件，具体的示意图如下：

从示意图可以看出，预处理部分包含了3个主要步骤：

map to reference
Mark Duplicates
Base(Quality Score) Recalibration

之前版本的最佳实践都是给出了具体的gatk的代码，但是GATK4 给出的是用wdl这种workflow 语言编写的流程。对于预处理部分，对应的链接如下：

https://github.com/gatk-workflows/gatk4-data-processing

共给出了3套流程用于参考：

对于没有编程基础的人来说，理解一个wdl 脚本的逻辑是比较头疼的事情。我从3条流程中抽取出最关键的几个部分，整理如下：

1. map to reference

将原始序列和参考基因组进行比对, 官方推荐使用bwa比对软件，主要针对 DNA测序的数据。原始数据的格式为ubam。命令如下

java -jar picard.jar SamToFastq \
       INPUT=${input_bam}  \
       FASTQ=/dev/stdout  \
       INTERLEAVE=true  \
       NON_PF=true  | \
bwa mem -K 100000000 -p -v 3 -t 16 -Y ${bash_ref_fasta}  /dev/stdin - | \
java -jar picard.jar MergeBamAlignment \
    VALIDATION_STRINGENCY=SILENT \
   EXPECTED_ORIENTATIONS=FR \
     ATTRIBUTES_TO_RETAIN=X0 \
    ATTRIBUTES_TO_REMOVE=NM \
    ATTRIBUTES_TO_REMOVE=MD \
    ALIGNED_BAM=/dev/stdin \
    UNMAPPED_BAM=${input_bam} \
    OUTPUT=${output_bam_basename}.bam \
    REFERENCE_SEQUENCE=${ref_fasta} \
    PAIRED_RUN=true \
    SORT_ORDER="unsorted" \
    IS_BISULFITE_SEQUENCE=false \
    ALIGNED_READS_ONLY=false \
    CLIP_ADAPTERS=false \
    MAX_RECORDS_IN_RAM=2000000 \
    ADD_MATE_CIGAR=true \
    MAX_INSERTIONS_OR_DELETIONS=-1 \
    PRIMARY_ALIGNMENT_STRATEGY=MostDistant \
    PROGRAM_RECORD_ID="bwamem" \
    PROGRAM_GROUP_VERSION="${bwa_version}" \
    PROGRAM_GROUP_COMMAND_LINE="bwa mem -K 100000000 -p -v 3 -t 16 -Y ${bash_ref_fasta}" \
    PROGRAM_GROUP_NAME="bwamem" \
    UNMAPPED_READ_STRATEGY=COPY_TO_TAG \
    ALIGNER_PROPER_PAIR_FLAGS=true \
    UNMAP_CONTAMINANT_READS=true \
    ADD_PG_TAG_TO_READS=false

比对参考基因组部分共涉及到3个步骤，第一步picard将ubam转换成fastq; 第二步bwa 比对参考基因组；第三步picard将原始数据ubam和比对产生的aligned bam 合并，产生一个最终的bam文件。每个样本都对应生成一个这样的bam文件。

上面的代码中，有几个需要替换的变量。${input_bam}代表输入的原始序列的ubam格式的文件，${bash_ref_fasta}代表参考基因组序列的bwa 索引，${ref_fasta}代表参考基因组的fasta格式的序列；${output_bam_basename}代表最终输出的bam文件的名字，${bwa_version}代表bwa软件的版本号。

2. Mark Duplicates

标记bam文件中的重复序列，使用picard的MarkDuplicates命令，代码如下：

java -jar picard.jar \
      MarkDuplicates \
      INPUT=${input_bams} \
      OUTPUT=${output_bam_basename}.bam \
      METRICS_FILE=${metrics_filename} \
      VALIDATION_STRINGENCY=SILENT \
      OPTICAL_DUPLICATE_PIXEL_DISTANCE=2500 \
      ASSUME_SORT_ORDER="queryname" \
      CREATE_MD5_FILE=true

${input_bams}代表比对参考基因组之后生成的所有的bam文件，${output_bam_basename}代表产生的bam文件的名字，${metrics_filename} 代表产生的metrics文件的名称。

标记完重复序列之后，需要对产生的bam文件进行排序，命令如下

java -jar picard.jar \
   SortSam \
   INPUT=${input_bam} \
   OUTPUT=/dev/stdout \
   SORT_ORDER="coordinate" \
   CREATE_INDEX=false \
   CREATE_MD5_FILE=false \
   MAX_RECORDS_IN_RAM=300000 | \
java -jar picard.jar \
   SetNmAndUqTags \
   INPUT=/dev/stdin \
   OUTPUT=${output_bam_basename}.bam \
   CREATE_INDEX=true \
   CREATE_MD5_FILE=true \
   REFERENCE_SEQUENCE=${ref_fasta}

${input_bam}代表生成标记重复序列生成的bam文件；{$output_bam_basename}代表输出的排序之后的bam文件的名称，

3. Base(Quality Score) Recalibration

运行BQSR包括三步，第一步的命令如下：

gatk BaseRecalibrator \
   -R ${ref_fasta} \
   -I ${input_bam} \
   --use-original-qualities \
   -O ${recalibration_report_filename} \
   --known-sites ${dbSNP_vcf} \
   --known-sites ${sep=” —known-sites “ known_indels_sites_VCFs}

${ref_fasta}代表参考基因组的fasta序列；${input_bam}代表mark duplicates生成的排序好的bam文件；${recalibration_report_filename}代表产生的report文件的名字；${dbSNP_vcf} 代表NCBI SNP数据库中物种对应的vcf文件；${known_indels_sites_VCFs}代表其他数据库中已知的indel位点的vcf文件。

第二步的命令如下：

gatk GatherBQSRReports \
  -I ${sep=' -I ' input_bqsr_reports} \
  -O ${output_report_filename}

${input_bqsr_reports}代表第一步生成的report文件，每个样本一个，多个样本就指定多次 , 比如 -I 1.bqsr.report -I 2.bqsr.report; ${output_report_filename}代表输出的report文件的名字。

第三步的命令如下：

gatk ApplyBQSR \
   -R ${ref_fasta} \
   -I ${input_bam} \
   -O ${output_bam_basename}.bam \
   -bqsr ${recalibration_report} \
   --static-quantized-quals 10 --static-quantized-quals 20 --static-quantized-quals 30 \
   --add-output-sam-program-record \
   --create-output-bam-md5 \
   --use-original-qualities

${ref_fasta}代表参考基因组的fasta序列；${input_bam}代表mark duplicates生成的排序好的bam文件；{$output_bam_basename}代表输出的bam文件的名称; ${recalibration_report} 代表第二步生成的report文件。

当然你也可以直接运行别人的写好的wdl文件，以下面的这个预处理流程为例

https://github.com/gatk-workflows/gatk4-data-processing

首先得到整个脚本所有的参数列表，命令如下

java -jar womtools.jar inputs processing-for-variant-discovery-gatk4.wdl > processing_inputs.json

这个流程内置了一个hg38的参数模板，内容如下

{
 "##_COMMENT1": "SAMPLE NAME AND UNMAPPED BAMS",
 "PreProcessingForVariantDiscovery_GATK4.sample_name": "NA12878",
 "PreProcessingForVariantDiscovery_GATK4.ref_name": "hg38",
 "PreProcessingForVariantDiscovery_GATK4.flowcell_unmapped_bams_list": "gs://gatk-test-data/wgs_ubam/NA12878_24RG/NA12878_24RG_small.txt",
 "PreProcessingForVariantDiscovery_GATK4.unmapped_bam_suffix": ".bam", "##_COMMENT2": "REFERENCE FILES",
 "PreProcessingForVariantDiscovery_GATK4.ref_dict": "gs://broad-references/hg38/v0/Homo_sapiens_assembly38.dict",
 "PreProcessingForVariantDiscovery_GATK4.ref_fasta": "gs://broad-references/hg38/v0/Homo_sapiens_assembly38.fasta",
 "PreProcessingForVariantDiscovery_GATK4.ref_fasta_index": "gs://broad-references/hg38/v0/Homo_sapiens_assembly38.fasta.fai",
 "PreProcessingForVariantDiscovery_GATK4.SamToFastqAndBwaMem.ref_alt": "gs://broad-references/hg38/v0/Homo_sapiens_assembly38.fasta.64.alt",
 "PreProcessingForVariantDiscovery_GATK4.SamToFastqAndBwaMem.ref_sa": "gs://broad-references/hg38/v0/Homo_sapiens_assembly38.fasta.64.sa",
 "PreProcessingForVariantDiscovery_GATK4.SamToFastqAndBwaMem.ref_amb": "gs://broad-references/hg38/v0/Homo_sapiens_assembly38.fasta.64.amb",
 "PreProcessingForVariantDiscovery_GATK4.SamToFastqAndBwaMem.ref_bwt": "gs://broad-references/hg38/v0/Homo_sapiens_assembly38.fasta.64.bwt",
 "PreProcessingForVariantDiscovery_GATK4.SamToFastqAndBwaMem.ref_ann": "gs://broad-references/hg38/v0/Homo_sapiens_assembly38.fasta.64.ann",
 "PreProcessingForVariantDiscovery_GATK4.SamToFastqAndBwaMem.ref_pac": "gs://broad-references/hg38/v0/Homo_sapiens_assembly38.fasta.64.pac", "##_COMMENT3": "KNOWN SITES RESOURCES",
 "PreProcessingForVariantDiscovery_GATK4.dbSNP_vcf": "gs://broad-references/hg38/v0/Homo_sapiens_assembly38.dbsnp138.vcf",
 "PreProcessingForVariantDiscovery_GATK4.dbSNP_vcf_index": "gs://broad-references/hg38/v0/Homo_sapiens_assembly38.dbsnp138.vcf.idx",
 "PreProcessingForVariantDiscovery_GATK4.known_indels_sites_VCFs": [
   "gs://broad-references/hg38/v0/Mills_and_1000G_gold_standard.indels.hg38.vcf.gz",
   "gs://broad-references/hg38/v0/Homo_sapiens_assembly38.known_indels.vcf.gz"
 ],
 "PreProcessingForVariantDiscovery_GATK4.known_indels_sites_indices": [
   "gs://broad-references/hg38/v0/Mills_and_1000G_gold_standard.indels.hg38.vcf.gz.tbi",
   "gs://broad-references/hg38/v0/Homo_sapiens_assembly38.known_indels.vcf.gz.tbi"
 ],
 "##_COMMENT4": "MISC PARAMETERS",
 "PreProcessingForVariantDiscovery_GATK4.bwa_commandline": "bwa mem -K 100000000 -p -v 3 -t 16 -Y $bash_ref_fasta",
 "PreProcessingForVariantDiscovery_GATK4.compression_level": 5,
 "PreProcessingForVariantDiscovery_GATK4.SamToFastqAndBwaMem.num_cpu": "16",
 "##_COMMENT5": "DOCKERS",
 "PreProcessingForVariantDiscovery_GATK4.gotc_docker": "broadinstitute/genomes-in-the-cloud:2.3.0-1501082129",
 "PreProcessingForVariantDiscovery_GATK4.gatk_docker": "broadinstitute/gatk:4.0.0.0",
 "PreProcessingForVariantDiscovery_GATK4.picard_docker": "broadinstitute/genomes-in-the-cloud:2.3.0-1501082129",
 "PreProcessingForVariantDiscovery_GATK4.python_docker": "python:2.7",
 "##_COMMENT6": "PATHS",  
 "PreProcessingForVariantDiscovery_GATK4.gotc_path": "/usr/gitc/",
 "PreProcessingForVariantDiscovery_GATK4.picard_path": "/usr/gitc/",
 "PreProcessingForVariantDiscovery_GATK4.gatk_path": "/gatk/gatk",
 "##_COMMENT7": "JAVA OPTIONS",
 "PreProcessingForVariantDiscovery_GATK4.SamToFastqAndBwaMem.java_opt": "-Xms3000m",
 "PreProcessingForVariantDiscovery_GATK4.MergeBamAlignment.java_opt": "-Xms3000m",
 "PreProcessingForVariantDiscovery_GATK4.MarkDuplicates.java_opt": "-Xms4000m",
 "PreProcessingForVariantDiscovery_GATK4.SortAndFixTags.java_opt_sort": "-Xms4000m",
 "PreProcessingForVariantDiscovery_GATK4.SortAndFixTags.java_opt_fix": "-Xms500m",
 "PreProcessingForVariantDiscovery_GATK4.BaseRecalibrator.java_opt": "-Xms4000m",
 "PreProcessingForVariantDiscovery_GATK4.GatherBqsrReports.java_opt": "-Xms3000m",
 "PreProcessingForVariantDiscovery_GATK4.ApplyBQSR.java_opt": "-Xms3000m",
 "PreProcessingForVariantDiscovery_GATK4.GatherBamFiles.java_opt": "-Xms2000m",
 "##_COMMENT8": "MEMORY ALLOCATION",
 "PreProcessingForVariantDiscovery_GATK4.GetBwaVersion.mem_size": "1 GB",
 "PreProcessingForVariantDiscovery_GATK4.SamToFastqAndBwaMem.mem_size": "14 GB",
 "PreProcessingForVariantDiscovery_GATK4.MergeBamAlignment.mem_size": "3500 MB",
 "PreProcessingForVariantDiscovery_GATK4.MarkDuplicates.mem_size": "7 GB",
 "PreProcessingForVariantDiscovery_GATK4.SortAndFixTags.mem_size": "5000 MB",
 "PreProcessingForVariantDiscovery_GATK4.CreateSequenceGroupingTSV.mem_size": "2 GB",
 "PreProcessingForVariantDiscovery_GATK4.BaseRecalibrator.mem_size": "6 GB",
 "PreProcessingForVariantDiscovery_GATK4.GatherBqsrReports.mem_size": "3500 MB",
 "PreProcessingForVariantDiscovery_GATK4.ApplyBQSR.mem_size": "3500 MB",
 "PreProcessingForVariantDiscovery_GATK4.GatherBamFiles.mem_size": "3 GB",
 "##_COMMENT9": "DISK SIZE ALLOCATION",
 "PreProcessingForVariantDiscovery_GATK4.agg_small_disk": 200,
 "PreProcessingForVariantDiscovery_GATK4.agg_medium_disk": 300,
 "PreProcessingForVariantDiscovery_GATK4.agg_large_disk": 400,
 "PreProcessingForVariantDiscovery_GATK4.flowcell_small_disk": 100,
 "PreProcessingForVariantDiscovery_GATK4.flowcell_medium_disk": 200,
 "##_COMMENT10": "PREEMPTIBLES",
 "PreProcessingForVariantDiscovery_GATK4.preemptible_tries": 3,
 "PreProcessingForVariantDiscovery_GATK4.agg_preemptible_tries": 3
}

大部分的参数还是很好配置的，但是要注意一点，DOCKERS下配置的是软件对应的docker镜像，是需要安装的。

参数列表配置好之后，运行下面命令即可

java -jar Cromwell.jar run processing-for-variant-discovery-gatk4.wdl —inputs processing_inputs.json

本文参与腾讯云自媒体分享计划，分享自微信公众号。

原始发表：2018-05-17，如有侵权请联系 cloudcommunity@tencent.com 删除

processing

本文分享自生信修炼手册微信公众号，前往查看

如有侵权，请联系 cloudcommunity@tencent.com 删除。

本文参与腾讯云自媒体分享计划，欢迎热爱写作的你一起参与！

processing

登录后参与评论

0 条评论

热度