Skip to content

Pipeline Overview

Jill V. Hagey, PhD edited this page Dec 5, 2024 · 6 revisions

Pipeline Summary:

Workflow Diagram

PhyloPHoeNIx SNVPhyl
Image 1 Image 2

While this image is neat, lets break down what is actually happening in PhyloPhoenix:

Input Prep

  1. Phylophoenix takes in either a directory and makes a samplesheet from it or takes a samplesheet directly as input. This file is used to find the files necessary to run PhyloPhoenix, specifically the reads, Pass/Fail determination and MLST. It also makes a GRiPHin file from these samples. 2 (optional). If a metadata file is provided with --metadata it will undergo a series of checks to ensure the file and its data is in the format that will work for Microreact. Specifically, dates are checked and if locations (country, state, county and city) are provided this information is used to get a center latitude and longitude for the lowest hierarchical location provided. 3 (optional). If --by_st is passed then GRiPHin is used to split the samples by their MLST (Primary_MLST column in Griphin) and run them through the exact same process that a set of all the samples goes through.
  2. The samplesheet is rearranged to construct the path to the assemblies created previously by PHoeNIx and confirms they exist and are in the expected format.
  3. Scaffolds are used to create pairwise comparisons for all the assemblies in the set and the centriod is determined and written to <ST>_centroid_info.txt.
  4. The sample reads and determined reference are passed to SNVPhyl for SNV Matrix and Tree construction.
  • If there are >3 samples in a set then centroid picking and tree construction are skipped, but a SNV Matrix will still be created.

SNVPhyl Steps (adapted from SNVPhyl paper:

  1. Repeat Identification: Identifies internal repeat regions on the reference genome and generates a masking file containing the locations of repetitive regions to exclude from analysis. This file is concatenated to the user-supplied masking file, if defined, and used in later analysis stages.

  2. Mapping/Variant Calling: Aligns the reads to the determined reference genome. Reference mapping is performed using SMALT (version 0.7.5), which outputs a read pileup.

  3. Mapping Quality: Each pileup is evaluated for the mean coverage across a user-defined proportion of the reference genome (e.g. 10× coverage across at least 80 % of the genome). Any sequenced genomes that do not meet the minimum mean coverage threshold are flagged for further assessment.

  4. Variant calling stage: Uses two independent variant callers, FreeBayes (version 0.9.20), and the SAMtools and BCFtools packages. FreeBayes is run using the haploid variant calling mode and the resulting variants are filtered to remove insertions/deletions and split complex variant calls. SAMtools and BCFtools are run independently of FreeBayes and are used to confirm the FreeBayes variant calls and generate base calls for non-variant positions.

  5. Variant Consolidation: Combines both sets of variant and non-variant (polymorphic and monomorphic) calls into a merged file, flagging mismatches between variant callers. Base calls below the defined minimum read coverage are identified and flagged. The merged base calls are scanned for positions that do not pass the minimum relative SNV abundance (proportion of reads supporting the SNV with respect to the depth of coverage at a site) and minimum mean mapping quality. These base calls are removed from the merged base calls file. The remaining base calls that pass all these criteria are defined as either a high-quality SNV (hqSNV) or a high-quality non-variant base call. The hqSNVs are scanned to identify high-density SNV regions. These regions are identified by passing a sliding window of a given size along the genome and counting the number of SNVs within the window that exceed a given SNV density threshold. The high-density SNV regions are recorded in a tab-delimited file and used to mask potential recombinant regions.

  6. SNV Alignment Generation: Examines the merged base calls to generate a table of identified variants, and an alignment of hqSNVs and the corresponding high-quality non-variant bases. The hqSNVs are evaluated and assigned a status using the base calls at the same reference genome position for every isolate. A status of ‘valid’ is assigned when the base calls from all isolates in the same position pass the minimum criteria (hqSNVs or high-quality non-variants). These base calls are incorporated into the SNV alignment used for phylogeny generation. A status of ‘filtered-coverage’ is assigned when one or more isolates fail the minimum base coverage threshold at a particular position and the failed isolates’ base calls are annotated as ‘−’ (indicating no nucleotide or a gap). A status of ‘filtered-mpileup’ is assigned when one or more isolates have conflicting base calls between FreeBayes and SAMtools/BCFtools and the conflicting isolates’ base calls are annotated as ‘N’ (indicating any nucleotide non-specifically). A status of ‘filtered-invalid’ is assigned when the identified hqSNV overlaps one of the masked locations. The hqSNVs, base calls and assigned status are recorded in the SNV table and saved for later inspection. The SNV table can be used to re-generate the downstream SNV alignment and phylogenetic tree without re-running the computationally intensive reference mapping and variant calling steps.

Reference mapping and variant calling to identify high-quality SNVs (hqSNVs) using SMALT, FreeBayes and SAMtools/BCFtools. Merging the identified hqSNVs to construct a multiple sequence alignment. Building a maximum likelihood tree with PhyML and generating other output files.