SliverWorkspace Documentation Center

English

简体中文

English

简体中文

Appearance

GATK Germline SNP INDEL V2.0 ​

This article is an upgraded version of Gatk Germline snps-indels Pipeline for Analyzing Genetic Diseases (Deafness), aiming to provide an out-of-the-box analysis pipeline with deployment and migration as simple as possible. ​

Updates are as follows: ​

  • The human reference genome and other referenced database files have been upgraded from GRCh37 (hg19) to GRCh38 (hg38)

  • The annotation software has been switched from annovar to Ensemble VEP (108.2); Annovar requires a commercial license, while VEP is under the Apache 2.0 license and can be used freely

  • The software used by the pipeline has been changed from pre-installed to docker+conda installation on first use; the initial run initializes the environment and downloads the necessary files, making migration easier

To allow the pipeline to run directly without deployment, a docker container is used here to ensure a consistent runtime environment: see the article Building a Bioinformatics Base Environment Image Based on Docker. The solution adopted here is a docker+conda environment with an SSH service, with the entire pipeline running in a single general-purpose container. ​

The code in this article is long and may be somewhat complex. If you want to run it directly, you can download the workflow file, import it into the sliverworkspace graphical bioinformatics platform and run it directly, then obtain and view the graphical analysis report ​

The related code and workflow files can be accessed at the author's GitHub project page

Workflow import operation

Workflow import operation

Pipeline overview diagram ​

Overview

Data preparation: ​

SRR9993255_R1.fastq.gz SRR9993255_R2.fastq.gz

Test data download

SRX6732499: Targeted NGS of human: blood sample

1 ILLUMINA (Illumina HiSeq X Ten) run: 679,472 spots, 203.8M bases, 80.8Mb downloads

If the downloaded data is in SRA format, you can use the official NCBI tool sra-toolkit to split it into fastq.gz files

fastq-dump SRR9993255 --split-3 --gzip

Obtaining SRR9993255_R1.fastq.gz SRR9993255_R2.fastq.gz

Environment setup: ​

This article uses docker + conda (mamba) as the base analysis environment, ​

Image acquisition: for the installation of docker/docker-compose and image building, see: ​

Building a Bioinformatics Base Environment Image Based on Docker, ​

bash
# Pull the docker image
docker     pull     doujiangbaozi/sliverworkspace:latest

# List docker images
docker     images

Base environment configuration, the docker-compose.yml config file, which can be modified as needed ​

yaml
version: "3"
services:
  SarsCov2:
    image: doujiangbaozi/sliverworkspace:latest
    container_name: Germline
    volumes:
      - /media/sliver/Data/data:/opt/data:rw                               #Mount the raw data, place it under the SC2 directory
      - /media/sliver/Manufacture/GL2/envs:/root/mambaforge-pypy3/envs:rw  #Mount the envs conda environment directory
      - /media/sliver/Manufacture/GL2/config:/opt/config:rw                #Mount the config conda config file directory
      - /media/sliver/Manufacture/GL2/ref:/opt/ref:rw                      #Mount the reference directory
      - /media/sliver/Manufacture/GL2/result:/opt/result:rw                #Mount the intermediate files and output results directory
    ports:
      - "9024:9024"                                                        #ssh connection port, can be modified as needed
    environment:
      - TZ=Asia/Shanghai                                                   #Set the time zone
      - PS=20191124                                                        #Change the default ssh password
      - PT=9024                                                            #Change the default ssh connection port

Running the base environment

bash
# Run in the directory containing docker-compose.yml
docker-compose up -d

# Or
docker-compose up -d -f /path/docker-compose.yaml

# Check whether docker is running normally, run in the docker-compose.yaml directory
docker-compose ps

# Or
docker ps

Using the docker container is similar to logging into a remote server

bash
# Log in to docker via the ssh service; it can be deployed locally or remotely
ssh root@192.168.6.6 -p9024

# If you see the following prompt, you have logged in successfully
(base) root@SliverWorkstation:~#

Software versions ​

Software nameSoftware versionRemarks
fastqc0.12.1
multiqc1.21
bwa0.7.17
sambamba1.0.1The static build runs about twice as fast as the mamba-installed version
samtools1.16.1
bedtools2.30.0
fastp0.23.2
gatk4.5.0.0
ensembl-vep108.2VEP software version
version_vep_cache108VEP cache version

Analysis results ​

File nameRemarks
report.pdfReport printed as a PDF file based on the report designer
SRR9993255/qc/multiqc_report.htmlmultiqc report file
SRR9993255/SRR9993255.result.qc.tsvData QC results
SRR9993255/SRR9993255.result.variants.tsvPathogenic variant results
SRR9993255/vcf/SRR9993255_filtered.vcfVCF file after hard filtering
SRR9993255/vcf/SRR9993255_filtered_vep.tsvFile annotated with VEP
SRR9993255/qc/SRR9993255_depth.tsvSequencing depth statistics obtained with samtools depth, used for plotting to display the overall sequencing depth

Analysis workflow ​

Variable setup: ​

