文档HOT
中文
注册
我要评分
文档获取效率
文档正确性
内容完整性
文档易理解
提交
在线提单
论坛求助

GATK4 Spark模式运行

前提条件

准备数据文件，GATK测试中需要两类数据文件：参考基因组数据文件（FASTA格式文件）和原始测序数据（FASTQ或BAM格式文件）。

在本测试算例中使用的数据文件如下：

  • 基因组数据文件（FASTA格式文件）

    以下为人类基因组信息。

    human_g1k_v37.fasta

    以下两个vcf来自于千人基因组和Mills项目，里面记录了那些项目中检测到的人群Indel区域。
    • Mills_and_1000G_gold_standard.indels.hg19.sites.vcf.gz
    • 1000G_phase1.indels.hg19.sites.vcf.gz

    以下文件汇集的是目前几乎所有的公开人群变异数据集。

    dbsnp132_20101103.vcf

    变异数据集根据测序目的选择一到三种数据。

  • 原始测序数据（FASTQ或BAM格式文件）
    • SRR742200_1.fastq
    • SRR742200_2.fastq
    1. 使用PuTTY工具，以root用户登录服务器。
    2. 执行以下命令创建“input”文件夹。
      mkdir projectDir
cd projectDir
mkdir input
    3. 执行以下命令解压算例文件。
      gzip -d SRR742200_1.fastq.gz
gzip -d SRR742200_1.fastq.gz
gzip -d human_g1k_v37.fasta.gz
gzip -d dbsnp132_20101103.vcf.gz
    4. 执行以下命令使用blat生成“2bit”格式文件。
      faToTwoBit human_g1k_v37.fasta human_g1k_v37.fasta.2bit
    5. 执行以下命令将所有数据文件放进“input”文件夹内。
      cp human_g1k_v37.fasta human_g1k_v37.fasta.2bit SRR742200_1.fastq SRR742200_2.fastq dbsnp132_20101103.vcf input

操作步骤

  1. 使用PuTTY工具，以root用户登录服务器。
  2. 执行以下命令，BWA中“index”对“fasta”文件构建索引。
    bwa index -a bwtsw human_g1k_v37.fasta

    回显显示如下信息：

    [bwa_index] Pack FASTA... 27.89 sec
[bwa_index] Construct BWT for the packed sequence...
[BWTIncCreate] textLength=6203609478, availableWord=448508744
[BWTIncConstructFromPacked] 10 iterations done. 99999990 characters processed.
[BWTIncConstructFromPacked] 20 iterations done. 199999990 characters processed.
[BWTIncConstructFromPacked] 30 iterations done. 299999990 characters processed.
[BWTIncConstructFromPacked] 40 iterations done. 399999990 characters processed.
[BWTIncConstructFromPacked] 50 iterations done. 499999990 characters processed.
[BWTIncConstructFromPacked] 60 iterations done. 599999990 characters processed.
[BWTIncConstructFromPacked] 70 iterations done. 699999990 characters processed.
[BWTIncConstructFromPacked] 80 iterations done. 799999990 characters processed.
[BWTIncConstructFromPacked] 90 iterations done. 899999990 characters processed.
<以下输出内容省略>
  3. 执行以下命令，查看生成索引文件。
    ls

    完成2和3后，会在当前路径查看到“human_g1k_v37.fasta.amb”、“human_g1k_v37.fasta.ann”、“human_g1k_v37.fasta.bwt”、“human_g1k_v37.fasta.pac”、“human_g1k_v37.fasta.sa”5个文件。生成的索引文件各pipeline通用，仅需构建一次，时间比较长，需要1个半小时左右，将会在后续测序步骤中使用生成的文件。

  4. 执行以下命令，BWA中“mapping”对“fastq”文件进行比对。
    1. 生成字典。
      gatk CreateSequenceDictionary -R human_g1k_v37.fasta -O human_g1k_v37.dict

      回显显示如下红框中所示信息，表示执行完成。

      生成的“dict”文件是为后面进行基因组的比对。

    2. 进行比对。
      bwa mem -M -t 96 human_g1k_v37.fasta SRR742200_1.fastq SRR742200_2.fastq > SRR7.sam
      表1 参数说明

      参数

      说明

      -t

      参数为线程数，一般与服务器核数一致。

      回显显示如下红框中所示信息，表示执行完成。

  5. 执行以下命令，Reorder对基因序列进行重新排序。

    BWA生成的“sam”文件默认按字典式排序，需要转换为按染色体组型进行排序，由picard工具完成，GATK4中已集成picard功能。

    gatk ReorderSam -I SRR7.sam -O SRR7_reorder.bam -R human_g1k_v37.fasta

    回显显示如下红框中所示信息，表示执行完成。

  6. Sort对基因序列按照坐标顺序从大到小排序。
    1. Spark运行于HDFS文件系统中，在linux系统下不可见，执行以下命令将“SRR7_reorder.bam”文件传入HDFS文件系统。
      hadoop fs -mkdir /seqdata
