GATK Somatic SNV+Indel+CNV+SV

Recently I was preparing to develop a report designer for the sliverworkspace graphical bioinformatics platform, which required a fairly complex pipeline as test data, so I dug out my previous [[GATK Somatic SNV+Indel+CNV+SV for Full-Mark External Quality Assessment (Part 2): Performance Optimization]]. After running it once, I found all sorts of problems, so I decided to rework the pipeline into a deployment-free version that initializes the environment on first run, so that it could be launched directly whenever needed. Thus this article came into being.

To describe it in one sentence: an out-of-the-box pipeline that automatically selects the reference genome version based on the sample's version_reference, selects the project BED based on the project_bed file, and automatically initializes the environment, installs the required software, and downloads the reference files and database versions

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

Overall Pipeline Overview

Overview diagram

Docker image pull and deployment configuration

bash

# Pull the docker image
docker     pull     doujiangbaozi/sliverworkspace:1.10

docker-compose.yml configuration file

yaml

version: "3"
services:
  GATK:
    image: doujiangbaozi/sliverworkspace:1.10
    container_name: GATK
    volumes:
      - /home/sliver/Data/data:/opt/data:rw                                #Mount the input data directory
      - /home/sliver/Manufacture/gatk/envs:/root/mambaforge-pypy3/envs     #Mount the envs directory
      - /home/sliver/Manufacture/sliver/ref:/opt/ref:rw                    #Mount the reference directory
      - /home/sliver/Manufacture/gatk/result:/opt/result:rw                #Mount the intermediate files and results directory
    environment:
      - TZ=Asia/Shanghai                                                   #Set the time zone
      - PS=20191124                                                        #Set the ssh password
      - PT=9024                                                            #Set the ssh connection port

Deploying and running the docker image

bash

# Run in the directory at the same level as docker-compose.yml
docker-compose up -d

# Or docker-compose -f with the directory containing docker-compose.yml
docker-compose -f somedir/docker-compose.yml up -d

# You can now connect to docker via ssh and run pipeline commands; the connection port and password are in the relevant fields of the docker-compose.yml config file
ssh root@127.0.0.1 -p9024

Test data

The test data comes from the BED files (downloaded automatically by the pipeline) and test data provided by the 2017 National Health and Family Planning Commission external quality assessment (EQA), renamed to match the pipeline's input. You can also replace them with your own data files. Since the EQA reference genome is still at hg19, this pipeline still uses hg19 (GRCH37). To switch to hg38, set the version_reference variable to hg38 and project_bed to Illumina_pt2_hg38.bed; the pipeline will then use the hg38 (GRCH38) version and the corresponding BED file.

File name (adjusted as needed)	File size	MD5
B1701_R1.fq.gz	4.85G	07d3cdccee41dbb3adf5d2e04ab28e5b
B1701_R2.fq.gz	4.77G	c2aa4a8ab784c77423e821b9f7fb00a7
B1701NC_R1.fq.gz	3.04G	4fc21ad05f9ca8dc93d2749b8369891b
B1701NC_R2.fq.gz	3.11G	bc64784f2591a27ceede1727136888b9

Variable names

bash

# Variable initialization and assignment
  sn=1701                               #Sample ID
  pn=GS03                               #Pipeline code
  data=/opt/data                        #Data input directory
  result=/opt/result                    #Data output and intermediate files directory
  project_bed=Illumina_pt2.bed          #BED file for the hg19 reference genome
  version_reference=hg19                #Human reference genome version: hg19 or hg38
  version_fastqc=0.12.1                 #fastqc version
  version_multiqc=1.21                  #multiqc version
  version_cnvkit=0.9.10                 #cnvkit version
  version_manta=1.6.0                   #manta version
  version_gatk=4.3.0.0                  #gatk version                                 
  version_sambamba=1.0.1                #sambamba version                             
  version_samtools=1.17                 #samtools version                             
  version_bwa=0.7.17                    #bwa version                                  
  version_fastp=0.23.2                  #fastp version                                
  version_vep=108.2                     #VEP software version
  version_vep_cache=108                 #vep-cache database version                               
  envs=/root/mambaforge-pypy3/envs	    #mamba envs directory                           
  threads=32                         	#Maximum number of threads                                   
  memory=32G                        	#Memory usage                                    
  scatter=8                          	#Number of interval list splits for Mutect2 parallel runs          
  event=2                          	    #Value of gatk FilterMutectCalls --max-events-in-region
  snv_tlod=16.00                      	#snv filter parameter: tlod value                        
  snv_vaf=0.01                       	#snv filter parameter: abundance/variant allele frequency                        
  snv_depth=500                        	#snv filter parameter: supporting reads count/sequencing depth            
  cnv_dep=1000                       	#cnv filter parameter: supporting reads count/sequencing depth            
  cnv_min=-0.5                       	#cnv filter parameter: minimum log2 value                        
  cnv_max=0.5                        	#cnv filter parameter: maximum log2 value                       
  sv_score=200                        	#sv filter parameter: score                           

## The above variables can be adjusted according to specific needs

Table:

Variable name	Value	Remarks
sn	1701	Sample ID
pn	GS03	Pipeline code, GATK Somatic version 03
project_bed	Illumina_pt2.bed	BED file for the hg19 reference genome
version_reference	hg19	Human reference genome version: hg19 / hg38
version_fastqc	0.12.1	fastqc version
version_multiqc	1.21	multiqc version
version_cnvkit	0.9.10	cnvkit version
version_manta	1.6.0	manta version
version_gatk	4.3.0.0	gatk version
version_sambamba	1.0.1	sambamba version
version_samtools	1.17	samtools version
version_bwa	0.7.17	bwa version
version_fastp	0.23.2	fastp version
version_vep	108.2	VEP software version
version_vep_cache	108	VEP cache version
envs	/root/mambaforge-pypy3/envs	mamba envs directory
threads	32	Maximum number of threads
memory	32G	Memory usage
scatter	8	Number of interval list splits for Mutect2 parallel runs
event	2	Value of gatk FilterMutectCalls --max-events-in-region
snv_tlod	16.00	snv filter parameter: tlod value
snv_vaf	0.01	snv filter parameter: abundance/variant allele frequency
snv_depth	500	snv filter parameter: supporting reads count/sequencing depth
cnv_dep	1000	cnv filter parameter: supporting reads count/sequencing depth
cnv_min	-0.5	cnv filter parameter: minimum log2 value
cnv_max	0.5	cnv filter parameter: maximum log2 value
sv_score	200	sv filter parameter: score