Variable nameValueRemarks
snSRR9993255Sample ID
pnGM02Pipeline code
project_bedhg38.exon.bedWhole-exome BED file for the hg38 reference genome
version_referencehg38Reference genome version hg
version_fastqc0.12.1
version_multiqc1.21
version_fastp0.23.4
version_bwa0.7.17
version_samtools1.20
version_bedtoos2.31.1
version_sambamba1.0.1
version_gatk4.5.0.0
version_vep108.2ensembl-vep software version
version_vep_cache108ensembl-vep database version
bash
#Variable initialization and assignment
sn=SRR9993255                           #Sample ID
pn=GM02                                 #Pipeline project code
data=/opt/data                          #Data input directory
result=/opt/result                      #Data output and intermediate files directory
envs=/root/mambaforge-pypy3/envs        #Directory of environments installed by conda
threads=32                              #Number of available threads
whitelist=/opt/ref/whitelist.csv        #File of deafness-related genes and variants
version_reference=hg38                  #Reference genome version
project_bed=hg38.exon.bed               #BED file used by the pipeline; if absent, one will be computed
version_fastqc=0.12.1                   #fastqc version
version_multiqc=1.21                    #multiqc version
version_fastp=0.23.4                    #fastp version
version_bwa=0.7.18                      #bwa version
version_samtools=1.20                   #samtools version
version_bedtools=2.31.1                 #bedtools version
version_sambamba=1.0.1                  #sambamba version
version_gatk=4.5.0.0                    #gatk version
version_vep=108.2                       #VEP software version
version_vep_cache=108                   #vep-cache version

1. Filtering with fastp using default parameters to check raw data quality, clean data: ​

bash
#conda checks whether the environment exists; on first run, if it does not exist, the environment is created and software installed
if [ ! -d "${envs}/${pn}.fastp" ]; then
    echo "Creating the environment ${pn}.fastp"
    mamba create -n ${pn}.fastp -y fastp=${version_fastp} fastqc=${version_fastqc} multiqc=${version_multiqc}
fi

mamba activate ${pn}.fastp

mkdir -p ${result}/${sn}/clean
mkdir -p ${result}/${sn}/qc

fastqc -t ${threads}\
	${data}/${pn}/${sn}_R1.fastq.gz \
	${data}/${pn}/${sn}_R2.fastq.gz \
	-o ${result}/${sn}/qc &

#Using default parameters
fastp \
	-w ${threads} \
	-i ${data}/${pn}/${sn}_R1.fastq.gz \
	-I ${data}/${pn}/${sn}_R2.fastq.gz \
	-o ${result}/${sn}/clean/${sn}_fastp_R1.fastq.gz \
	-O ${result}/${sn}/clean/${sn}_fastp_R2.fastq.gz \
	-j ${result}/${sn}/qc/${sn}_fastp.json \
	-h ${result}/${sn}/qc/${sn}_fastp.html &

wait

fastqc -t ${threads}\
	${result}/${sn}/clean/${sn}_fastp_R1.fastq.gz \
	${result}/${sn}/clean/${sn}_fastp_R2.fastq.gz \
	-o ${result}/${sn}/qc

mamba deactivate

2. Alignment, sorting, duplicate marking; QC ​

bash
#conda checks whether the environment exists; on first run, if it does not exist, the environment is created and software installed
if [ ! -d "${envs}/${pn}.align" ]; then
  mamba create -n ${pn}.align -y bwa=${version_bwa} samtools=${version_samtools} bedtools=${version_bedtools}
fi
mamba activate ${pn}.align

#Download the sambamba static build from GitHub; it is more than twice as fast as the mamba-installed version.
if [ ! -f "${envs}/${pn}.align/bin/sambamba" ]; then
	aria2c https://github.com/biod/sambamba/releases/download/v${version_sambamba}/sambamba-${version_sambamba}-linux-amd64-static.gz -d ${envs}/${pn}.align/bin
	gzip -cdf ${envs}/${pn}.align/bin/sambamba-${version_sambamba}-linux-amd64-static.gz  >  ${envs}/${pn}.align/bin/sambamba 
	chmod a+x ${envs}/${pn}.align/bin/sambamba
fi
#Install bamdst from GitHub source
if [ ! -f "${envs}/bamdst/bamdst" ]; then
	apt-get update && apt-get install -y git make gcc zlib1g-dev
	git clone https://github.com/shiquan/bamdst.git "${envs}/bamdst"
	cd ${envs}/bamdst
	make
	cd ${result}
fi

if [ "${version_reference}" == "hg19" ]; then
    echo	"USE reference Version : ${version_reference}"
    mkdir	-p /opt/ref/hg19
    #If the reference genome sequence is not detected, download the sequence and create the index with bwa
    if [ ! -f "/opt/ref/hg19/hg19.fasta" ]; then
    	aria2c -x 16 -s 32 ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg19/ucsc.hg19.fasta.gz -d /opt/ref/hg19 -o hg19.fasta.gz
    	cd /opt/ref/hg19 && gzip -d /opt/ref/hg19/hg19.fasta.gz
	else
        if [ -f "/opt/ref/hg19/ucsc.hg19.fasta.gz.aria2" ]; then
            aria2c -x 16 -s 32 ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg19/ucsc.hg19.fasta.gz -c -d /opt/ref/hg19 -o hg19.fasta.gz
            cd /opt/ref/hg19 && gzip -d /opt/ref/hg19/hg19.fasta.gz
        fi
	fi
	#Check whether index files exist; if not, create the index
    if  [ ! -f /opt/ref/hg19/hg19.fasta.amb ] ||
        [ ! -f /opt/ref/hg19/hg19.fasta.ann ] ||
        [ ! -f /opt/ref/hg19/hg19.fasta.bwt ] ||
        [ ! -f /opt/ref/hg19/hg19.fasta.pac ] ||
        [ ! -f /opt/ref/hg19/hg19.fasta.sa ]; then
        bwa index /opt/ref/hg19/hg19.fasta
    fi
    #Check whether the samtools index exists; if not, create the reference genome index with samtools
    if [ ! -f "/opt/ref/hg19/hg19.fasta.fai" ]; then
    	samtools faidx /opt/ref/hg19/hg19.fasta
    fi
