search

TRON-Bioinformatics / covigator-ngs-pipeline

A Nextflow pipeline for NGS variant calling on SARS-CoV-2. From FASTQ files to normalized and annotated VCF files from GATK, BCFtools, LoFreq and iVar.
MIT License
17 stars 7 forks source link
bioinformatics coronavirus covid-19 genomics nextflow ngs sars-cov-2 variant-callers

readme

CoVigator logo

CoVigator pipeline: variant detection pipeline for Sars-CoV-2

Nextflow run with conda DOI Run tests License

The Covigator pipeline processes SARS-CoV-2 FASTQ or FASTA files into annotated and normalized analysis ready VCF files. It also classifies samples into lineages using pangolin. The pipeline is implemented in the Nextflow framework (Di Tommaso, 2017), it is a stand-alone pipeline that can be used independently of the CoVigator dashboard and knowledge base.

Although it is configured by default for SARS-CoV-2 it can be employed for the analysis of other microbial organisms if the required references are provided.

The result of the pipeline is one or more annotated VCFs with the list of SNVs and indels ready for analysis.

The results from the CoVigator pipeline populate our CoVigator dashboard https://covigator.tron-mainz.de

Table of Contents

  1. Two pipelines in one
  2. Implementation
  3. How to run
  4. Understanding the output
  5. Annotation resources
  6. Future work
  7. Bibliography

Two pipelines in one

In CoVigator we analyse samples from two different formats, FASTQ files (e.g.: as provided by the European Nucleotide Archive) and FASTA files containing a consensus assembly. While from the first we get the raw reads, from the second we obtain already assembled genomes. Each of these formats has to be analysed differently. Also, the output data that we can obtain from each of these is different.

CoVigator pipeline

Pipeline for FASTQ files

When FASTQ files are provided the pipeline includes the following steps:

Both single end and paired end FASTQ files are supported.

Pipeline for FASTA files

When a FASTA file is provided with a single assembly sequence the pipeline includes the following steps:

The FASTA file is expected to contain a single assembly sequence. Bear in mind that only clonal variants can be called on the assembly.

Pipeline for VCF files

When a VCF file is provided the pipeline includes the following steps:

Implementation

The pipeline is implemented as a Nextflow workflow with its DSL2 syntax. The dependencies are managed through a conda environment to ensure version traceability and reproducibility. The references for SARS-CoV-2 are embedded in the pipeline. The pipeline is based on a number of third-party tools, plus a custom implementation based on biopython (Cock, 2009) for the alignment and subsequent variant calling over a FASTA file.

All code is open sourced in GitHub https://github.com/TRON-Bioinformatics/covigator-ngs-pipeline and made available under the MIT license. We welcome any contribution. If you have troubles using the CoVigator pipeline or you find an issue, we will be thankful if you would report a ticket in GitHub.

The alignment, BAM preprocessing and variant normalization pipelines are based on the implementations in additional Nextflow pipelines within the TronFlow initiative https://tronflow-docs.readthedocs.io/.

Variant annotations