Detailed Pipeline/workflow steps:

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}/trimmed
mkdir -p ${result}/${sn}/qc

fastqc -t ${threads}\
	${data}/GS03/${sn}_tumor_R1.fq.gz \
	${data}/GS03/${sn}_tumor_R2.fq.gz \
	-o ${result}/${sn}/qc &

fastqc -t ${threads}\
	${data}/GS03/${sn}_normal_R1.fq.gz \
	${data}/GS03/${sn}_normal_R2.fq.gz \
	-o ${result}/${sn}/qc &

fastp -w 16 \
    -i ${data}/GS03/${sn}_tumor_R1.fq.gz  \
    -I ${data}/GS03/${sn}_tumor_R2.fq.gz  \
    -o ${result}/${sn}/trimmed/${sn}_tumor_R1_trimmed.fq.gz \
    -O ${result}/${sn}/trimmed/${sn}_tumor_R2_trimmed.fq.gz \
    -h ${result}/${sn}/trimmed/${sn}_tumor_fastp.html \
    -j ${result}/${sn}/trimmed/${sn}_tumor_fastp.json &

fastp -w 16 \
    -i ${data}/GS03/${sn}_normal_R1.fq.gz  \
    -I ${data}/GS03/${sn}_normal_R2.fq.gz  \
    -o ${result}/${sn}/trimmed/${sn}_normal_R1_trimmed.fq.gz \
    -O ${result}/${sn}/trimmed/${sn}_normal_R2_trimmed.fq.gz \
    -h ${result}/${sn}/trimmed/${sn}_normal_fastp.html \
    -j ${result}/${sn}/trimmed/${sn}_normal_fastp.json &
wait

fastqc -t ${threads}\
	${result}/${sn}/trimmed/${sn}_tumor_R1_trimmed.fq.gz \
	${result}/${sn}/trimmed/${sn}_tumor_R2_trimmed.fq.gz \
	-o ${result}/${sn}/qc &

fastqc -t ${threads}\
	${result}/${sn}/trimmed/${sn}_normal_R1_trimmed.fq.gz \
	${result}/${sn}/trimmed/${sn}_normal_R2_trimmed.fq.gz \
	-o ${result}/${sn}/qc &

wait

mamba	deactivate

2. Align the normal fastq files to the reference genome, sort, mark duplicates to obtain marked.bam

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} 
fi

#Download the sambamba static build from GitHub; it is more than twice as fast as the mamba-installed version. This is quite a strange phenomenon
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

mamba	activate ${pn}.align

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
	
    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 the reference genome sequence is not detected, download the sequence and create the index with bwa
    if [ ! -f "/opt/ref/hg38/hg38.fasta" ]; then
    	aria2c -x 16 -s 32 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
	else
        if [ -f "/opt/ref/hg38/hg38.fasta.gz.aria2" ]; then
            aria2c -x 16 -s 32 ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg38/Homo_sapiens_assembly38.fasta.gz -d /opt/ref/hg38 -d /opt/ref/hg38 -o hg38.fasta.gz
            cd /opt/ref/hg38 && gzip -d /opt/ref/hg38/hg38.fasta.gz
        fi
    fi
    
    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




mkdir -p ${result}/${sn}/aligned
#Pipe the genome alignment output into a BAM file, sorted via the pipe
bwa mem \
    -t ${threads} -M \
    -R "@RG\\tID:${sn}_normal\\tLB:${sn}_normal\\tPL:Illumina\\tPU:Miseq\\tSM:${sn}_normal" \
    /opt/ref/${version_reference}/${version_reference}.fasta  ${result}/${sn}/trimmed/${sn}_normal_R1_trimmed.fq.gz ${result}/${sn}/trimmed/${sn}_normal_R2_trimmed.fq.gz \
    | sambamba view -S -f bam -l 0 /dev/stdin \
    | sambamba sort -t ${threads} -m 2G --tmpdir=${result}/${sn}/aligned -o ${result}/${sn}/aligned/${sn}_normal_sorted.bam /dev/stdin

#Prevent the number of open file handles on Linux from exceeding the limit and causing an error exit
ulimit -n 10240

#Mark duplicates on the sorted BAM file with sambamba
sambamba markdup \
    --tmpdir ${result}/${sn}/aligned \
    -t ${threads} ${result}/${sn}/aligned/${sn}_normal_sorted.bam ${result}/${sn}/aligned/${sn}_normal_marked.bam 

#Rename the marked BAM file index; the index file names for gatk and sambamba need to be consistent
mv  ${result}/${sn}/aligned/${sn}_normal_marked.bam.bai ${result}/${sn}/aligned/${sn}_normal_marked.bai
#Delete the sorted BAM file
rm -f ${result}/${sn}/aligned/${sn}_normal_sorted.bam*

mamba	deactivate

3. Align the tumor fastq files to the reference genome, sort, mark duplicates to obtain marked.bam, similar to step 2

bash

if [ ! -d "${envs}/${pn}.align" ]; then
	mamba create -n ${pn}.align -y bwa=${version_bwa} samtools=${version_samtools} 
fi

#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

mamba	activate ${pn}.align

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
	
    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 the reference genome sequence is not detected, download the sequence and create the index with bwa
    if [ ! -f "/opt/ref/hg38/hg38.fasta" ]; then
    	aria2c -x 16 -s 32 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
	else
        if [ -f "/opt/ref/hg38/hg38.fasta.gz.aria2" ]; then
            aria2c -x 16 -s 32 ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg38/Homo_sapiens_assembly38.fasta.gz -d /opt/ref/hg38 -d /opt/ref/hg38 -o hg38.fasta.gz
            cd /opt/ref/hg38 && gzip -d /opt/ref/hg38/hg38.fasta.gz
        fi
    fi
    
    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

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