hadoop fs -put SRR7_reorder.bam /seqdata
    2. 执行以下命令确认“sparklog”目录已手动创建。
      hadoop fs -mkdir /sparklog
hadoop fs -ls /
    3. 执行以下命令进行排序。
      gatk SortReadFileSpark -I hdfs://$(hostname):9000/seqdata/SRR7_reorder.bam -O hdfs://$(hostname):9000/seqdata/SRR7_sorted.bam -- --spark-runner SPARK --spark-master spark://$(hostname):7077

      根据Hadoop和Spark所布置的master节点来修改此处的Hostname节点，详细请参考《Hadoop 3.1.2+Spark 2.4.4 移植指南（CentOS 7.6）》。

      回显显示如下红框中所示信息，表示执行完成。

    4. 执行以下命令将生成的文件拷贝到数据目录。
      hadoop fs -get /seqdata/SRR7_sorted.bam ./
    5. 执行以下命令查看数据目录。
      ll

      回显显示如下红框中所示信息，表示执行完成。

  7. 执行以下命令，Add head对基因序列进行加头处理。

    GATK2.0以上已不支持无头文件检测，该步骤可以在BWA比对时加-r进行，也可以使用GATK4中AddOrReplaceReadGroups工具完成。

    gatk AddOrReplaceReadGroups -I SRR7_sorted.bam -O SRR7_header.bam -LB lib1 -PL illumina -PU unit1 -SM 20

    回显显示如下红框中所示信息，表示执行完成。

  8. Mark Duplicates对基因序列去重。
    1. 执行以下命令将加头文件传入HDFS文件系统。
      hadoop fs -put SRR7_header.bam /seqdata
    2. 执行以下命令查看HDFS文件系统。
      hadoop fs -ls /seqdata

    3. 执行以下命令进行去重操作。
      gatk MarkDuplicatesSpark -I hdfs://$(hostname):9000/seqdata/SRR7_header.bam -O hdfs://$(hostname):9000/seqdata/SRR7_markdup.bam -- --spark-runner SPARK --spark-master spark://$(hostname):7077

      根据Hadoop和Spark所布置的master节点来修改此处的Hostname节点，详细请参考《Hadoop 3.1.2+Spark 2.4.4 移植指南（CentOS 7.6）》。

      回显显示如下红框中所示信息，表示执行完成。

  9. BQSR进行重新校准。
    1. 执行以下命令将原始文件传入HDFS文件系统。
      hadoop fs -put dbsnp132_20101103.vcf /seqdata
hadoop fs -put human_g1k_v37.fasta.2bit /seqdata
hadoop fs -ls /seqdata

      回显显示如红框中所示信息，表示执行完成。

    2. 执行BQSR（Base Quality Score Recalibration）。
      gatk BQSRPipelineSpark -I hdfs://$(hostname):9000/seqdata/SRR7_markdup.bam -R hdfs://$(hostname):9000/seqdata/human_g1k_v37.fasta.2bit -O hdfs://$(hostname):9000/seqdata/SRR7_bqsr.bam --known-sites hdfs://$(hostname):9000/seqdata/dbsnp132_20101103.vcf --disable-sequence-dictionary-validation true -- --spark-runner SPARK --spark-master spark://$(hostname):7077

      根据Hadoop和Spark所布置的master节点来修改此处的Hostname节点，详细请参考《Hadoop 3.1.2+Spark 2.4.4 移植指南（CentOS 7.6）》。

      显示如下红框中所示信息表示执行完成。

  10. 执行以下命令，进行HaplortypeCaller变异流程检测，生成变异体“vcf”文件。
    gatk HaplotypeCallerSpark -R hdfs://$(hostname):9000/seqdata/human_g1k_v37.fasta.2bit -I hdfs://$(hostname):9000/seqdata/SRR7_bqsr.bam -O hdfs://$(hostname):9000/seqdata/SRR7.snp.raw.vcf -- --spark-runner SPARK --spark-master spark://$(hostname):7077

    根据Hadoop和Spark所布置的master节点来修改此处的Hostname节点，详细请参考《Hadoop 3.1.2+Spark 2.4.4 移植指南（CentOS 7.6）》。

    最终生成的数据信息如下所示。

父主题： 运行和验证