The variants derived from a FASTQ file are annotated on the FILTER column using the VAFator (https://github.com/TRON-Bioinformatics/vafator) variant allele frequency (VAF) into LOW_FREQUENCY, SUBCLONAL, LOW_QUALITY_CLONAL and finally PASS variants correspond to clonal variants. By default, variants with a VAF < 2 % are considered LOW_FREQUENCY, variants with a VAF >= 2 % and < 50 % are considered SUBCLONAL and variants with a VAF >= 50 % and < 80 % are considered LOW_QUALITY_CLONAL. This thresholds can be changed with the parameters --low_frequency_variant_threshold, --subclonal_variant_threshold and --low_quality_clonal_variant_threshold respectively.

VAFator technical annotations:

SnpEff provides the functional annotations. And all mutations are additionally annotated with the following SARS-CoV-2 specific annotations:

Biological annotations:

According to DeMaio et al. (2020), mutations at the beginning (ie: POS <= 50) and end (ie: POS >= 29,804) of the genome are filtered out

This is an example of biological annotations of a missense mutation in the spike protein on the N-terminal subunit 1 domain.

ANN=A|missense_variant|MODERATE|S|gene-GU280_gp02|transcript|TRANSCRIPT_gene-GU280_gp02|protein_coding|1/1|c.118G>A|
p.D40N|118/3822|118/3822|40/1273||;CONS_HMM_SARS_COV_2=0.57215;CONS_HMM_SARBECOVIRUS=0.57215;CONS_HMM_VERTEBRATE_COV=0;
PFAM_NAME=bCoV_S1_N;PFAM_DESCRIPTION=Betacoronavirus-like spike glycoprotein S1, N-terminal

Phasing limitations

The phasing implementation is applicable only to clonal mutations. It assumes all clonal mutations are in phase and hence it merges those occurring in the same amino acid. In order to phase intrahost mutations we would need to implement a read-backed phasing approach such as in WhatsHap or GATK's ReadBackedPhasing. Unfortunately these tools do not support the scenario of a haploid organism with an undefined number of subclones. For this reason, phasing is implemented with custom Python code at bin/phasing.py.

Primers trimming

With some library preparation protocols such as ARTIC it is recommended to trim the primers from the reads. We have observed that if primers are not trimmed spurious mutations are being called specially SNVs with lower frequencies and long deletions. Also the variant allele frequencies of clonal mutations are underestimated.

The BED files containing the primers for each ARTIC version can be found at https://github.com/artic-network/artic-ncov2019/tree/master/primer_schemes/nCoV-2019.

If the adequate BED file is provided to the CoVigator pipeline with --primers the primers will be trimmed with iVar. This affects the output of every variant caller, not only iVar.

Reference data

The default SARS-CoV-2 reference files correspond to Sars_cov_2.ASM985889v3 and were downloaded from Ensembl servers. No additional parameter needs to be provided to use the default SARS-CoV-2 reference genome.

Using a custom reference genome

These references can be customised to use a different SARS-CoV-2 reference or to analyse a different virus.
Two files need to be provided:

As well as an organism name for the SnpEff annotation:

NOTE: beware that for Nextflow to find these files the reference needs to be passed as an absolute path.

When a custom reference genome is provided, the pipeline automatically generates a bwa-mem2 index in the case of fastq input and a custom SnpEff database for functional annotations using the steps described below. The steps are listed here for documentation purposes and should not be performed manually. If the indexes should be reused for multiple runs, they can be generated permanently with the --reference_generate option. Afterwards the generated references can either be set in the nextflow.config file or set via command line arguments of the workflow.

The reference genome FASTA needs to be indexed for different components of the pipeline. bwa-mem2 indices, .fai index and a .dict index can be generated with the following three commands:

bwa-mem2 index reference.fasta
samtools faidx reference.fasta
gatk CreateSequenceDictionary --REFERENCE your.fasta

In order to have SnpEff functional annotations available you will need to prepare the new reference with SnpEff.

When running CoVigator you will also need to provide three parameters:

Limitations

Intrahost mutations

Some mutations may be observed in a subset of the virus sample, this may arise through intrahost virus evolution or co-infection. Intrahost mutations can only be detected when analysing the raw reads (ie: the FASTQs) as in the assembly (ie: the FASTA file) a single virus consensus sequence is represented. BCFtools and GATK do not normally capture intrahost mutations; on the other hand LoFreq and iVar both capture mutations that deviate from a clonal-like VAF. Nevertheless, mutations with lower variant allele frequency (VAF) are challenging to distinguish from sequencing and analytical errors.

Mutations are annotated on the FILTER column using the VAF into three categories:

Other low quality mutations are removed from the output.

The VAF thresholds can be changed with the parameters --low_frequency_variant_threshold, --subclonal_variant_threshold and --low_quality_clonal_variant_threshold.

Lineage only mode

At some point, it might be desirable to run only the lineage determination to obtain the pangolin calls of a sample. To solve this problem, we have implemented the lineage-only mode, which can be executed by specifying the option "--lineage_mode" option. Please note that this mode can only be used with FASTA or VCF as input.

How to run

Requirements

Testing

To run the workflow on a test assembly dataset run:

nextflow run tron-bioinformatics/covigator-ngs-pipeline -profile conda,test_fasta

Find the output in the folder covigator_test_fasta.

To run the workflow on a test raw reads dataset run:

nextflow run tron-bioinformatics/covigator-ngs-pipeline -profile conda,test_fastq

Find the output in the folder covigator_test_fastq.

The above commands are useful to create the conda environments beforehand.

NOTE: pangolin is the most time-consuming step of the whole pipeline. To make it faster, locate the conda environment that Nextflow created with pangolin (eg: find $YOUR_NEXTFOW_CONDA_ENVS_FOLDER -name pangolin) and run pangolin --decompress-model.

Running

For paired end reads:

nextflow run tron-bioinformatics/covigator-ngs-pipeline \
[-r v0.10.0] \
[-profile conda] \
--fastq1 <FASTQ_FILE> \
--fastq2 <FASTQ_FILE> \
--name example_run \
--output <OUTPUT_FOLDER> \
[--reference <path_to_reference>/Sars_cov_2.ASM985889v3.fa] \
[--gff <path_to_reference>/Sars_cov_2.ASM985889v3.gff3]

For single end reads:

nextflow run tron-bioinformatics/covigator-ngs-pipeline \
[-r v0.10.0] \
[-profile conda] \
--fastq1 <FASTQ_FILE> \
--name example_run \
--output <OUTPUT_FOLDER> \
[--reference <path_to_reference>/Sars_cov_2.ASM985889v3.fa] \
[--gff <path_to_reference>/Sars_cov_2.ASM985889v3.gff3]

For assembly:

nextflow run tron-bioinformatics/covigator-ngs-pipeline \
[-r v0.10.0] \
[-profile conda] \
--fasta <FASTA_FILE> \
--name example_run \
--output <OUTPUT_FOLDER> \
[--reference <path_to_reference>/Sars_cov_2.ASM985889v3.fa] \
[--gff <path_to_reference>/Sars_cov_2.ASM985889v3.gff3]

For VCF:

nextflow run tron-bioinformatics/covigator-ngs-pipeline \
[-r v0.10.0] \
[-profile conda] \
--vcf <VCF_FILE> \
--name example_run \
--output <OUTPUT_FOLDER> \
[--reference <path_to_reference>/Sars_cov_2.ASM985889v3.fa] \
[--gff <path_to_reference>/Sars_cov_2.ASM985889v3.gff3]

As an optional input when processing directly VCF files you can provide BAM files to annotate VAFs:

nextflow run tron-bioinformatics/covigator-ngs-pipeline \
[-r v0.10.0] \
[-profile conda] \
--vcf <VCF_FILE> \
--bam <BAM_FILE> \
--bai <BAI_FILE> \
--name example_run \
--output <OUTPUT_FOLDER> \
[--reference <path_to_reference>/Sars_cov_2.ASM985889v3.fa] \
[--gff <path_to_reference>/Sars_cov_2.ASM985889v3.gff3]

For batch processing of reads use --input_fastqs_list and --name.

nextflow run tron-bioinformatics/covigator-ngs-pipeline [-profile conda] --input_fastqs_list <TSV_FILE> --library <paired|single> --output <OUTPUT_FOLDER> [--reference <path_to_reference>/Sars_cov_2.ASM985889v3.fa] [--gff <path_to_reference>/Sars_cov_2.ASM985889v3.gff3]

where the TSV file contains two or three columns tab-separated columns without header. Columns: sample name, path to FASTQ 1 and optionally path to FASTQ 2.

Sample FASTQ 1 FASTQ 2 (optional column)
sample1 /path/to/sample1_fastq1.fastq /path/to/sample1_fastq2.fastq
sample2 /path/to/sample2_fastq1.fastq /path/to/sample2_fastq2.fastq
... ... ...

For batch processing of assemblies use --input_fastas_list.

nextflow run tron-bioinformatics/covigator-ngs-pipeline [-profile conda] --input_fastas_list <TSV_FILE> --library <paired|single> --output <OUTPUT_FOLDER> [--reference <path_to_reference>/Sars_cov_2.ASM985889v3.fa] [--gff <path_to_reference>/Sars_cov_2.ASM985889v3.gff3]

where the TSV file contains two columns tab-separated columns without header. Columns: sample name and path to FASTA.

Sample FASTA
sample1 /path/to/sample1.fasta
sample2 /path/to/sample2.fasta
... ...

For batch processing of VCFs use --input_vcfs_list.

nextflow run tron-bioinformatics/covigator-ngs-pipeline [-profile conda] --input_vcfs_list <TSV_FILE> --output <OUTPUT_FOLDER> [--reference <path_to_reference>/Sars_cov_2.ASM985889v3.fa] [--gff <path_to_reference>/Sars_cov_2.ASM985889v3.gff3]

where the TSV file contains two columns tab-separated columns without header. Columns: sample name and path to VCF.

Sample FASTA
sample1 /path/to/sample1.vcf
sample2 /path/to/sample2.vcf
... ...

Optionally, provide BAM files for batch processing of VCFs using --input_bams_list.

nextflow run tron-bioinformatics/covigator-ngs-pipeline [-profile conda] \
  --input_vcfs_list <TSV_FILE> \
  --input_bams_list <TSV_FILE> \
  --output <OUTPUT_FOLDER> \
  [--reference <path_to_reference>/Sars_cov_2.ASM985889v3.fa] \
  [--gff <path_to_reference>/Sars_cov_2.ASM985889v3.gff3]

where the BAMs TSV file contains three columns tab-separated columns without header. Columns: sample name, path to BAM and path to BAI.

Sample BAM BAI
sample1 /path/to/sample1.bam /path/to/sample1.bai
sample2 /path/to/sample2.bam /path/to/sample2.bai
... ... ...

For lineage-only:

nextflow run tron-bioinformatics/covigator-ngs-pipeline \
[-r v0.10.0] \
[-profile conda] \
--fasta <FASTA_FILE> \
--name example_run \
--output <OUTPUT_FOLDER> \
--lineage_mode

For preparation of custom references:

nextflow run tron-bioinformatics/covigator-ngs-pipeline \
[-r v0.10.0] \
[-profile conda] \
--reference <path_to_reference>/<custom_genome>.fasta \
--gff <path_to_reference>/<custom_genome_annotation>.gff3
--output <OUTPUT_FOLDER> \
--snpeff_organism <name_of_organism>
--reference_generate

Getting help

You can always contact us directly or create a GitHub issue, otherwise see all available options using --help:

$ nextflow run tron-bioinformatics/covigator-ngs-pipeline -profile conda --help

Usage:
    nextflow run tron-bioinformatics/covigator-ngs-pipeline -profile conda --help

Input:
    * --fastq1: the first input FASTQ file (not compatible with --fasta, nor --vcf)
    * --fasta: the FASTA file containing the assembly sequence (not compatible with --fastq1, nor --vcf)
    * --vcf: the VCF file containing mutations to analyze (not compatible with --fastq1, nor --fasta)
    * --bam: the BAM file containing reads to annotate VAFs on a VCF (not compatible with --fastq1, nor --fasta)
    * --bai: the BAI index for a BAM file (not compatible with --fastq1, nor --fasta)
    * --name: the sample name, output files will be named after this name
    * --output: the folder where to publish output
    * --input_fastqs_list: alternative to --name and --fastq1 for batch processing
    * --library: required only when using --input_fastqs
    * --input_fastas_list: alternative to --name and --fasta for batch processing
    * --input_vcfs_list: alternative to --name and --vcf for batch processing
    * --input_bams_list: alternative to --name, --vcf, --bam and --bai for batch processing

Optional input only required to use a custom reference:
    * --reference: the reference genome FASTA file, *.fai, *.dict and bwa indexes are required.
    * --gff: the GFFv3 gene annotations file (required to run iVar and to phase mutations from all variant callers)    
    * --snpeff_data: path to the SnpEff data folder, it will be useful to use the pipeline on other virus than SARS-CoV-2
    * --snpeff_config: path to the SnpEff config file, it will be useful to use the pipeline on other virus than SARS-CoV-2
    * --snpeff_organism: organism to annotate with SnpEff, it will be useful to use the pipeline on other virus than SARS-CoV-2
    * --reference_generate: Run reference generate to prepare reusable custom reference

Optional input:
    * --fastq2: the second input FASTQ file
    * --primers: a BED file containing the primers used during library preparation. If provided primers are trimmed from the reads.
    * --min_base_quality: minimum base call quality to take a base into account for variant calling (default: 20)
    * --min_mapping_quality: minimum mapping quality to take a read into account for variant calling (default: 20)
    * --vafator_min_base_quality: minimum base call quality to take a base into account for VAF annotation (default: 0)
    * --vafator_min_mapping_quality: minimum mapping quality to take a read into account for VAF annotation (default: 0)
    * --low_frequency_variant_threshold: VAF threshold to mark a variant as low frequency (default: 0.02)
    * --subclonal_variant_threshold: VAF superior threshold to mark a variant as subclonal  (default: 0.5)
    * --lq_clonal_variant_threshold: VAF superior threshold to mark a variant as loq quality clonal (default: 0.8)
    * --memory: the ammount of memory used by each job (default: 3g)
    * --cpus: the number of CPUs used by each job (default: 1)
    * --skip_lofreq: skips calling variants with LoFreq
    * --skip_gatk: skips calling variants with GATK
    * --skip_bcftools: skips calling variants with BCFTools
    * --skip_ivar: skips calling variants with iVar
    * --skip_pangolin: skips lineage determination with pangolin
    * --match_score: global alignment match score, only applicable for assemblies (default: 2)
    * --mismatch_score: global alignment mismatch score, only applicable for assemblies (default: -1)
    * --open_gap_score: global alignment open gap score, only applicable for assemblies (default: -3)
    * --extend_gap_score: global alignment extend gap score, only applicable for assemblies (default: -0.1)
    * --skip_sarscov2_annotations: skip some of the SARS-CoV-2 specific annotations (default: false)
    * --keep_intermediate: keep intermediate files (ie: BAM files and intermediate VCF files)
    * --args_bcftools_mpileup: additional arguments for bcftools mpileup command (eg: --args_bcftools_mpileup='--ignore-overlaps')
    * --args_bcftools_call: additional arguments for bcftools call command (eg: --args_bcftools_call='--something')
    * --args_lofreq: additional arguments for lofreq command (eg: --args_lofreq='--something')
    * --args_gatk: additional arguments for gatk command (eg: --args_gatk='--something')
    * --args_ivar_samtools: additional arguments for ivar samtools mpileup command (eg: --args_ivar_samtools='--ignore-overlaps')
    * --args_ivar: additional arguments for ivar command (eg: --args_ivar='--something')

Output:
    * Output a VCF file for each of BCFtools, GATK, LoFreq and iVar when FASTQ files are
    provided or a single VCF obtained from a global alignment when a FASTA file is provided.
    * A pangolin results file for each of the VCF files.
    * Only when FASTQs are provided:
      * FASTP statistics
      * Depth and breadth of coverage analysis results

Understanding the output

Although the VCFs are normalized for both pipelines, the FASTQ pipeline runs four variant callers, while the FASTA pipeline runs a single variant caller. Also, there are several metrics in the FASTQ pipeline that are not present in the output of the FASTA pipeline. Here we will describe these outputs.

FASTQ pipeline output

Find in the table below a description of each of the expected files and a link to a sample file for the FASTQ pipeline. The VCF files will be described in more detail later.

Name Description Sample file
$NAME.fastp_stats.json Output metrics of the fastp trimming process in JSON format ERR4145453.fastp_stats.json
$NAME.fastp_stats.html Output metrics of the fastp trimming process in HTML format ERR4145453.fastp_stats.html
$NAME.deduplication_metrics.txt Deduplication metrics ERR4145453.deduplication_metrics.txt
$NAME.coverage.tsv Coverage metrics (eg: mean depth, % horizontal coverage) ERR4145453.coverage.tsv
$NAME.depth.tsv Depth of coverage per position ERR4145453.depth.tsv
$NAME.bcftools.vcf.gz Bgzipped, tabix-indexed and annotated output VCF from BCFtools ERR4145453.bcftools.normalized.annotated.vcf.gz
$NAME.gatk.vcf.gz Bgzipped, tabix-indexed and annotated output VCF from GATK ERR4145453.gatk.normalized.annotated.vcf.gz
$NAME.lofreq.vcf.gz Bgzipped, tabix-indexed and annotated output VCF from LoFreq ERR4145453.lofreq.normalized.annotated.vcf.gz
$NAME.ivar.vcf.gz Bgzipped, tabix-indexed and annotated output VCF from LoFreq ERR4145453.ivar.tsv
$NAME.lofreq.pangolin.csv Pangolin CSV output file derived from LoFreq mutations ERR4145453.lofreq.pangolin.csv

FASTA pipeline output

The FASTA pipeline returns a single VCF file. The VCF files will be described in more detail later.

Name Description Sample file
$NAME.assembly.vcf.gz Bgzipped, tabix-indexed and annotated output VCF ERR4145453.assembly.normalized.annotated.vcf.gz

Annotations resources

SARS-CoV-2 ASM985889v3 references were downloaded from Ensembl on 6th of October 2020:

ConsHMM mutation depletion scores downloaded on 1st of July 2021:

Gene annotations including Pfam domains downloaded from Ensembl on 25th of February 2021 from:

Future work

Bibliography