elif [ "${version_reference}" == "hg38" ]; then
    echo	"USE reference Version : ${version_reference}"
    mkdir	-p /opt/ref/hg38
    if [ ! -f "/opt/ref/hg38/hg38.fasta" ]; then
        aria2c ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg38/Homo_sapiens_assembly38.fasta.gz -d /opt/ref/hg38 -o hg38.fasta.gz	
        cd /opt/ref/hg38 && gzip -d /opt/ref/hg38/hg38.fasta.gz
    elif [ -f "/opt/ref/hg38/hg38.fasta.aria2" ]; then
        aria2c ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg38/Homo_sapiens_assembly38.fasta.gz -c -d /opt/ref/hg38 -o hg38.fasta.gz
        cd /opt/ref/hg38 && gzip -d /opt/ref/hg38/hg38.fasta.gz
    fi
    #Check whether index files exist; if not, create the index
    if  [ ! -f /opt/ref/hg38/hg38.fasta.amb ] ||
        [ ! -f /opt/ref/hg38/hg38.fasta.ann ] ||
        [ ! -f /opt/ref/hg38/hg38.fasta.bwt ] ||
        [ ! -f /opt/ref/hg38/hg38.fasta.pac ] ||
        [ ! -f /opt/ref/hg38/hg38.fasta.sa ]; then
        bwa index /opt/ref/hg38/hg38.fasta
    fi
    #Check whether the samtools index exists; if not, create the reference genome index with samtools
    if [ ! -f "/opt/ref/hg38/hg38.fasta.fai" ]; then
        samtools faidx /opt/ref/hg38/hg38.fasta
    fi
fi


cd       ${result}
mkdir -p ${result}/${sn}/align

bwa mem \
    -t ${threads} -M \
    -R "@RG\\tID:${sn}\\tLB:${sn}\\tPL:Illumina\\tPU:NextSeq550\\tSM:${sn}" \
    /opt/ref/${version_reference}/${version_reference}.fasta  \
    ${result}/${sn}/clean/${sn}_fastp_R1.fastq.gz \
    ${result}/${sn}/clean/${sn}_fastp_R2.fastq.gz \
    | sambamba view -S -f bam -l 0 /dev/stdin \
    | sambamba sort -t ${threads} -m 2G \
    	--tmpdir=${result}/${sn} \
        -o ${result}/${sn}/align/${sn}_sorted.bam /dev/stdin

#Temporarily increase the number of open files, otherwise errors easily occur
ulimit -n 10240

sambamba  markdup \
	--tmpdir ${result}/${sn} \
	-t ${threads} \
    ${result}/${sn}/align/${sn}_sorted.bam \
    ${result}/${sn}/align/${sn}_marked.bam 
mv  ${result}/${sn}/align/${sn}_marked.bam.bai ${result}/${sn}/align/${sn}_marked.bai
rm  -f ${result}/${sn}/align/${sn}_sorted.bam


if [ "${project_bed}" == "None" ]; then
    #The data file does not come with a BED file, so compute one with bedtools (reverse complement)
    bedtools genomecov 	-ibam ${result}/${sn}/align/${sn}_marked.bam -bg > ${result}/${sn}/align/${sn}_bedtools.bed
    bedtools merge     	-i    ${result}/${sn}/align/${sn}_bedtools.bed   > ${result}/${sn}/align/${sn}_merged.bed
	samtools flagstat --threads ${threads} \
		${result}/${sn}/align/${sn}_marked.bam > ${result}/${sn}/qc/${sn}_marked.flagstat &
	${envs}/bamdst/bamdst \
		-p ${result}/${sn}/align/${sn}_merged.bed \
		-o ${result}/${sn}/qc \
		--cutoffdepth 20 \
		${result}/${sn}/align/${sn}_marked.bam &
else
	if [ ! -f "/opt/ref/${version_reference}/${project_bed}" ]; then
		mkdir  -p /opt/ref/${version_reference}
		aria2c https://raw.githubusercontent.com/doujiangbaozi/sliverworkspace-util/main/germline/projects/${project_bed} -d /opt/ref/${version_reference}
	fi
	samtools flagstat --threads ${threads} \
		${result}/${sn}/align/${sn}_marked.bam > ${result}/${sn}/qc/${sn}_marked.flagstat &
	${envs}/bamdst/bamdst \
		-p /opt/ref/${version_reference}/${project_bed} \
		-o ${result}/${sn}/qc \
		--cutoffdepth 20 \
		${result}/${sn}/align/${sn}_marked.bam &
fi
wait

mamba deactivate

3. GATK: generate the base quality score recalibration table ​

bash
#conda checks whether the environment exists; on first run, if it does not exist, the environment is created and software installed
if [ ! -f "${envs}/${pn}.gatk/bin/gatk" ]; then
mamba create -n ${pn}.gatk -y gatk4=${version_gatk} \
    r-base=3.6.2 \
    r-data.table=1.12.8 \
    r-dplyr=0.8.5 \
    r-getopt=1.20.3 \
    r-ggplot2=3.3.0 \
    r-gplots=3.0.3 \
    r-gsalib=2.1 \
    r-optparse=1.6.4 \
    r-backports=1.1.10 \
    tabix \
    pandas
