美文网首页全基因组/外显子组测序分析高通量测序数据处理WES测序分析
GATK4最佳实践-数据预处理篇

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

作者: 生信修炼手册 | 来源:发表于2018-05-26 14:24 被阅读970次

    欢迎关注"生信修炼手册"！

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

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

    1. map to reference

    2. Mark Duplicates

    3. 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

    扫描关注微信号，更多精彩内容等着你！

    相关文章

      网友评论

      • 50e2e169eafe:--useOriginalQualities 命令 是使用原始的质量分数的意思么？ 为什么要用这个呢？

      本文标题：GATK4最佳实践-数据预处理篇

      本文链接：https://www.haomeiwen.com/subject/xlsmjftx.html

      延伸阅读

      深度阅读

      • 您也可以注册成为美文阅读网的作者，发表您的原创作品、分享您的心情！

      栏目导航

      热点阅读

      全基因组/外显子组测序分析

      高通量测序数据处理

      WES测序分析

      NGS

      单核苷酸多态性（SNP）

      关于我们|服务条款|联系我们|GATK4最佳实践-数据预处理篇|投稿指南|网站地图|RSS订阅|排版工具|手机版

      提供经典美文摘抄,优美散文欣赏,现代诗歌精选,短篇小说,心情随笔,表白情书范文,故事会在线阅读欣赏

      Copyright © 2014-2023 Haomeiwen.com All Rights Reserved. 好美文阅读网 版权所有

      备案信息：桂公网安备 45052102000051号 · 桂ICP备13007215号-3

      本站所收录作品、热点评论等信息部分来源互联网，目的只是为了系统归纳学习和传递资讯

      所有作品版权归原创作者所有，与本站立场无关，如不慎侵犯了你的权益，请联系我们告知，我们将做删除处理！