背景信息

生命科學領域內基因測序技術的飛速發展，使得人類發現的基因序列以指數級增長，對于如此數量龐大的基因進行同源性搜尋、比對、變異檢查等，往往伴隨著巨大的數據處理和并行計算。高性能計算可以提供強大的算力支持，使用多種調度器提高并發效率，使用GPU進行計算加速等。

本文以經典及普及的二代全基因組測序WGS（Whole Genome Sequencing）流程為例，結合二代測序軟件GATK，介紹人類全基因組測序的通用流程。在實際生信分析中，需要結合不同的業務需求進行相應的調整，如增加變異檢測質控及過濾等，您可以結合實際場景進行調整。

在進行基因測序時，您可以使用BWA構建索引及比對記錄，再使用Samtools對比對記錄進行排序，然后使用GATK去除重復序列、重新校正堿基質量值、變異檢查。

• BWA（Burrows-Wheeler-Alignment Tool）是一款將DNA序列映射到參考基因組上的軟件，例如比對人類基因組。

• GATK（The Genome Analysis Toolkit）是一款二代重測序數據分析軟件，是基因分析的工具集。主要用于去除重復序列、重新校正堿基質量值、變異檢查等。

• Samtools是用于處理sam和bam格式的工具軟件，能夠查看二進制文件、轉換文件格式、對文件排序及合并，可以結合sam格式文件中的flag、tag等信息，對比對結果進行統計匯總。

準備工作

1. 創建E-HPC集群。具體操作，請參見使用向導創建集群。

   配置集群時，軟硬件參數配置如下：

   參數	說明
   硬件參數	部署方式為精簡，包含1個管控節點和1個計算節點，規格如下： 管控節點：采用ecs.c7.large實例規格，該規格配置為2 vCPU，4 GiB內存。 計算節點：采用ecs.g7.16xlarge實例規格，該規格配置為64 vCPU、256 GiB內存。
   軟件參數	鏡像選擇CentOS 7.6公共鏡像，調度器選擇slurm，打開VNC開關。

2. 創建集群用戶。具體操作，請參見創建用戶。

   集群用戶用于登錄集群，進行編譯軟件、提交作業等操作。本文創建的用戶示例如下：

   • 用戶名：testuser

   • 用戶組：sudo權限組

3. 安裝軟件。

   待安裝軟件詳情如下：

   軟件名稱	版本	下載地址
   BWA	4.2.0.0	BWA
   GATK	0.7.17	GATK
   Samtools	1.14	Samtools

步驟一：連接集群

1. 選擇以下一種方式連接集群：

   • 通過客戶端

     該方式僅支持使用PBS調度器的集群。操作前，請確保您已下載安裝E-HPC客戶端，且已配置客戶端所需環境。具體操作，請參見配置客戶端所需環境。

     1. 打開并登錄E-HPC客戶端。

     2. 在客戶端左側導航欄，單擊會話管理。

     3. 在會話管理頁面的右上角，單擊terminal，打開Terminal窗口。

   • 通過控制臺

     1. 登錄彈性高性能計算控制臺。

     2. 在頂部菜單欄左上角處，選擇地域。

     3. 在左側導航欄，單擊集群。

     4. 在集群頁面，找到目標集群，單擊遠程連接。

     5. 在遠程連接頁面，輸入集群用戶名、登錄密碼和端口，單擊ssh連接。

2. 執行以下命令，創建作業執行目錄，本文以/home/data/human為例。

   說明

   為確保有足夠的存儲空間用于下載和處理數據，建議將作業執行目錄配置在集群掛載的共享存儲目錄/home下。

   sudo mkdir -p /home/data/human

步驟二：下載數據并進行預處理

1. 執行以下命令，創建并打開作業腳本文件，腳本文件命名為data_pre.slurm。

   sudo vim data_pre.slurm

   腳本內容示例如下：

   #!/bin/bash -l
   #SBATCH -J data_pre
   #SBATCH --partition comp
   #SBATCH -N 1
   #SBATCH -n 1
   #SBATCH --output=data_pre.log
   
   cd /home/data/human/
   
   # 下載并解壓人類參考基因組、人類基因測序樣本1、人類基因測序樣本2
   wget https://public-ehpc-package.oss-cn-hangzhou.aliyuncs.com/lifescience/b37_human_g1k_v37.fasta
   wget https://public-ehpc-package.oss-cn-hangzhou.aliyuncs.com/lifescience/b37_1000G_phase1.indels.b37.vcf
   wget https://public-ehpc-package.oss-cn-hangzhou.aliyuncs.com/lifescience/b37_Mills_and_1000G_gold_standard.indels.b37.vcf
   wget https://public-ehpc-package.oss-cn-hangzhou.aliyuncs.com/lifescience/b37_dbsnp_138.b37.vcf.zip
   unzip b37_dbsnp_138.b37.vcf.zip
   
   wget https://public-ehpc-package.oss-cn-hangzhou.aliyuncs.com/lifescience/gatk-examples_example1_NA20318_ERR250971_1.filt.fastq.gz
   gunzip gatk-examples_example1_NA20318_ERR250971_1.filt.fastq.gz
   
   wget https://public-ehpc-package.oss-cn-hangzhou.aliyuncs.com/lifescience/gatk-examples_example1_NA20318_ERR250971_2.filt.fastq.gz
   gunzip gatk-examples_example1_NA20318_ERR250971_2.filt.fastq.gz
   
   # 使用bgzip對文件進行壓縮
   bgzip -f b37_1000G_phase1.indels.b37.vcf
   bgzip -f b37_Mills_and_1000G_gold_standard.indels.b37.vcf
   bgzip -f b37_dbsnp_138.b37.vcf
   
   # 使用tabix建立索引
   tabix -p vcf  b37_1000G_phase1.indels.b37.vcf.gz
   tabix -p vcf  b37_Mills_and_1000G_gold_standard.indels.b37.vcf.gz
   tabix -p vcf  b37_dbsnp_138.b37.vcf.gz

