WES Pipeline ReadMe

Whole Exome Sequencing (WES) Documentation

Whole Exome Sequencing Pipeline (--workflow wes)

•	Step 1: FastP read and adapter trimming   
•	Step 2: Get Read Group Information  
•	Step 3: BWA-MEM Alignment  
•	Step 4: Variant Preprocessing (Picard sortsam/mark duplicates)  

If Mouse:

•	Step 6: Variant Calling (GATK Haplotype Caller)  
•	Step 7: Variant Filtration (GATK)  
•	Step 8: Post Variant Calling Processing (spnEff, GATK VariantAnnotator, snpSift Extract Fields)  

If Human:

•	Step 5: Variant Pre-Processing – Part 2 (GATK Base Recalibrator Apply BQSR)  
•	Step 6: Variant Calling (GATK Haplotype Caller)  
•	Step 7: Variant Filtration (GATK Select Variants and Variant Filtration)  SNP and Indel  
•	Step 8: Post Variant Calling Processing - Part 1 (Cosmic, SnpEff, SnpSift) SNP and Indel  
•	Step 9: Post Variant Calling Processing - Part 2 (GATK Merge VCF)  

Both Mouse and Human:

•	Step 10: Picard Collect HS Metrics  
•	Step 11: Aggregate Stats  
•	Step 12: MultiQC report generation 

Mouse WES Flowchart

flowchart TD
	p1((Sample))
	p2[FASTP]
	p3[FASTQC]
	p6[BWA_MEM]
    o1([Genomic BAM]):::output
    p7[PICARD_SORTSAM]
	p8[PICARD_MARKDUPLICATES]
	p11[PICARD_COLLECTHSMETRICS]
	p14[GATK_HAPLOTYPECALLER]
    o3([Raw Variant Calls]):::output
    p16[SNPSIFT_ANNOTATE_DBSNP]
	p17[GATK_VARIANTFILTRATION]
    p18[GATK_SELECTVARIANTS_SNP]
	p19[GATK_SELECTVARIANTS_INDEL]
	p20[SNPEFF]
	p21[SNPEFF_ONEPERLINE]
	p22[SNPSIFT_EXTRACTFIELDS]
    o4([Filtered Unannotated SNPs]):::output
    o5([Filtered Unannotated InDELs]):::output
    o6([Filtered Annotated VCF]):::output
    o7([Variant Table]):::output
	p33[MULTIQC]
    o8([MultiQC Report]):::output
	
    p1 --> |Raw Reads| p2
    subgraph alignment [  ]
    p2 --> p6
    p6 --> p7
    p7 --> p8
    p8 --> o1
    end

    subgraph calling [  ]
    o1 --> p14
    p14 --> o3
    o3 --> p16
    p16 --> p17
    p17 --> p18
    p17 --> p19
    p18 --> o4
    p19 --> o5
    p17 --> p20
    p20 --> p21
    p21 --> p22
    p20 --> o6
    p22 --> o7

    end
    subgraph qc [  ]
    p2 --> p3
    o1 --> p11
    p2 --> p33
    p3 --> p33
    p8 --> p33
    p11 --> p33
    p33 --> o8
    end

classDef output fill:#90aaff,stroke:#6c8eff,stroke-width:2px,color:#000000
style alignment stroke:#333,stroke-width:2px
style calling stroke:#333,stroke-width:2px
style qc stroke:#333,stroke-width:2px

Loading

Human WES Flowchart