fi

mamba activate ${pn}.gatk

if [ "${version_reference}" == hg19 ]; then
    echo	"USE reference Version : ${version_reference}"
    #Here is a huge pitfall: for the hg19 reference files provided in the bundle directory of the Broad Institute FTP site, running with the default format throws errors saying there is no index, and creating the index with GATK still fails; in fact, the gz files need to be recompressed with bgzip and indexed with tabix
    if [ ! -f "/opt/ref/hg19/dbsnp_138.hg19.vcf.gz" ]; then
    	aria2c -x 16 -s 32 ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg19/dbsnp_138.hg19.vcf.gz -d /opt/ref/hg19
    else 
        if [ -f "/opt/ref/hg19/dbsnp_138.hg19.vcf.gz.aria2" ]; then
            echo 'download continue...'
            aria2c -x 16 -s 32 ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg19/dbsnp_138.hg19.vcf.gz -c -d /opt/ref/hg19
        fi
    fi

    if [ ! -f "/opt/ref/hg19/dbsnp_138.hg19.vcf.gz.tbi" ]; then
        gzip -k -f -d /opt/ref/hg19/dbsnp_138.hg19.vcf.gz > /opt/ref/hg19/dbsnp_138.hg19.vcf
        bgzip -f --threads ${threads} /opt/ref/hg19/dbsnp_138.hg19.vcf > /opt/ref/hg19/dbsnp_138.hg19.vcf.gz
        tabix -f /opt/ref/hg19/dbsnp_138.hg19.vcf.gz
    fi

    if [ ! -f "/opt/ref/hg19/Mills_and_1000G_gold_standard.indels.hg19.sites.vcf.gz" ]; then
        aria2c -x 16 -s 32 ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg19/Mills_and_1000G_gold_standard.indels.hg19.sites.vcf.gz -d /opt/ref/hg19
    else 
        if [ -f "/opt/ref/hg19/Mills_and_1000G_gold_standard.indels.hg19.sites.vcf.gz.aria2" ]; then
            echo 'download continue...'
            aria2c -x 16 -s 32 ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg19/Mills_and_1000G_gold_standard.indels.hg19.sites.vcf.gz -c -d /opt/ref/hg19
        fi
    fi

    if [ ! -f "/opt/ref/hg19/Mills_and_1000G_gold_standard.indels.hg19.sites.vcf.gz.tbi" ]; then
        gzip -k -f -d /opt/ref/hg19/Mills_and_1000G_gold_standard.indels.hg19.sites.vcf.gz > /opt/ref/hg19/Mills_and_1000G_gold_standard.indels.hg19.sites.vcf
        bgzip -f --threads ${threads} /opt/ref/hg19/Mills_and_1000G_gold_standard.indels.hg19.sites.vcf > /opt/ref/hg19/Mills_and_1000G_gold_standard.indels.hg19.sites.vcf.gz
        tabix -f /opt/ref/hg19/Mills_and_1000G_gold_standard.indels.hg19.sites.vcf.gz
    fi

    if [ ! -f "/opt/ref/hg19/1000G_phase1.snps.high_confidence.hg19.sites.vcf.gz" ]; then
        aria2c -x 16 -s 32 ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg19/1000G_phase1.snps.high_confidence.hg19.sites.vcf.gz -d /opt/ref/hg19
    else 
        if [ -f "/opt/ref/hg19/1000G_phase1.snps.high_confidence.hg19.sites.vcf.gz.aria2" ]; then
            echo 'download continue...'
            aria2c -x 16 -s 32 ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg19/1000G_phase1.snps.high_confidence.hg19.sites.vcf.gz -c -d /opt/ref/hg19
        fi
    fi

    if [ ! -f "/opt/ref/hg19/1000G_phase1.snps.high_confidence.hg19.sites.vcf.gz.tbi" ]; then
        gzip -k -f -d /opt/ref/hg19/1000G_phase1.snps.high_confidence.hg19.sites.vcf.gz > /opt/ref/hg19/1000G_phase1.snps.high_confidence.hg19.sites.vcf
        bgzip -f --threads ${threads} /opt/ref/hg19/1000G_phase1.snps.high_confidence.hg19.sites.vcf > /opt/ref/hg19/1000G_phase1.snps.high_confidence.hg19.sites.vcf.gz
        tabix -f /opt/ref/hg19/1000G_phase1.snps.high_confidence.hg19.sites.vcf.gz
    fi

    if [ ! -f "/opt/ref/hg19/1000G_phase1.indels.hg19.sites.vcf.gz" ]; then
        aria2c -x 16 -s 32 ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg19/1000G_phase1.indels.hg19.sites.vcf.gz -d /opt/ref/
    else 
        if [ -f "/opt/ref/hg19/1000G_phase1.indels.hg19.sites.vcf.gz.aria2" ]; then
            echo 'download continue...'
            aria2c -x 16 -s 32 ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg19/1000G_phase1.indels.hg19.sites.vcf.gz -c -d /opt/ref/hg19
        fi
    fi

    if [ ! -f "/opt/ref/hg19/1000G_phase1.indels.hg19.sites.vcf.gz.tbi" ]; then
        gzip -k -f -d /opt/ref/hg19/1000G_phase1.indels.hg19.sites.vcf.gz > /opt/ref/hg19/1000G_phase1.indels.hg19.sites.vcf
        bgzip -f --threads ${threads} /opt/ref/hg19/1000G_phase1.indels.hg19.sites.vcf > /opt/ref/hg19/1000G_phase1.indels.hg19.sites.vcf.gz
        tabix -f /opt/ref/hg19/1000G_phase1.indels.hg19.sites.vcf.gz
    fi
    #Create the dict dictionary file for the hg19 reference sequence
    if [ ! -f "/opt/ref/hg19/hg19.dict" ]; then
        gatk CreateSequenceDictionary -R /opt/ref/hg19/hg19.fasta -O /opt/ref/hg19/hg19.dict
    fi

    if [ "${project_bed}" == "None" ]; then
        if [ ! -f "${result}/${sn}/align/${sn}_merged.bed.interval_list" ]; then
            gatk BedToIntervalList \
                -I	${result}/${sn}/align/${sn}_merged.bed \
                -SD	/opt/ref/${version_reference}/${version_reference}.dict \
                -O	${result}/${sn}/align/${sn}_merged.bed.interval_list
        fi
        mkdir -p ${result}/${sn}/recal
        gatk BaseRecalibrator \
            --known-sites /opt/ref/hg19/dbsnp_138.hg19.vcf.gz \
            --known-sites /opt/ref/hg19/Mills_and_1000G_gold_standard.indels.hg19.sites.vcf.gz \
            --known-sites /opt/ref/hg19/1000G_phase1.indels.hg19.sites.vcf.gz \
            -L ${result}/${sn}/align/${sn}_merged.bed.interval_list \
            -R /opt/ref/${version_reference}/${version_reference}.fasta \
            -I ${result}/${sn}/align/${sn}_marked.bam \
            -O ${result}/${sn}/align/${sn}_recal.table
    else
        if [ ! -f "/opt/ref/${version_reference}/${project_bed}" ]; then
            #Create the interval list file based on the downloaded Illumina_pt2.bed file, converting coordinates from 0-based to 1-based
            #sed 's/chr//; s/\t/ /g' /opt/ref/hg19/Illumina_pt2.bed > /opt/ref/hg19/Illumina_pt2.processed.bed
            mkdir  -p /opt/ref/${version_reference}
            aria2c https://raw.githubusercontent.com/doujiangbaozi/sliverworkspace-util/main/germline/projects/${project_bed} -d /opt/ref/${version_reference}
        fi
        if [ ! -f "/opt/ref/${version_reference}/${project_bed}.interval_list" ]; then
            gatk BedToIntervalList \
                -I	/opt/ref/${version_reference}/${project_bed} \
                -SD	/opt/ref/${version_reference}/${version_reference}.dict \
                -O	/opt/ref/${version_reference}/${project_bed}.interval_list
        fi
        mkdir -p ${result}/${sn}/recal
        gatk BaseRecalibrator \
            --known-sites /opt/ref/hg19/dbsnp_138.hg19.vcf.gz \
            --known-sites /opt/ref/hg19/Mills_and_1000G_gold_standard.indels.hg19.sites.vcf.gz \
            --known-sites /opt/ref/hg19/1000G_phase1.indels.hg19.sites.vcf.gz \
            -L /opt/ref/${version_reference}/${project_bed}.interval_list \
            -R /opt/ref/${version_reference}/${version_reference}.fasta \
            -I ${result}/${sn}/align/${sn}_marked.bam \
            -O ${result}/${sn}/align/${sn}_recal.table
    fi