2. 執行以下命令，提交作業。

   sbatch data_pre.slurm

   返回示例如下，表示已成功提交作業。

   

3. 執行以下命令，查看作業狀態。

   squeue

4. 待作業執行完成后，執行以下命令，查看data_pre.log文件結果。

   cat data_pre.log

步驟三：構建參考基因組索引

1. 執行以下命令，創建并打開作業腳本文件，腳本文件命名為bwa_index.slurm。

   sudo vim bwa_index.slurm

   腳本內容示例如下：

   #!/bin/bash -l
   #SBATCH -J bwa_index
   #SBATCH --partition comp
   #SBATCH -N 1
   #SBATCH -n 1
   #SBATCH --output=bwa_index.log
   
   cd /home/data/human/
   
   # 構建BWA比對所需的參考基因組的index數據庫
   time bwa index b37_human_g1k_v37.fasta
   
   # 創建fasta序列格式索引
   samtools faidx b37_human_g1k_v37.fasta

2. 執行以下命令，提交作業。

   sbatch --dependency=afterok:jobid1 bwa_index.slurm

3. 執行以下命令，查看作業狀態。

   squeue

4. 待作業執行完成后，執行以下命令，查看bwa_index.log文件結果。

   cat bwa_index.log

步驟四：執行基因測序分析

1. 執行以下命令，創建并打開作業腳本文件，腳本文件命名為wgs_main_process.slurm。

   sudo vim wgs_main_process.slurm

   腳本內容示例如下：

   重要

   請將第六行#SBATCH -w compute000字段中compute000內容替換為集群實際運行節點。

   #!/bin/bash -l
   #SBATCH -J wgs_main_process
   #SBATCH --partition comp
   #SBATCH -N 1
   #SBATCH -n 64
   #SBATCH -w compute000
   #SBATCH --output=wgs_main_process.log
   
   cd /home/data/human/
   
   # 將樣本測序數據reads與人類參考基因組進行比對，并將輸出文件轉化為bam格式，以節省磁盤空間
   time bwa mem -t 64 \
   -R '@RG\tID:ehpc\tPL:illumina\tLB:library\tSM:b37' \
   b37_human_g1k_v37.fasta \
   gatk-examples_example1_NA20318_ERR250971_1.filt.fastq \
   gatk-examples_example1_NA20318_ERR250971_2.filt.fastq | samtools view -S -b - > ERR250971.bam
   
   # 將比對記錄按照順序從小到大進行排序
   time samtools sort -@ 64 -O bam -o ERR250971.sorted.bam ERR250971.bam
   
   # 去除DNA序列PCR擴增產生的重復reads序列
   time gatk MarkDuplicates \
   -I ERR250971.sorted.bam -O ERR250971.markdup.bam \
   -M ERR250971.sorted.markdup_metrics.txt
   
   # 構建markdup測序bam文件索引
   samtools index ERR250971.markdup.bam
   
   # 生成參考基因組dict文件
   time gatk CreateSequenceDictionary  \
   -R b37_human_g1k_v37.fasta \
   -O b37_human_g1k_v37.dict
   
   # 利用已有的snp及indels數據庫，建立相關模型，構建重校準表
   time gatk BaseRecalibrator \
   -R b37_human_g1k_v37.fasta \
   -I ERR250971.markdup.bam  \
   --known-sites b37_1000G_phase1.indels.b37.vcf.gz \
   --known-sites b37_Mills_and_1000G_gold_standard.indels.b37.vcf.gz \
   --known-sites b37_dbsnp_138.b37.vcf.gz -O ERR250971.BQSR.table
   
   # 對原始堿基質量值進行調整
   time gatk ApplyBQSR \
   --bqsr-recal-file ERR250971.BQSR.table \
   -R b37_human_g1k_v37.fasta \
   -I ERR250971.markdup.bam \
   -O ERR250971.BQSR.bam
   
   # 構建BQSR測序bam文件索引
   time samtools index ERR250971.BQSR.bam
   
   # 輸出變異檢測結果vcf文件
   time gatk HaplotypeCaller \
   -R b37_human_g1k_v37.fasta  \
   -I ERR250971.BQSR.bam \
   -O ERR250971.HC.vcf.gz

2. 執行以下命令，提交作業。

   sbatch --dependency=afterok:jobid2 wgs_main_process.slurm

3. 執行以下命令，查看作業狀態。

   squeue

4. 待作業執行完成后，執行以下命令，查看wgs_main_process.log文件結果。

   cat wgs_main_process.log

測序耗時說明

本文使用的計算節點為ecs.g7.16xlarge，僅為演示全基因組測序流程，整體大約耗時6~7小時，下表展示了基因測序各流程大約使用的時間供您參考，并非最優性能。