flowchart TD
	p1((Sample))
	p2[FASTP]
	p3[FASTQC]
	%% p5[READ_GROUPS]
	p6[BWA_MEM]
    o1([Genomic BAM]):::output
    p7[PICARD_SORTSAM]
	p8[PICARD_MARKDUPLICATES]
	p9[GATK_BASERECALIBRATOR]
	p10[GATK_APPLYBQSR]
	p11[PICARD_COLLECTHSMETRICS]
	%% p13[GATK_GETSAMPLENAME]
	p14[GATK_HAPLOTYPECALLER]
    o3([Raw Variant Calls]):::output
    p16[GATK_SELECTVARIANTS_SNP]
	p17[GATK_VARIANTFILTRATION_SNP]
	p18[GATK_SELECTVARIANTS_INDEL]
	p19[GATK_VARIANTFILTRATION_INDEL]
	p20[SNPSIFT_ANNOTATE_SNP_DBSNP]
	p21[SNPSIFT_ANNOTATE_SNP_COSMIC]
	p22[SNPEFF_SNP]
	p23[SNPSIFT_DBNSFP_SNP]
	p24[SNPEFF_ONEPERLINE_SNP]
	p25[SNPSIFT_ANNOTATE_INDEL_DBSNP]
	p26[SNPSIFT_ANNOTATE_INDEL_COSMIC]
	p27[SNPEFF_INDEL]
	p28[SNPSIFT_DBNSFP_INDEL]
	p29[SNPEFF_ONEPERLINE_INDEL]
	p30[GATK_MERGEVCF_UNANNOTATED]
    o4([Filtered Unannoated VCF]):::output
	p31[GATK_MERGEVCF_ANNOTATED]
    o5([Filtered Annotated VCF]):::output
	p32[SNPSIFT_EXTRACTFIELDS]
    o6([Variant Table]):::output
	p33[MULTIQC]
    o7([MultiQC Report]):::output
    p1 --> |Raw Reads| p2
    subgraph alignment [  ]
    p2 --> p6
    %% p2 --> p5
    %% p5 --> p6
    p6 --> p7
    p7 --> p8
    p8 --> p9
    p9 --> p10
    p10 --> o1
    end

    subgraph calling [  ]
    %% o1 --> p13
    o1 --> p14
    %% p13 --> p14
    p14 --> o3
    o3 --> p16
    o3 --> p18
    p16 --> p17
    p18 --> p19
    p17 --> p20
    p20 --> p21
    p21 --> p22
    p22 --> p23
    p23 --> p24
    p19 --> p25
    p25 --> p26
    p26 --> p27 
    p27 --> p28
    p28 --> p29
    p17 --> p30
    p19 --> p30
    p30 --> o4
    p24 --> p31
    p29 --> p31
    p31 --> o5
    o4 --> p32
    p32 --> o6
    end
    subgraph qc [  ]
    p2 --> p3
    o1 --> p11
    p2 --> p33
    p3 --> p33
    p8 --> p33
    p9 --> p33
    p11 --> p33
    p33 --> o7
    end

classDef output fill:#90aaff,stroke:#6c8eff,stroke-width:2px,color:#000000
style alignment stroke:#333,stroke-width:2px
style calling stroke:#333,stroke-width:2px
style qc stroke:#333,stroke-width:2px

Loading

Parameters for WES Pipeline

• --pubdir

  • Default: /<PATH>
  • Comment: The directory that the saved outputs will be stored.

• --organize_by

  • Default: sample
  • Comment: How to organize the output folder structure. Options: sample or analysis.

• --cacheDir

  • Default: /projects/omics_share/meta/containers
  • Comment: This is directory that contains cached Singularity containers. JAX users should not change this parameter.

• -w

  • Default: /<PATH>
  • Comment: The directory that all intermediary files and nextflow processes utilize. This directory can become quite large. This should be a location on /fastscratch or other directory with ample storage.

• --sample_folder

  • Default: /<PATH>
  • Comment: The path to the folder that contains all the samples to be run by the pipeline. The files in this path can also be symbolic links.

• --extension

  • Default: .fastq.gz
  • Comment: The expected extension for the input read files.

• --pattern

  • Default: "*_R{1,2}*"
  • Comment: The expected R1 / R2 matching pattern. The default value will match reads with names like this READ_NAME_R1_MoreText.fastq.gz or READ_NAME_R1.fastq.gz