elif [ "${version_reference}" == "hg38" ]; then
	echo	"USE reference Version : ${version_reference}"
    if [ ! -f "/opt/ref/hg38/dbsnp_146.hg38.vcf.gz" ]; then
    	aria2c -x 16 -s 32 ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg38/dbsnp_146.hg38.vcf.gz -d /opt/ref/hg38
    else 
    	if [ -f "/opt/ref/hg38/dbsnp_146.hg38.vcf.gz.aria2" ]; then
    		echo 'download continue...'
    		aria2c -x 16 -s 32 ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg38/dbsnp_146.hg38.vcf.gz -c -d /opt/ref/hg38
    	fi
	fi
    if [ ! -f "/opt/ref/hg38/dbsnp_146.hg38.vcf.gz.tbi" ]; then
        aria2c -x 16 -s 32 ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg38/dbsnp_146.hg38.vcf.gz.tbi -d /opt/ref/hg38
    else 
        if [ -f "/opt/ref/hg38/dbsnp_146.hg38.vcf.gz.tbi.aria2" ]; then
            echo 'download continue...'
            aria2c -x 16 -s 32 ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg38/dbsnp_146.hg38.vcf.gz.tbi -c -d /opt/ref/hg38
        fi
    fi

    if [ ! -f "/opt/ref/hg38/Mills_and_1000G_gold_standard.indels.hg38.vcf.gz" ]; then
        aria2c -x 16 -s 32 ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg38/Mills_and_1000G_gold_standard.indels.hg38.vcf.gz -d /opt/ref/hg38
    else 
        if [ -f "/opt/ref/hg38/Mills_and_1000G_gold_standard.indels.hg38.vcf.gz.aria2" ]; then
            echo 'download continue...'
            aria2c -x 16 -s 32 ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg38/Mills_and_1000G_gold_standard.indels.hg38.vcf.gz -c -d /opt/ref/hg38
        fi
    fi
    if [ ! -f "/opt/ref/hg38/Mills_and_1000G_gold_standard.indels.hg38.vcf.gz.tbi" ]; then
        aria2c -x 16 -s 32 ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg38/Mills_and_1000G_gold_standard.indels.hg38.vcf.gz.tbi -d /opt/ref/hg38
    else 
        if [ -f "/opt/ref/hg38/Mills_and_1000G_gold_standard.indels.hg38.vcf.gz.tbi.aria2" ]; then
            echo 'download continue...'
            aria2c -x 16 -s 32 ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg38/Mills_and_1000G_gold_standard.indels.hg38.vcf.gz.tbi -c -d /opt/ref/hg38
        fi
    fi

    if [ ! -f "/opt/ref/hg38/1000G_phase1.snps.high_confidence.hg38.vcf.gz" ]; then
        aria2c -x 16 -s 32 ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg38/1000G_phase1.snps.high_confidence.hg38.vcf.gz -d /opt/ref/hg38
    else 
        if [ -f "/opt/ref/hg38/1000G_phase1.snps.high_confidence.hg38.vcf.gz.aria2" ]; then
            echo 'download continue...'
            aria2c -x 16 -s 32 ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg38/1000G_phase1.snps.high_confidence.hg38.vcf.gz -c -d /opt/ref/hg38
        fi
    fi
    if [ ! -f "/opt/ref/hg38/1000G_phase1.snps.high_confidence.hg38.vcf.gz.tbi" ]; then
        aria2c -x 16 -s 32 ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg38/1000G_phase1.snps.high_confidence.hg38.vcf.gz.tbi -d /opt/ref/hg38
    else 
        if [ -f "/opt/ref/hg38/1000G_phase1.snps.high_confidence.hg38.vcf.gz.tbi.aria2" ]; then
            echo 'download continue...'
            aria2c -x 16 -s 32 ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg38/1000G_phase1.snps.high_confidence.hg38.vcf.gz.tbi -c -d /opt/ref/hg38
        fi
    fi
    #Create the dict dictionary file for the hg38 reference sequence
    if [ ! -f "/opt/ref/hg38/hg38.dict" ]; then
        gatk CreateSequenceDictionary -R /opt/ref/hg38/hg38.fasta -O /opt/ref/hg38/hg38.dict
    fi

    if [ "${project_bed}" == "None" ]; then
        if [ ! -f "${result}/${sn}/align/${sn}_merged.bed.interval_list" ]; then
            gatk BedToIntervalList \
                -I	${result}/${sn}/align/${sn}_merged.bed \
                -SD	/opt/ref/${version_reference}/${version_reference}.dict \
                -O	${result}/${sn}/align/${sn}_merged.bed.interval_list
        fi
        mkdir -p ${result}/${sn}/recal
        gatk BaseRecalibrator \
            --known-sites /opt/ref/hg38/dbsnp_146.hg38.vcf.gz \
            --known-sites /opt/ref/hg38/Mills_and_1000G_gold_standard.indels.hg38.vcf.gz \
            --known-sites /opt/ref/hg38/1000G_phase1.snps.high_confidence.hg38.vcf.gz \
            -L ${result}/${sn}/align/${sn}_merged.bed.interval_list \
            -R /opt/ref/${version_reference}/${version_reference}.fasta \
            -I ${result}/${sn}/align/${sn}_marked.bam \
            -O ${result}/${sn}/align/${sn}_recal.table
    else
        if [ ! -f "/opt/ref/${version_reference}/${project_bed}.interval_list" ]; then
            gatk BedToIntervalList \
                -I	/opt/ref/${version_reference}/${project_bed} \
                -SD	/opt/ref/${version_reference}/${version_reference}.dict \
                -O	/opt/ref/${version_reference}/${project_bed}.interval_list
        fi
        mkdir -p ${result}/${sn}/recal
        gatk BaseRecalibrator \
            --known-sites /opt/ref/hg38/dbsnp_146.hg38.vcf.gz \
            --known-sites /opt/ref/hg38/Mills_and_1000G_gold_standard.indels.hg38.vcf.gz \
            --known-sites /opt/ref/hg38/1000G_phase1.snps.high_confidence.hg38.vcf.gz \
            -L /opt/ref/${version_reference}/${project_bed}.interval_list \
            -R /opt/ref/${version_reference}/${version_reference}.fasta \
            -I ${result}/${sn}/align/${sn}_marked.bam \
            -O ${result}/${sn}/align/${sn}_recal.table
    fi
