Variant calling using GATK - Biomedical Data Science Laboratory

Introduction

Germline and somatic variants can reveal the genetic factors underlying a disease or phenotype. These variants are identified by mapping sequencing reads to a reference genome and detecting differences. However, it is necessary to distinguish disease-associated variants from natural genetic variation, such as single nucleotide polymorphisms (SNPs), which are widely distributed across population rather than enriched in specific disease groups.

Therefore, accurate identification of phenotype-associated variants requires not only read mapping but also the incorporation of external knowledge and filtering strategies. Mapping to a reference genome followed by systematic filtering enables the detection of candidate variants. Because these procedures are well-defined, they can be standardized to improve reproducibility and analytical convenience.

In this exercise, we will use a GATK workflow to perform variant discovery. This material is intended for a graduate-level class in the Department of Molecular Biology at Jeonbuk National University and is provided for educational purposes only. The code examples are designed specifically for the exercise server and are not recommended for direct use in real research settings.

Installation of gatk using conda

Create a local directory on the exercise server to install gatk. For manual execution, specify the path to the gatk executable, for example ~/local/gatk-4.3.0.0/gatk.

cd ~
mkdir –p local/gatk 
cd local/gatk
wget https://github.com/broadinstitute/gatk/releases/download/4.3.0.0/gatk-4.3.0.0.zip
unzip gatk-4.3.0.0.zip
cd gatk-4.3.0.0
conda env create –n gatk –f gatkcondaenv.yml

Copy the exercise directory

In this activity, we will use a GATK workflow, which requires several input files. First, inspect the file system in the gatk-workflow directory. You will find FASTQ files that are part of a publicly available dataset.

cd ~
cp /home/shared/gatk-workflows/ . -r

The source of the FASTQ input files
-       https://zenodo.org/record/2582555/files/SLGFSK-N_231335_r1_chr5_12_17.fastq.gz
        https://zenodo.org/record/2582555/files/SLGFSK-N_231335_r2_chr5_12_17.fastq.gz
        https://zenodo.org/record/2582555/files/SLGFSK-T_231336_r1_chr5_12_17.fastq.gz
        https://zenodo.org/record/2582555/files/SLGFSK-T_231336_r2_chr5_12_17.fastq.gz

Step 1. Convert FASTQ to unmapped SAM file

Paired short read data files (FASTQ) can be merged into a single SAM file.

The script (1_fastq2usam.sh) can be found in the gatk-workflow directory. You do not need to execute the code as the bam file was created in advance.

Script: 1_fastq2usam.sh

~/local/gatk-4.3.0.0/gatk FastqToSam -F1 inputs/1M/normal_1M_r1.fastq -F2 inputs/1M/normal_1M_r2.fastq -SM n1 -RG HT33CBBXX -PL ILLUMINA -O inputs/normal.1M.unmapped.bam
~/local/gatk-4.3.0.0/gatk FastqToSam -F1 inputs/1M/tumor_1M_r1.fastq -F2 inputs/1M/tumor_1M_r2.fastq -SM t1 -RG HT33CBBXX -PL ILLUMINA -O inputs/tumor.1M.unmapped.bam

Step 2. Align to the reference genome (normal tissue)

The sequencing reads will be mapped to the human reference genome. This process internally uses BWA, and therefore the corresponding BWT index files for the reference genome must be provided. The input data configuration is defined in the file gatk-workflows/gatk-workflows/gatk4-data-processing/processing-for-variant-discovery-gatk4.hg38.wgs.normal.inputs.json.

Script: 2_gatk4_data_processing_normal.sh

java -jar gatk-workflows/cromwell-47.jar run gatk-workflows/gatk4-data-processing/processing-for-variant-discovery-gatk4.wdl --inputs gatk-workflows/gatk4-data-processing/processing-for-variant-discovery-gatk4.hg38.wgs.normal.inputs.json

Step 3. Align to the reference genome (tumor tissue)

Same procedure but for the tumor samples.

Script: 3_gatk4_data_processing_tumor.sh

java -jar gatk-workflows/cromwell-47.jar run gatk-workflows/gatk4-data-processing/processing-for-variant-discovery-gatk4.wdl --inputs gatk-workflows/gatk4-data-processing/processing-for-variant-discovery-gatk4.hg38.wgs.tumor.inputs.json

Step 4. Collect the results

In the GATK workflow, results are generated under the cromwell-executions directory. Each workflow run is stored in a separate directory identified by a unique UUID (e.g., 2dbed345-2ae8-4eb2-837d-28c98825685e).

The alignment step outputs this UUID, and any scripts used to collect or process results should be updated accordingly to reference the correct workflow directory.

Script: 4_collect_results.sh

# uid=2dbed345-2ae8-4eb2-837d-28c98825685e
uid=f01518a9-caba-4f77-ad7e-189b179b05cc
cp cromwell-executions/PreProcessingForVariantDiscovery_GATK4/$uid/call-GatherBamFiles/execution/normal.hg38.bam.md5 intermediates
cp cromwell-executions/PreProcessingForVariantDiscovery_GATK4/$uid/call-GatherBamFiles/execution/normal.hg38.bam intermediates
cp cromwell-executions/PreProcessingForVariantDiscovery_GATK4/$uid/call-GatherBamFiles/execution/normal.hg38.bai intermediates

# uid=e05ce202-2903-47fa-8d0e-883bff8877d3
uid=28a5b370-cf8f-44e7-af8f-391cdcd4b326
cp cromwell-executions/PreProcessingForVariantDiscovery_GATK4/$uid/call-GatherBamFiles/execution/tumor.hg38.bam.md5 intermediates
cp cromwell-executions/PreProcessingForVariantDiscovery_GATK4/$uid/call-GatherBamFiles/execution/tumor.hg38.bam intermediates
cp cromwell-executions/PreProcessingForVariantDiscovery_GATK4/$uid/call-GatherBamFiles/execution/tumor.hg38.bai intermediates

Step 5. Variant calling

Variants can be identified by comparing two BAM files representing the control (normal) and case (tumor) samples. Mutect2 is used to detect somatic variants specific to cancer. To accurately identify tumor-specific variants, the sample name corresponding to the normal sample must be provided.

Script: 5_variant_call_pair.sh

mkdir outputs
~/local/gatk-4.3.0.0/gatk Mutect2 \
-R /home/shared/GATK/Homo_sapiens_assembly38.fasta \
-I ./intermediates/tumor.hg38.bam \
-I ./intermediates/normal.hg38.bam \
-normal n1 \
-O ./outputs/somatic_pair.vcf.gz

6. Function annotation (Funcotator)

Finally, the functions including genes are annotated to each variant using Funcotator.

Script: 6_funcotator.sh

~/local/gatk-4.3.0.0/gatk Funcotator \
  --variant ./outputs/somatic_pair.vcf.gz \
  --reference /home/shared/GATK/Homo_sapiens_assembly38.fasta \
  --ref-version hg38 \
  --data-sources-path /home/shared/funcotator_dataSources.v1.7.20200521s \
  -O tumor_somatic_variants.funcotated.vcf \
  --output-file-format VCF

Inspect the final VCF file to identify tumor-specific variants. You may need to use dedicated software or tools for VCF analysis.