bwa	mem \
    -t ${threads} -M \
    -R "@RG\\tID:${sn}_tumor\\tLB:${sn}_tumor\\tPL:Illumina\\tPU:Miseq\\tSM:${sn}_tumor" \
    /opt/ref/${version_reference}/${version_reference}.fasta  ${result}/${sn}/trimmed/${sn}_tumor_R1_trimmed.fq.gz ${result}/${sn}/trimmed/${sn}_tumor_R2_trimmed.fq.gz \
    | sambamba view -S -f bam -l 0 /dev/stdin \
    | sambamba sort -t ${threads} -m 2G --tmpdir=${result}/${sn}/aligned -o ${result}/${sn}/aligned/${sn}_tumor_sorted.bam /dev/stdin
ulimit -n 10240
sambamba  markdup \
    --tmpdir ${result}/${sn}/aligned \
    -t ${threads} ${result}/${sn}/aligned/${sn}_tumor_sorted.bam ${result}/${sn}/aligned/${sn}_tumor_marked.bam
    mv  ${result}/${sn}/aligned/${sn}_tumor_marked.bam.bai ${result}/${sn}/aligned/${sn}_tumor_marked.bai
    rm  -f ${result}/${sn}/aligned/${sn}_tumor_sorted.bam*

mamba	deactivate

