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）》。

    

    

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