• --read_type

  • Default: PE
  • Comment: Options: PE and SE. Default: PE. Type of reads: paired end (PE) or single end (SE).

• --concat_lanes

  • Default: false
  • Comment: Options: false and true. Default: false. If this boolean is specified, FASTQ files will be concatenated by sample. Used in cases where samples are divided across individual sequencing lanes.

• --csv_input

  • Default: null
  • Comment: Provide a CSV manifest file with the header: "sampleID,lane,fastq_1,fastq_2". See below for an example file. Fastq_2 is optional and used only in PE data. Fastq files can either be absolute paths to local files, or URLs to remote files. If remote URLs are provided, * --download_data can be specified.

• --download_data

  • Default: null
  • Comment: Requires * --csv_input. When specified, read data in the CSV manifest will be downloaded from provided URLs with Aria2.

• --gen_org

  • Default: mouse
  • Comment: Options: mouse and human.

• --genome_build

  • Default: GRCm38
  • Comment: Mouse specific. Options: GRCm38 or GRCm39. If gen_org == human, build defaults to GRCh38.

• --ref_fa

  • Default: '/projects/omics_share/mouse/GRCm38/genome/sequence/ensembl/v102/Mus_musculus.GRCm38.dna.toplevel.fa'
  • Comment: The reference fasta to be used throughout the process for alignment as well as any downstream analysis, points to human reference when * --gen_org human. JAX users should not change this parameter.

• --ref_fa_indices

  • Default: '/projects/omics_share/mouse/GRCm38/genome/indices/ensembl/v102/bwa/Mus_musculus.GRCm38.dna.toplevel.fa'
  • Comment: Pre-compiled BWA index files, points to human reference when * --gen_org human. JAX users should not change this parameter.

• --quality_phred

  • Default: 15
  • Comment: The quality value that is required for a base to pass. Default: 15 which is a phred quality score of >=Q15.

• --unqualified_perc

  • Default: 40
  • Comment: Percent of bases that are allowed to be unqualified (0~100). Default: 40 which is 40%.

• --detect_adapter_for_pe

  • Default: false
  • Comment: If true, adapter auto-detection is used for paired end data. By default, paired-end data adapter sequence auto-detection is disabled as the adapters can be trimmed by overlap analysis. However, --detect_adapter_for_pe will enable it. Fastp will run a little slower if you specify the sequence adapters or enable adapter auto-detection, but usually result in a slightly cleaner output, since the overlap analysis may fail due to sequencing errors or adapter dimers.

• --target_gatk

  • Default: '/projects/omics_share/mouse/GRCm38/supporting_files/mm10Exome_v4_12-19.1.mm10.baits_and_targets.bare.bed'
  • Comment: A bed file with WES target intervals as defined in the capture array used in the data, points to human bed file when * --gen_org human. NOTE: This file MUST reflect the capture array used to generate your data.

• --target_picard

  • Default: '/projects/omics_share/mouse/GRCm38/supporting_files/mm10Exome_v4_12-19.1.mm10.baits_and_targets.bare.picard.interval_list'
  • Comment: A GATK interval file covering WES target intervals. Used in calculating coverage metrics. NOTE: This file MUST reflect the capture array used to generate your data.

• --bait_picard

  • Default: '/projects/omics_share/mouse/GRCm38/supporting_files/mm10Exome_v4_12-19.1.mm10.baits_and_targets.bare.picard.interval_list'
  • Comment: A GATK interval file covering WES target intervals. Used in calculating coverage metrics. This file can be the same as the interval file, NOTE: This file MUST reflect the capture array used to generate your data.

• --mismatch_penalty

  • Default: -B 8
  • Comment: The BWA penalty for a mismatch.

• --call_val

  • Default: 50
  • Comment: The minimum phred-scaled confidence threshold at which variants should be called.