fi

mamba deactivate

​

4. GATK: apply base quality score recalibration and obtain insert size statistics ​

bash
mamba activate ${pn}.gatk
    
if [ "${project_bed}" == "None" ]; then
    gatk SplitIntervals \
    	-L ${result}/${sn}/align/${sn}_merged.bed.interval_list \
    	-R /opt/ref/${version_reference}/${version_reference}.fasta \
    	-O ${result}/${sn}/interval \
    	--scatter-count ${scatter} &
    
    gatk ApplyBQSR \
    	--bqsr-recal-file ${result}/${sn}/align/${sn}_recal.table \
    	-L ${result}/${sn}/align/${sn}_merged.bed.interval_list \
    	-R /opt/ref/${version_reference}/${version_reference}.fasta \
    	-I ${result}/${sn}/align/${sn}_marked.bam \
    	-O ${result}/${sn}/align/${sn}_bqsr.bam &
    
    gatk CollectInsertSizeMetrics \
    	 -I ${result}/${sn}/align/${sn}_marked.bam \
    	 -O ${result}/${sn}/qc/${sn}_insertsize_metrics.txt \
    	 -H ${result}/${sn}/qc/${sn}_insertsize_histogram.pdf &
else
    rm -f ${result}/${sn}/interval/*.interval_list
    gatk SplitIntervals \
    	-L /opt/ref/${version_reference}/${project_bed}.interval_list \
    	-R /opt/ref/${version_reference}/${version_reference}.fasta \
    	-O ${result}/${sn}/interval \
    	--scatter-count ${scatter} &
    
    gatk ApplyBQSR \
    	--bqsr-recal-file ${result}/${sn}/align/${sn}_recal.table \
    	-L /opt/ref/${version_reference}/${project_bed}.interval_list \
    	-R /opt/ref/${version_reference}/${version_reference}.fasta \
    	-I ${result}/${sn}/align/${sn}_marked.bam \
    	-O ${result}/${sn}/align/${sn}_bqsr.bam &
    
    gatk CollectInsertSizeMetrics \
    	-I ${result}/${sn}/align/${sn}_marked.bam \
    	-O ${result}/${sn}/qc/${sn}_insertsize_metrics.txt \
    	-H ${result}/${sn}/qc/${sn}_insertsize_histogram.pdf &
fi
    
wait
    
mamba deactivate

​

5. GATK HaplotypeCaller to call SNP/indel variants: ​

bash
mamba activate ${pn}.gatk

mkdir -p ${result}/${sn}/vcf

rm -f ${result}/${sn}/vcf/vcf-file.list
touch ${result}/${sn}/vcf/vcf-file.list

if [ "${version_reference}" == hg19 ]; then
	for i in `ls ${result}/${sn}/interval/*.interval_list`;
	do
		rm -f ${i%.*}_bqsr.vcf.gz
		gatk HaplotypeCaller \
            -R /opt/ref/${version_reference}/${version_reference}.fasta \
            -L $i \
            -I ${result}/${sn}/align/${sn}_bqsr.bam \
            -D /opt/ref/${version_reference}/dbsnp_138.hg19.vcf.gz \
            -O ${i%.*}.vcf
		echo ${i%.*}.vcf >> ${result}/${sn}/vcf/vcf-file.list
	done
elif [ "${version_reference}" == "hg38" ]; then
	for i in `ls ${result}/${sn}/interval/*.interval_list`;
	do
		rm -f ${i%.*}_bqsr.vcf.gz
		gatk HaplotypeCaller \
			-R /opt/ref/${version_reference}/${version_reference}.fasta \
			-L $i \
			-I ${result}/${sn}/align/${sn}_bqsr.bam \
			-D /opt/ref/${version_reference}/dbsnp_146.hg38.vcf.gz \
			-O ${i%.*}.vcf
		echo ${i%.*}.vcf >> ${result}/${sn}/vcf/vcf-file.list
	done
fi

wait

gatk MergeVcfs \
    -I ${result}/${sn}/vcf/vcf-file.list \
    -O ${result}/${sn}/vcf/${sn}_bqsr.vcf
	
mamba deactivate

​

6. A single sample does not provide enough data to run VQSR, so hard filtering is applied directly; see reference values at #https://gatk.broadinstitute.org/hc/en-us/articles/360035532412?id=11097 ​

bash
mamba activate ${pn}.gatk
#https://gatk.broadinstitute.org/hc/en-us/articles/360035532412?id=11097
gatk VariantFiltration \
    -R /opt/ref/${version_reference}/${version_reference}.fasta \
    -V ${result}/${sn}/vcf/${sn}.vcf \
    -O ${result}/${sn}/vcf/${sn}_snp.vcf \
    --filter-name "SNP_QD" \
    --filter-expression "QD < 2.0" \
    --filter-name "SNP_MQ" \
    --filter-expression "MQ <40.0" \
    --filter-name "SNP_FS" \
    --filter-expression "FS > 60.0" \
    --filter-name "SNP_SOR" \
    --filter-expression "SOR > 3.0" \
    --filter-name "SNP_MQRankSum" \
    --filter-expression "MQRankSum < -12.5" \
    --filter-name "SNP_ReadPosRankSum" \
    --filter-expression "ReadPosRankSum < -8.0"

gatk VariantFiltration \
    -R /opt/ref/${version_reference}/${version_reference}.fasta \
    -V ${result}/${sn}/vcf/${sn}_snp.vcf \
    -O ${result}/${sn}/vcf/${sn}_filtered.vcf \
    --filter-name "INDEL_QD" \
    --filter-expression "QD < 2.0" \
    --filter-name "INDEL_FS" \
    --filter-expression "FS > 200.0" \
    --filter-name "INDEL_SOR" \
    --filter-expression "SOR > 10.0" \
    --filter-name "INDEL_ReadPosRankSum" \
    --filter-expression "ReadPosRankSum < -20.0" \
    --filter-name "INDEL_InbreedingCoeff" \
    --filter-expression "InbreedingCoeff < -0.8"
	
mamba deactivate

7. Annotate the variant results with VEP: ​

bash
#conda checks whether the environment exists; on first run, if it does not exist, the environment is created and software installed
if [ ! -d "${envs}/${pn}.vep" ]; then
  echo "Creating the environment ${pn}.vep"
  mamba create -n ${pn}.vep -y ensembl-vep=${version_vep}
fi

mamba activate ${pn}.vep

mkdir -p /opt/result/${sn}/vcf

if [ "${version_reference}" == "hg19" ]; then
	if [ ! -d "/opt/ref/vep-cache/homo_sapiens/${version_vep_cache}_GRCh37" ]; then
        aria2c -x 16 -s 48 https://ftp.ensembl.org/pub/release-${version_vep_cache}/variation/indexed_vep_cache/homo_sapiens_vep_${version_vep_cache}_GRCh37.tar.gz -d /opt/ref/
        tar -zxvf /opt/ref/homo_sapiens_vep_${version_vep_cache}_GRCh37.tar.gz -C /opt/ref/vep-cache/
    elif [ -f "/opt/ref/homo_sapiens_vep_${version_vep_cache}_GRCh37.tar.gz.aria2" ]; then
        aria2c -x 16 -s 48 https://ftp.ensembl.org/pub/release-${version_vep_cache}/variation/indexed_vep_cache/homo_sapiens_vep_${version_vep_cache}_GRCh37.tar.gz -c -d /opt/ref/
        tar -zxvf /opt/ref/homo_sapiens_vep_${version_vep_cache}_GRCh37.tar.gz -C /opt/ref/vep-cache/
    fi
	vep \
	  -i ${result}/${sn}/vcf/${sn}_filtered.vcf  \
	  -o ${result}/${sn}/vcf/${sn}_filtered_vep.tsv \
	  --offline \
	  --cache \
	  --cache_version ${version_vep_cache} \
	  --everything \
	  --dir_cache /opt/ref/vep-cache \
	  --dir_plugins /opt/ref/vep-cache/Plugins \
	  --species homo_sapiens \
	  --assembly GRCh37 \
	  --hgvs \
	  --fasta /opt/ref/${version_reference}/${version_reference}.fasta \
	  --force_overwrite \
	  --format vcf \
	  --tab \
	  --fork ${threads} \
	  --offline
elif [ "${version_reference}" == "hg38" ]; then
	if [ ! -d "/opt/ref/vep-cache/homo_sapiens/${version_vep_cache}_GRCh38" ]; then
		aria2c -x 16 -s 48 https://ftp.ensembl.org/pub/release-${version_vep_cache}/variation/indexed_vep_cache/homo_sapiens_vep_${version_vep_cache}_GRCh38.tar.gz -d /opt/ref/
		tar -zxvf /opt/ref/homo_sapiens_vep_${version_vep_cache}_GRCh38.tar.gz -C /opt/ref/vep-cache/
	elif [ -f "/opt/ref/homo_sapiens_vep_${version_vep_cache}_GRCh38.tar.gz.aria2" ]; then
		aria2c -x 16 -s 48 https://ftp.ensembl.org/pub/release-${version_vep_cache}/variation/indexed_vep_cache/homo_sapiens_vep_${version_vep_cache}_GRCh38.tar.gz -c -d /opt/ref/
		tar -zxvf /opt/ref/homo_sapiens_vep_${version_vep_cache}_GRCh38.tar.gz -C /opt/ref/vep-cache/
	fi
	vep \
	  -i ${result}/${sn}/vcf/${sn}_filtered.vcf  \
	  -o ${result}/${sn}/vcf/${sn}_filtered_vep.tsv \
	  --offline \
	  --cache \
	  --cache_version ${version_vep_cache} \
	  --everything \
	  --dir_cache /opt/ref/vep-cache \
	  --dir_plugins /opt/ref/vep-cache/Plugins \
	  --species homo_sapiens \
	  --assembly GRCh38 \
	  --hgvs \
	  --fasta /opt/ref/${version_reference}/${version_reference}.fasta \
	  --force_overwrite \
	  --format vcf \
	  --tab \
	  --fork ${threads} \
	  --offline
fi

mamba deactivate

8. Write a script to filter the variant VCF file (annotated with VEP) and obtain pathogenic variants with ClinVar annotation ​

​

bash
mamba activate ${pn}.gatk 

if [ ! -f "${envs}/GermlineVepAnnotationUtil.py" ]; then
    aria2c https://raw.githubusercontent.com/doujiangbaozi/sliverworkspace-util/main/somatic/GermlineVepAnnotationUtil.py -d ${envs}/
    chmod a+x ${envs}/GermlineVepAnnotationUtil.py
fi

${envs}/GermlineVepAnnotationUtil.py \
    -v ${result}/${sn}/vcf/${sn}_filtered.vcf \
    -a ${result}/${sn}/vcf/${sn}_filtered_vep.tsv \
    -o ${result}/${sn}/${sn}.result.variants.tsv

mamba deactivate

​

9. Write a script to obtain QC information from the outputs of fastp, bamdst, and GATK CollectInsertSize ​

bash
mamba	activate ${pn}.gatk

if [ ! -f "${envs}/GermlineQcUtil.py" ]; then
    aria2c https://raw.githubusercontent.com/doujiangbaozi/sliverworkspace-util/main/somatic/GermlineQcUtil.py -d ${envs}/
    chmod a+x ${envs}/GermlineQcUtil.py
fi

python ${envs}/GermlineQcUtil.py \
    --out=${result}/${sn}/${sn}.result.qc.tsv \
    --sample-fastp=${result}/${sn}/qc/${sn}_fastp.json  \
    --sample-bamdst=${result}/${sn}/qc/coverage.report \
    --sample-insertsize=${result}/${sn}/qc/${sn}_insertsize_metrics.txt

mamba	deactivate

mamba activate ${pn}.fastp

multiqc ${result}/${sn}/ -f -o ${result}/${sn}/qc

mamba deactivate

  1. Result confirmation: bam file and variant positions in IGV:

  2. Report generated from a Word docx template:

    Output report

© 2019 - SliverWorkspace. All rights reserved.