4. Generate a recalibration table for the BAM files above for later use by BQSR; Generates recalibration table for Base Quality Score Recalibration (BQSR)

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
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 [ ! -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/somatic/projects/${project_bed} -d /opt/ref/hg19
        if [ ! -f "/opt/ref/hg19/${project_bed}.interval_list" ]; then
            gatk BedToIntervalList \
                -I	/opt/ref/hg19/${project_bed} \
                -SD	/opt/ref/hg19/hg19.dict \
                -O	/opt/ref/hg19/${project_bed}.interval_list
        fi
    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/hg19/${project_bed}.interval_list \
            -R /opt/ref/hg19/hg19.fasta \
            -I ${result}/${sn}/aligned/${sn}_tumor_marked.bam \
            -O ${result}/${sn}/recal/${sn}_tumor_recal.table &

    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/hg19/${project_bed}.interval_list \
            -R /opt/ref/hg19/hg19.fasta \
            -I ${result}/${sn}/aligned/${sn}_normal_marked.bam \
            -O ${result}/${sn}/recal/${sn}_normal_recal.table &

    wait
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 [ ! -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/hg38/Illumina_pt2.bed > /opt/ref/hg38/Illumina_pt2.processed.bed
        mkdir  -p /opt/ref/hg38
        aria2c https://raw.githubusercontent.com/doujiangbaozi/sliverworkspace-util/main/somatic/projects/${project_bed} -d /opt/ref/hg38 -o ${project_bed}
        if [ ! -f "/opt/ref/hg38/${project_bed}.interval_list" ]; then
            gatk BedToIntervalList \
                -I	/opt/ref/hg38/${project_bed} \
                -SD	/opt/ref/hg38/hg38.dict \
                -O	/opt/ref/hg38/${project_bed}.interval_list
        fi
    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/hg38/${project_bed}.interval_list \
            -R /opt/ref/hg38/hg38.fasta \
            -I ${result}/${sn}/aligned/${sn}_tumor_marked.bam \
            -O ${result}/${sn}/recal/${sn}_tumor_recal.table &

    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/hg38/${project_bed}.interval_list \
            -R /opt/ref/hg38/hg38.fasta \
            -I ${result}/${sn}/aligned/${sn}_normal_marked.bam \
            -O ${result}/${sn}/recal/${sn}_normal_recal.table &

    wait
fi


    
mamba deactivate

5. Apply base quality score recalibration to the BAM files using the recalibration table. Since GATK's efficiency at this step is disappointing, insert size is also computed at the same time, the interval list is split (needed later for Mutect2 parallel runs), cnvkit batch is run, samtools depth is run to compute sequencing depth, and samtools flagstat is run to compute mapping rate and quality statistics

bash

mkdir -p ${result}/${sn}/bqsr
mkdir -p ${result}/${sn}/stat
mkdir -p ${result}/${sn}/cnv
mkdir -p ${result}/${sn}/interval

mamba activate ${pn}.gatk

gatk ApplyBQSR \
    --bqsr-recal-file ${result}/${sn}/recal/${sn}_tumor_recal.table \
    -L /opt/ref/${version_reference}/${project_bed}.interval_list \
    -R /opt/ref/${version_reference}/${version_reference}.fasta \
    -I ${result}/${sn}/aligned/${sn}_tumor_marked.bam \
    -O ${result}/${sn}/bqsr/${sn}_tumor_bqsr.bam &


    
gatk ApplyBQSR \
    --bqsr-recal-file ${result}/${sn}/recal/${sn}_normal_recal.table \
    -L /opt/ref/${version_reference}/${project_bed}.interval_list \
    -R /opt/ref/${version_reference}/${version_reference}.fasta \
    -I ${result}/${sn}/aligned/${sn}_normal_marked.bam \
    -O ${result}/${sn}/bqsr/${sn}_normal_bqsr.bam &


    
gatk CollectInsertSizeMetrics \
    -I ${result}/${sn}/aligned/${sn}_tumor_marked.bam \
    -O ${result}/${sn}/stat/${sn}_tumor_insertsize_metrics.txt \
    -H ${result}/${sn}/stat/${sn}_tumor_insertsize_histogram.pdf &


    
gatk CollectInsertSizeMetrics \
    -I ${result}/${sn}/aligned/${sn}_normal_marked.bam \
    -O ${result}/${sn}/stat/${sn}_normal_insertsize_metrics.txt \
    -H ${result}/${sn}/stat/${sn}_normal_insertsize_histogram.pdf &

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} &

mamba deactivate

if [ ! -d "${envs}/${pn}.cnvkit" ]; then
    mamba create -n ${pn}.cnvkit -y cnvkit=${version_cnvkit}
fi

if [ ! -f "/opt/ref/${version_reference}/refFlat.txt" ]; then
    aria2c -x 16 -s 16 http://hgdownload.soe.ucsc.edu/goldenPath/${version_reference}/database/refFlat.txt.gz -d /opt/ref/${version_reference}
    cd /opt/ref/${version_reference} && gzip -d refFlat.txt.gz
fi

mamba activate ${pn}.cnvkit

rm -f ${result}/${sn}/cnv/${sn}_reference.cnn

cnvkit.py batch \
    ${result}/${sn}/aligned/${sn}_tumor_marked.bam \
    --normal ${result}/${sn}/aligned/${sn}_normal_marked.bam \
    --method hybrid \
    --fasta  /opt/ref/${version_reference}/${version_reference}.fasta \
    --targets /opt/ref/${version_reference}/${project_bed} \
    --annotate /opt/ref/${version_reference}/refFlat.txt \
    --output-reference ${result}/${sn}/cnv/${sn}_reference.cnn \
    --output-dir ${result}/${sn}/cnv/ \
    --diagram \
    -p ${threads} &

mamba deactivate

mamba activate ${pn}.align

samtools depth -a -b /opt/ref/${version_reference}/Illumina_pt2.bed  --threads ${threads} \
    ${result}/${sn}/aligned/${sn}_tumor_marked.bam > \
    ${result}/${sn}/stat/${sn}_tumor_marked.depth  &

samtools depth -a -b /opt/ref/${version_reference}/Illumina_pt2.bed  --threads ${threads} \
    ${result}/${sn}/aligned/${sn}_normal_marked.bam > \
    ${result}/${sn}/stat/${sn}_normal_marked.depth &

samtools flagstat --threads ${threads} \
    ${result}/${sn}/aligned/${sn}_tumor_marked.bam  > \
    ${result}/${sn}/stat/${sn}_tumor_marked.flagstat   &

samtools flagstat --threads ${threads} \
    ${result}/${sn}/aligned/${sn}_normal_marked.bam > \
    ${result}/${sn}/stat/${sn}_normal_marked.flagstat &
    
mamba deactivate

wait

6. Compute pileup metrics for subsequent contamination estimation; they can also be processed in parallel using the split interval lists and merged afterwards.

bash

#A major official pitfall: the chromosome coordinates of the default-provided small_exac_common_3_b37.vcf.gz are not prefixed with "chr" but are plain numbers
mamba activate ${pn}.gatk

if [ "${version_reference}" == hg19 ]; then
    echo	"USE reference Version : ${version_reference}"
    #Here is a huge pitfall: the reference files downloaded from the Mutect2 directory of the Broad Institute FTP bundle site are inconsistent with the coordinate system of the reference genome sequence downloaded from the same site; the reference genome sequence uses the chr1 style, whereas this af-only-gnomad uses the 1, 2, 3 style, so a script is needed to process it
    if [ ! -f "/opt/ref/hg19/small_exac_common_3_b37.processed.vcf.gz" ]; then
        aria2c ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/Mutect2/GetPileupSummaries/small_exac_common_3_b37.vcf.gz -d /opt/ref/hg19
    else
        if [ -f "/opt/ref/hg19/small_exac_common_3_b37.vcf.gz.aria2" ]; then
            aria2c ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/Mutect2/GetPileupSummaries/small_exac_common_3_b37.vcf.gz -c -d /opt/ref/hg19
        fi
    fi
    if [ ! -f "/opt/ref/${version_reference}/small_exac_common_3_b37.processed.vcf.gz" ]; then
		if [ ! -f "${envs}/VcfProcessUtil.py" ]; then
			aria2c https://raw.fgit.cf/doujiangbaozi/sliverworkspace-util/main/somatic/VcfProcessUtil.py -d ${envs}/
			#aria2c https://raw.githubusercontent.com/doujiangbaozi/sliverworkspace-util/main/somatic/VcfProcessUtil.py -d ${envs}/
			chmod a+x ${envs}/VcfProcessUtil.py
		fi
		mamba activate ${pn}.cnvkit
		${envs}/VcfProcessUtil.py \
			-f /opt/ref/${version_reference}/small_exac_common_3_b37.vcf.gz \
			-o /opt/ref/${version_reference}/small_exac_common_3_b37.processed.vcf
		mamba deactivate
		mamba activate ${pn}.gatk
		cd /opt/ref/${version_reference}
		bgzip -f --threads ${threads} small_exac_common_3_b37.processed.vcf
		tabix -f small_exac_common_3_b37.processed.vcf.gz
		mamba deactivate
    fi
	mamba activate ${pn}.gatk
    for i in `ls ${result}/${sn}/interval/*.interval_list`;
    do
        echo $i
        gatk GetPileupSummaries \
            -R /opt/ref/${version_reference}/${version_reference}.fasta \
            -I ${result}/${sn}/bqsr/${sn}_tumor_bqsr.bam \
            -O ${i%.*}-tumor-pileups.table \
            -V /opt/ref/${version_reference}/small_exac_common_3_b37.processed.vcf.gz \
            -L $i \
            --interval-set-rule INTERSECTION &

        gatk GetPileupSummaries \
            -R /opt/ref/${version_reference}/${version_reference}.fasta \
            -I ${result}/${sn}/bqsr/${sn}_normal_bqsr.bam \
            -O ${i%.*}-normal-pileups.table \
            -V /opt/ref/${version_reference}/small_exac_common_3_b37.processed.vcf.gz \
            -L $i \
            --interval-set-rule INTERSECTION &

    done
    wait
    mamba deactivate
elif [ "${version_reference}" == "hg38" ]; then
    echo	"USE reference Version : ${version_reference}"
    if [ ! -f "/opt/ref/${version_reference}/small_exac_common_3.hg38.vcf.gz" ]; then
        aria2c ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/Mutect2/GetPileupSummaries/small_exac_common_3.hg38.vcf.gz -d /opt/ref/${version_reference}
    else
        if [ -f "/opt/ref/hg19/small_exac_common_3_b37.vcf.gz.aria2" ]; then
            aria2c ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/Mutect2/GetPileupSummaries/small_exac_common_3.hg38.vcf.gz -c -d /opt/ref/${version_reference}
        fi
    fi
    if [ ! -f "/opt/ref/${version_reference}/small_exac_common_3.hg38.vcf.gz.tbi" ]; then
        aria2c ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/Mutect2/GetPileupSummaries/small_exac_common_3.hg38.vcf.gz.tbi -d /opt/ref/${version_reference}
    else
        if [ -f "/opt/ref/hg19/small_exac_common_3_b37.vcf.gz.tbi.aria2" ]; then
            aria2c ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/Mutect2/GetPileupSummaries/small_exac_common_3.hg38.vcf.gz.tbi -c -d /opt/ref/${version_reference}
        fi
    fi
    mamba activate ${pn}.gatk
    for i in `ls ${result}/${sn}/interval/*.interval_list`;
    do
        echo $i
        gatk GetPileupSummaries \
            -R /opt/ref/${version_reference}/${version_reference}.fasta \
            -I ${result}/${sn}/bqsr/${sn}_tumor_bqsr.bam \
            -O ${i%.*}-tumor-pileups.table \
            -V /opt/ref/${version_reference}/small_exac_common_3.hg38.vcf.gz \
            -L $i \
            --interval-set-rule INTERSECTION &

        gatk GetPileupSummaries \
            -R /opt/ref/${version_reference}/${version_reference}.fasta \
            -I ${result}/${sn}/bqsr/${sn}_normal_bqsr.bam \
            -O ${i%.*}-normal-pileups.table \
            -V /opt/ref/${version_reference}/small_exac_common_3.hg38.vcf.gz \
            -L $i \
            --interval-set-rule INTERSECTION &

    done
    wait
    mamba deactivate
fi


mamba activate ${pn}.gatk
tables=
for i in `ls ${result}/${sn}/interval/*-tumor-pileups.table`;
do
    tables="$tables -I $i"
done

gatk GatherPileupSummaries \
    --sequence-dictionary /opt/ref/${version_reference}/${version_reference}.dict \
    $tables \
    -O ${result}/${sn}/stat/${sn}_tumor_pileups.table

nctables=
for i in `ls ${result}/${sn}/interval/*-normal-pileups.table`;
do
    nctables="$nctables -I $i"
done

gatk GatherPileupSummaries \
    --sequence-dictionary /opt/ref/${version_reference}/${version_reference}.dict \
    $nctables \
    -O ${result}/${sn}/stat/${sn}_normal_pileups.table
    
mamba deactivate

7. Use the GetPileupSummaries results to estimate cross-sample contamination; the results are used later by FilterMutectCalls to filter the Mutect2 output

bash

mamba activate ${pn}.gatk

gatk CalculateContamination \
    -matched ${result}/${sn}/stat/${sn}_normal_pileups.table \
    -I ${result}/${sn}/stat/${sn}_tumor_pileups.table \
    -O ${result}/${sn}/stat/${sn}_contamination.table
    
mamba deactivate

8. Call variants with Mutect2 using the split interval lists, then merge the results; meanwhile run manta in parallel to call SVs

bash

mkdir -p ${result}/${sn}/sv
mkdir -p ${result}/${sn}/snv

mamba activate ${pn}.gatk

if [ "${version_reference}" == hg19 ]; then
    echo	"USE reference Version : ${version_reference}"
    #Here is a huge pitfall: the reference files downloaded from the Mutect2 directory of the Broad Institute FTP bundle site are inconsistent with the coordinate system of the reference genome sequence downloaded from the same site; the reference genome sequence uses the chr1 style, whereas this af-only-gnomad uses the 1, 2, 3 style, so a script is needed to process it; hg38 does not seem to have this problem—are the hg19 data simply no longer maintained?
    if [ ! -f "/opt/ref/hg19/af-only-gnomad.raw.sites.b37.vcf.gz" ]; then
        aria2c -x 16 -s 32 ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/Mutect2/af-only-gnomad.raw.sites.b37.vcf.gz -d /opt/ref/hg19
    else
        if [ -f "/opt/ref/hg19/af-only-gnomad.raw.sites.b37.vcf.gz.aria2" ]; then
            aria2c -x 16 -s 32 ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/Mutect2/af-only-gnomad.raw.sites.b37.vcf.gz -c -d /opt/ref/hg19
        fi
        if [ ! -f "${envs}/VcfProcessUtil.py" ]; then
            aria2c https://raw.fgit.cf/doujiangbaozi/sliverworkspace-util/main/somatic/VcfProcessUtil.py -d ${envs}/
            #aria2c https://raw.githubusercontent.com/doujiangbaozi/sliverworkspace-util/main/somatic/VcfProcessUtil.py -d ${envs}/
            chmod a+x ${envs}/VcfProcessUtil.py
        fi
        mamba activate ${pn}.cnvkit
        ${envs}/VcfProcessUtil.py \
            -f /opt/ref/hg19/af-only-gnomad.raw.sites.b37.vcf.gz \
            -o /opt/ref/hg19/af-only-gnomad.raw.sites.b37.processed.vcf
        mamba deactivate
        mamba activate ${pn}.gatk
        cd /opt/ref/hg19
        bgzip -f --threads ${threads} af-only-gnomad.raw.sites.b37.processed.vcf
        tabix -f af-only-gnomad.raw.sites.b37.processed.vcf.gz
        mamba deactivate
    fi
elif [ "${version_reference}" == "hg38" ]; then
    echo	"USE reference Version : ${version_reference}"
    if [ ! -f "/opt/ref/${version_reference}/af-only-gnomad.hg38.vcf.gz" ]; then
        aria2c -x 16 -s 32 ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/Mutect2/af-only-gnomad.hg38.vcf.gz -d /opt/ref/${version_reference}
    else
        if [ -f "/opt/ref/${version_reference}/af-only-gnomad.hg38.vcf.gz.aria2" ]; then
            aria2c -x 16 -s 32 ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/Mutect2/af-only-gnomad.hg38.vcf.gz -c -d /opt/ref/${version_reference}
        fi
    fi
    if [ ! -f "/opt/ref/${version_reference}/af-only-gnomad.hg38.vcf.gz.tbi" ]; then
        aria2c -x 16 -s 32 ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/Mutect2/af-only-gnomad.hg38.vcf.gz.tbi -d /opt/ref/${version_reference}
    else
        if [ -f "/opt/ref/${version_reference}/af-only-gnomad.hg38.vcf.gz.tbi.aria2" ]; then
            aria2c -x 16 -s 32 ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/Mutect2/af-only-gnomad.hg38.vcf.gz.tbi -c -d /opt/ref/${version_reference}
        fi
    fi
fi


    
mamba activate ${pn}.gatk
if [ ! -f "/opt/ref/${version_reference}/${project_bed}.gz" ]; then
    bgzip -f -c /opt/ref/${version_reference}/${project_bed} > /opt/ref/hg19/${project_bed}.gz
    tabix -f -p bed /opt/ref/${version_reference}/${project_bed}.gz
else
    if [ ! -f "/opt/ref/${version_reference}/${project_bed}.gz.tbi" ]; then
        tabix -f -p bed /opt/ref/${version_reference}/${project_bed}.gz
    fi
fi
mamba deactivate

if [ ! -d "${envs}/${pn}.manta" ]; then
    mamba create -n ${pn}.manta -y manta=1.6.0
fi

mamba activate ${pn}.manta

rm -f ${result}/${sn}/sv/runWorkflow.py*
configManta.py  \
    --normalBam ${result}/${sn}/bqsr/${sn}_normal_bqsr.bam \
    --tumorBam  ${result}/${sn}/bqsr/${sn}_tumor_bqsr.bam \
    --referenceFasta /opt/ref/${version_reference}/${version_reference}.fasta \
    --exome \
    --callRegions /opt/ref/${version_reference}/${project_bed}.gz \
    --runDir ${result}/${sn}/sv

rm -rf ${result}/${sn}/sv/workspace
python ${result}/${sn}/sv/runWorkflow.py -m local -j ${threads} &

mamba deactivate

mamba activate ${pn}.gatk

if [ "${version_reference}" == hg19 ]; then
    echo	"USE reference Version : ${version_reference}"
    rm -f ${result}/${sn}/snv/vcf-file.list
    touch ${result}/${sn}/snv/vcf-file.list
    for i in `ls ${result}/${sn}/interval/*.interval_list`;
    do
        rm -f ${i%.*}_bqsr.vcf.gz
        gatk Mutect2 \
            -R /opt/ref/${version_reference}/${version_reference}.fasta \
            -I ${result}/${sn}/bqsr/${sn}_tumor_bqsr.bam  -tumor  ${sn}_tumor  \
            -I ${result}/${sn}/bqsr/${sn}_normal_bqsr.bam -normal ${sn}_normal \
            -L $i \
            -O ${i%.*}_bqsr.vcf.gz \
            --max-mnp-distance 10 \
            --germline-resource /opt/ref/${version_reference}/af-only-gnomad.raw.sites.b37.processed.vcf.gz \
            --native-pair-hmm-threads ${threads} &
        echo ${i%.*}_bqsr.vcf.gz >> ${result}/${sn}/snv/vcf-file.list
    done
elif [ "${version_reference}" == "hg38" ]; then
    echo	"USE reference Version : ${version_reference}"
    rm -f ${result}/${sn}/snv/vcf-file.list
    touch ${result}/${sn}/snv/vcf-file.list
    for i in `ls ${result}/${sn}/interval/*.interval_list`;
    do
        rm -f ${i%.*}_bqsr.vcf.gz
        gatk Mutect2 \
            -R /opt/ref/${version_reference}/${version_reference}.fasta \
            -I ${result}/${sn}/bqsr/${sn}_tumor_bqsr.bam  -tumor  ${sn}_tumor  \
            -I ${result}/${sn}/bqsr/${sn}_normal_bqsr.bam -normal ${sn}_normal \
            -L $i \
            -O ${i%.*}_bqsr.vcf.gz \
            --max-mnp-distance 10 \
            --germline-resource /opt/ref/${version_reference}/af-only-gnomad.hg38.vcf.gz \
            --native-pair-hmm-threads ${threads} &
        echo ${i%.*}_bqsr.vcf.gz >> ${result}/${sn}/snv/vcf-file.list
    done
fi
wait

    
rm -f ${result}/${sn}/snv/${sn}_bqsr.vcf.gz.stats
stats=
for z in `ls ${result}/${sn}/interval/*_bqsr.vcf.gz.stats`;
do
    stats="$stats -stats $z"
done

gatk MergeMutectStats $stats \
    -O ${result}/${sn}/snv/${sn}_bqsr.vcf.gz.stats

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

9. FilterMutectCalls to filter the Mutect2 variant results

bash

mamba activate ${pn}.gatk

gatk FilterMutectCalls \
    --max-events-in-region ${event} \
    --contamination-table ${result}/${sn}/stat/${sn}_contamination.table \
    -R /opt/ref/${version_reference}/${version_reference}.fasta \
    -V ${result}/${sn}/snv/${sn}_bqsr.vcf.gz \
    -O ${result}/${sn}/snv/${sn}_filtered.vcf.gz
    
mamba deactivate

10. Annotate the filtered 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

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

if [ "${version_reference}" == hg19 ]; then
echo	"USE reference Version : ${version_reference}"
    #Check whether the VEP annotation database exists; if not, download it first
    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

    if [ ! -d "/opt/ref/vep-cache/homo_sapiens_refseq/${version_vep_cache}_GRCh37" ]; then
        aria2c -x 16 -s 48 http://ftp.ensembl.org/pub/release-${version_vep_cache}/variation/vep/homo_sapiens_refseq_vep_${version_vep_cache}_GRCh37.tar.gz -d /opt/ref/
        tar -zxvf /opt/ref/homo_sapiens_refseq_vep_${version_vep_cache}_GRCh37.tar.gz -C /opt/ref/vep-cache/
    elif [ -f "/opt/ref/homo_sapiens_refseq_vep_${version_vep_cache}_GRCh37.tar.gz.aria2" ]; then
        aria2c -x 16 -s 48 http://ftp.ensembl.org/pub/release-${version_vep_cache}/variation/vep/homo_sapiens_refseq_vep_${version_vep_cache}_GRCh37.tar.gz -c -d /opt/ref/
        tar -zxvf /opt/ref/homo_sapiens_refseq_vep_${version_vep_cache}_GRCh37.tar.gz -C /opt/ref/vep-cache/
    fi
    
    mamba activate ${pn}.vep

    mkdir -p ${result}/${sn}/annotation
    vep \
        -i ${result}/${sn}/snv/${sn}_filtered.vcf.gz  \
        -o ${result}/${sn}/annotation/${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 \
        --fasta /opt/ref/${version_reference}/${version_reference}.fasta \
        --refseq \
        --force_overwrite \
        --format vcf \
        --tab \
        --shift_3prime 1  \
        --offline

    mamba deactivate
    
elif [ "${version_reference}" == "hg38" ]; then
    echo	"USE reference Version : ${version_reference}"
    #Check whether the VEP annotation database exists; if not, download it first
    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

    if [ ! -d "/opt/ref/vep-cache/homo_sapiens_refseq/${version_vep_cache}_GRCh38" ]; then
        aria2c -x 16 -s 48 http://ftp.ensembl.org/pub/release-${version_vep_cache}/variation/vep/homo_sapiens_refseq_vep_${version_vep_cache}_GRCh38.tar.gz -d /opt/ref/
        tar -zxvf /opt/ref/homo_sapiens_refseq_vep_${version_vep_cache}_GRCh38.tar.gz -C /opt/ref/vep-cache/
    elif [ -f "/opt/ref/homo_sapiens_refseq_vep_${version_vep_cache}_GRCh38.tar.gz.aria2" ]; then
        aria2c -x 16 -s 48 http://ftp.ensembl.org/pub/release-${version_vep_cache}/variation/vep/homo_sapiens_refseq_vep_${version_vep_cache}_GRCh38.tar.gz -c -d /opt/ref/
        tar -zxvf /opt/ref/homo_sapiens_refseq_vep_${version_vep_cache}_GRCh38.tar.gz -C /opt/ref/vep-cache/
    fi
    
    mamba activate ${pn}.vep

    mkdir -p ${result}/${sn}/annotation
    vep \
        -i ${result}/${sn}/snv/${sn}_filtered.vcf.gz  \
        -o ${result}/${sn}/annotation/${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 \
        --fasta /opt/ref/${version_reference}/${version_reference}.fasta \
        --refseq \
        --force_overwrite \
        --format vcf \
        --tab \
        --shift_3prime 1  \
        --offline

    mamba deactivate
fi

11. Use a script to process the annotation results and the filtered VCF results, outputting a data table matching the external quality assessment (EQA) required format

bash

mamba activate ${pn}.cnvkit

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

${envs}/MatchedSnvVepAnnotationFilter.py \
    -e normal_artifact   \
    -e germline   \
    -i strand_bias   \
    -i clustered_events   \
    --min-vaf=${snv_vaf} \
    --min-tlod=${snv_tlod} \
    --min-depth=${snv_depth} \
    -v ${result}/${sn}/snv/${sn}_filtered.vcf.gz   \
    -a ${result}/${sn}/annotation/${sn}_filtered_vep.tsv   \
    -o ${result}/${sn}/annotation/${sn}.result.SNV.tsv

mamba deactivate

12. Use cnvkit's provided tools to output scatter plots and heatmaps

bash

mamba activate ${pn}.cnvkit
    
    cnvkit.py scatter ${result}/${sn}/cnv/${sn}_tumor_marked.cnr \
    -s ${result}/${sn}/cnv/${sn}_tumor_marked.cns \
    -i ' ' \
    -n ${sn}_normal \
    -o ${result}/${sn}/cnv/${sn}_cnv_scatter.png -t  &&

cnvkit.py heatmap ${result}/${sn}/cnv/${sn}_tumor_marked.cns \
-o ${result}/${sn}/cnv/${sn}_cnv_heatmap.png

mamba deactivate

13. Use cnvkit call to infer copy numbers from the cnvkit batch output results

bash

mamba activate ${pn}.cnvkit

cnvkit.py call ${result}/${sn}/cnv/${sn}_tumor_marked.cns \
	-o ${result}/${sn}/cnv/${sn}_tumor_marked.call.cns

mamba deactivate

14. Write a script to process the cnvkit output, computing CNV genes, exon positions, gain/loss, and copy number

bash

mamba activate ${pn}.cnvkit

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

if [ ! -f "/opt/ref/${version_reference}/refGene.txt" ]; then
    aria2c -x 16 -s 16 http://hgdownload.cse.ucsc.edu/goldenPath/${version_reference}/database/refGene.txt.gz -d /opt/ref/${version_reference} -o refGene.txt.gz
    cd /opt/ref/${version_reference} && gzip -d refGene.txt.gz
fi

python ${envs}/CnvAnnotationFilter.py  \
    -r /opt/ref/${version_reference}/refGene.txt \
    -i ${cnv_min} \
    -x ${cnv_max} \
    -D ${cnv_dep} \
    -f ${result}/${sn}/cnv/${sn}_tumor_marked.call.cns \
    -o ${result}/${sn}/cnv/${sn}.result.CNV.tsv

mamba deactivate

15. Write a script to process the manta output, obtaining the final SV results with start positions, genes, frequencies, etc.

bash

mamba activate ${pn}.cnvkit

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

if [ ! -f "/opt/ref/${version_reference}/refGene.txt" ]; then
    aria2c -x 16 -s 16 http://hgdownload.cse.ucsc.edu/goldenPath/${version_reference}/database/refGene.txt.gz -d /opt/ref/${version_reference} -o refGene.txt.gz
    cd /opt/ref/${version_reference} && gzip -d refGene.txt.gz
fi

${envs}/SvAnnotationFilter.py \
    -r /opt/ref/${version_reference}/refGene.txt \
    -s ${sv_score} \
    -f ${result}/${sn}/sv/results/variants/somaticSV.vcf.gz \
    -o ${result}/${sn}/sv/${sn}.result.SV.tsv

mamba deactivate

16. Based on the previous outputs of fastp, samtools depth, samtools flagstat, gatk CollectInsertSizeMetrics, etc., provide a comprehensive QC summary

bash

mamba activate ${pn}.cnvkit

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

${envs}/MatchedQcProcessor.py  --bed /opt/ref/${version_reference}/${project_bed} \
    --out ${result}/${sn}/stat/${sn}.result.QC.tsv \
    --sample-fastp=${result}/${sn}/trimmed/${sn}_tumor_fastp.json \
    --sample-depth=${result}/${sn}/stat/${sn}_tumor_marked.depth \
    --sample-flagstat=${result}/${sn}/stat/${sn}_tumor_marked.flagstat \
    --sample-insertsize=${result}/${sn}/stat/${sn}_tumor_insertsize_metrics.txt \
    --normal-fastp=${result}/${sn}/trimmed/${sn}_normal_fastp.json \
    --normal-depth=${result}/${sn}/stat/${sn}_normal_marked.depth \
    --normal-flagstat=${result}/${sn}/stat/${sn}_normal_marked.flagstat  \
    --normal-insertsize=${result}/${sn}/stat/${sn}_normal_insertsize_metrics.txt

mamba deactivate

mamba activate ${pn}.fastp

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

mamba deactivate

Final output

File name	Remarks
report.pdf	Graphical analysis report designed by sliverworkspace, printed and exported as PDF
1701/qc/multiqc_report.html	multiqc aggregated report
1701/1701.result.SNV.tsv	Final SNV variant results
1701/1701/cnv/1701_cnv_heatmap.png	CNV results heatmap
1701/cnv/1701_cnv_scatter.png	CNV results scatter plot
1701/cnv/1701.result.CNV.tsv	Final CNV results
1701.result.SV.tsv	Final SV results
1701.result.QC.tsv	Final QC results
1701/qc/multiqc_report.html	multiqc report

GATK Somatic SNV+Indel+CNV+SV ​

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 ​

Overall Pipeline Overview ​

Docker image pull and deployment configuration ​

docker-compose.yml configuration file ​

Deploying and running the docker image ​

Test data ​

Variable names ​

Detailed Pipeline/workflow steps: ​

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

2. Align the normal fastq files to the reference genome, sort, mark duplicates to obtain marked.bam ​

3. Align the tumor fastq files to the reference genome, sort, mark duplicates to obtain marked.bam, similar to step 2 ​

4. Generate a recalibration table for the BAM files above for later use by BQSR; Generates recalibration table for Base Quality Score Recalibration (BQSR) ​

6. Compute pileup metrics for subsequent contamination estimation; they can also be processed in parallel using the split interval lists and merged afterwards. ​

7. Use the GetPileupSummaries results to estimate cross-sample contamination; the results are used later by FilterMutectCalls to filter the Mutect2 output ​

8. Call variants with Mutect2 using the split interval lists, then merge the results; meanwhile run manta in parallel to call SVs ​

9. FilterMutectCalls to filter the Mutect2 variant results ​

10. Annotate the filtered results with VEP ​

11. Use a script to process the annotation results and the filtered VCF results, outputting a data table matching the external quality assessment (EQA) required format ​

12. Use cnvkit's provided tools to output scatter plots and heatmaps ​

13. Use cnvkit call to infer copy numbers from the cnvkit batch output results ​

14. Write a script to process the cnvkit output, computing CNV genes, exon positions, gain/loss, and copy number ​

15. Write a script to process the manta output, obtaining the final SV results with start positions, genes, frequencies, etc. ​

16. Based on the previous outputs of fastp, samtools depth, samtools flagstat, gatk CollectInsertSizeMetrics, etc., provide a comprehensive QC summary ​

Final output ​

GATK Somatic SNV+Indel+CNV+SV

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

Overall Pipeline Overview

Docker image pull and deployment configuration

docker-compose.yml configuration file

Deploying and running the docker image

Test data

Variable names

Detailed Pipeline/workflow steps:

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

2. Align the normal fastq files to the reference genome, sort, mark duplicates to obtain marked.bam

3. Align the tumor fastq files to the reference genome, sort, mark duplicates to obtain marked.bam, similar to step 2

4. Generate a recalibration table for the BAM files above for later use by BQSR; Generates recalibration table for Base Quality Score Recalibration (BQSR)

6. Compute pileup metrics for subsequent contamination estimation; they can also be processed in parallel using the split interval lists and merged afterwards.

7. Use the GetPileupSummaries results to estimate cross-sample contamination; the results are used later by FilterMutectCalls to filter the Mutect2 output

8. Call variants with Mutect2 using the split interval lists, then merge the results; meanwhile run manta in parallel to call SVs

9. FilterMutectCalls to filter the Mutect2 variant results

10. Annotate the filtered results with VEP

11. Use a script to process the annotation results and the filtered VCF results, outputting a data table matching the external quality assessment (EQA) required format

12. Use cnvkit's provided tools to output scatter plots and heatmaps

13. Use cnvkit call to infer copy numbers from the cnvkit batch output results

14. Write a script to process the cnvkit output, computing CNV genes, exon positions, gain/loss, and copy number

15. Write a script to process the manta output, obtaining the final SV results with start positions, genes, frequencies, etc.

16. Based on the previous outputs of fastp, samtools depth, samtools flagstat, gatk CollectInsertSizeMetrics, etc., provide a comprehensive QC summary

Final output