• --ploidy_val

  • Default: "-ploidy 2"
  • Comment: Sample ploidy

• --dbSNP

  • Default: '/projects/omics_share/mouse/GRCm38/genome/annotation/snps_indels/GCA_000001635.6_current_ids.vcf.gz'
  • Comment: The dbSNP database contains known single nucleotide polymorphisms, and is used in the annotation of known variants. Points to human dbSNP when * --gen_org human. JAX users should not change this parameter.

• --gen_ver

  • Default: "GRCm38.99"
  • Comment: snpEff genome version. Sets to 'hg38' when * --gen_org human

• --snpEff_config

  • Default: "/projects/omics_share/mouse/GRCm38/genome/indices/snpEff_5_1/snpEff.config"
  • Comment: The configuration file used while running snpEff, points to human snpEff file when * --gen_org human. JAX users should not change this parameter.

• --gold_std_indels

  • Default: '/projects/omics_share/human/GRCh38/genome/annotation/snps_indels/Mills_and_1000G_gold_standard.indels.hg38.vcf.gz’
  • Comment: Human Only - Used in GATK BaseRecalibrator. JAX users should not change this parameter.

• --phase1_1000G

  • Default: '/projects/omics_share/human/GRCh38/genome/annotation/snps_indels/1000G_phase1.snps.high_confidence.hg38.vcf.gz'
  • Comment: Human Only - Used in GATK BaseRecalibrator. JAX users should not change this parameter.

• --dbNSFP

  • Default: '/projects/omics_share/mouse/GRCm38/genome/annotation/snps_indels/GCA_000001635.6_current_ids.vcf.gz'
  • Comment: Used in variant annotation.

• --dbSNP_index

  • Default: '/projects/omics_share/mouse/GRCm38/genome/annotation/snps_indels/GCA_000001635.6_current_ids.vcf.gz.tbi'
  • Comment: Index file for dbSNP

• --cosmic

  • Default: '/projects/omics_share/human/GRCh38/genome/annotation/function/COSMICv95_Coding_Noncoding.gatk_formatted.vcf'
  • Comment: Human Only - Used in variant annotation.

• --cosmic_index

  • Default: '/projects/omics_share/human/GRCh38/genome/annotation/function/COSMICv95_Coding_Noncoding.gatk_formatted.vcf.gz.tbi'
  • Comment: Human Only - Index file for cosmic annotation file.

Pipeline Default Outputs

NOTE: * Represents a wild card that is a placeholder for values that will be filled by input file names and/or parameters when the pipeline is run.

Naming Convention	Description
wes_report.html	Nextflow autogenerated report
trace.txt	Nextflow trace of processes
*/bam/*.bam	Final BAM alignment file
*/bam/*.bai	Final BAM index file
*/*_variants_raw.vcf	VCF from GATK_HAPLOTYPECALLER
*/*_variants_raw.gvcf	GVCF output from GATK_HAPLOTYPECALLER. When run_gvcf = true.
*/*_INDEL_filtered_dbsnpID.vcf	Filtered unannotated INDELs only
*/*_SNP_filtered_dbsnpID.vcf	Filtered unannotated SNPs only
*/*SNP_INDEL_filtered_unannotated_final.vcf	Final VCF file, with filtered unannotated INDEL and SNP calls
*/*SNP_INDEL_filtered_annotated_final.vcf	Final VCF file, with filtered SNPeff annotated INDEL and SNP calls. See SNPEff notes below
*/*_snpsift_finalTable.txt	Extracted fields from final VCF, in tabular format. From SNPSIFT_EXTRACTFIELDS
*/stats/*fastp_report.html	FastP trimming metrics
*/stats/*_dup_metrics.txt	Duplicate metrics from PICARD_MARKDUPLICATES
*/stats/*_recal_data.table	(human only) recalibrated table for BQSR from GATK_BASERECALIBRATOR
*/stats/*_CoverageMetrics.txt	Coverage metrics from PICARD_COLLECTHSMETRICS

NOTE: In the final VCF file *SNP_INDEL_filtered_unannotated_final.vcf, the number of variants, will not match un-annotated variant counts (e.g., *SNP_INDEL_filtered_annotated_final.vcf). This difference in variant count is a function of SNPeff annotation.

From the SNPeff documentation:

Counting variants / annotations

It is important to remember that the VCF format specification allows having multiple variants in a single line. Also, a single variant can have more than one annotation, due to:

  * Multiple transcripts (isoforms) of a gene (e.g. the human genome has on average 8.8 transcripts per gene)
  * Multiple (overlapping) genes in the genomic location of the variant.
  * A variant spanning multiple genes (e.g. a translocation, large deletion, etc.)

When you count the number of variants, you must keep all these in mind to count them properly. Obviously, SnpEff does take all this into account when counting the variants for the summary HTML.

Typical counting mistake

Many people who claim that there is a mismatch between the number of variants in the summary (HTML) file and the number of variants in the VCF file, are just making mistakes when counting the variants because they forget one or more of these previous items.

A typical scenario is, for example, that people are "counting missense variants" using something like this:

grep missense file.vcf | wc -l

This is counting "lines in a VCF file that have at least one missense variants", as opposed to counting "missense annotations" and, as mentioned previously, the number of lines in a VCF file is not the same as the number of annotations or the number of variants.

Pipeline Optional Outputs

These output will only be saved when --keep_intermediate true is specified.

Naming Convention	Description (--keep_intermediate true)
*/*_read_group.txt	Read groups for fastq files from READ_GROUPS
*/*.sam	Sam alignment file from BWA_MEM
*/bam/*_sortsam.bam	Sorted bam from PICARD_SORTSAM
*/bam/*_realigned_BQSR.bam	(human only) output from GATK_APPLYBQSR
*/*_selectedvariants_INDEL.vcf	(human only) selected variants from GATK_SELECTVARIANTS
*/*_selectedvariants_SNP.vcf	(human only) selected variants from GATK_SELECTVARIANTS
*/*_snpeff_INDEL.vcf	(human only) from SNPEFF_HUMAN
*/*_snpeff_SNP.vcf	(human only) from SNPEFF_HUMAN

CSV Input Sample Sheet

The required input header is: sampleID,lane,fastq_1,fastq_2. Samples can be provided either paired or un-paired.

• The sampleID column is a unique identifies for each individual sample, which is associated with other samples based on status and patient ID.
• The lane column contains lane information for individual samples. If a single sample ID is provided with multiple lanes, the sequences from each lane will be concatenated prior to analysis.
• The fastq_1 and fastq_2 columns must contain absolute paths or URLs to read 1 and read 2 from an Illumina paired-end sequencing run.

Basic examples:

An example PE csv file:

sampleID,lane,fastq_1,fastq_2
Sample_42,Lane_1,/path/to/sample_42_001_R1.fastq.gz,/path/to/sample_42_001_R2.fastq.gz
Sample_42,Lane_2,/path/to/sample_42_002_R1.fastq.gz,/path/to/sample_42_002_R2.fastq.gz
Sample_101,Lane_1,/path/to/sample_101_001_R1.fastq.gz,/path/to/sample_101_001_R2.fastq.gz
Sample_10191,Lane_1,/path/to/sample_10191_001_R1.fastq.gz,/path/to/sample_10191_001_R2.fastq.gz

An example SE csv file:

sampleID,lane,fastq_1,fastq_2
Sample_42,Lane_1,/path/to/sample_42_001_R1.fastq.gz
Sample_42,Lane_2,/path/to/sample_42_002_R1.fastq.gz
Sample_101,Lane_1,/path/to/sample_101_001_R1.fastq.gz
Sample_10191,Lane_1,/path/to/sample_10191_001_R1.fastq.gz