Harvest¶

Quickstart¶

parsnp -g <reference_replicon1,reference_replicon2,..> -d <genome_dir> -p <threads>

parsnp -r <reference_genome> -d <genome_dir> -p <threads>

parsnp -q <draft_assembly> -d <genome_db> -p <threads>

-c = <flag>: (c)urated genome directory, use all genomes in dir and ignore MUMi? (default = NO)
-d = <path>: (d)irectory containing genomes/contigs/scaffolds
-r = <path>: (r)eference genome (set to ! to pick random one from genome dir)
-g = <string>: Gen(b)ank file(s) (gbk), comma separated list (default = None)
-o = <string>: output directory? default [./P_CURRDATE_CURRTIME]
-q = <path>: (optional) specify (assembled) query genome to use, in addition to genomes found in genome dir (default = NONE)

-U = <float>: max MUMi distance value for MUMi distribution
-M = <flag>: calculate MUMi and exit? overrides all other choices! (default: NO)
-i = <float>: max MUM(i) distance (default: autocutoff based on distribution of MUMi values)

-a = <int>: min (a)NCHOR length (default = 1.1*Log(S))
-C = <int>: maximal cluster D value? (default=100)
-z = <path>: min LCB si(z)e? (default = 25)

-D = <float>: maximal diagonal difference? Either percentage (e.g. 0.2) or bp (e.g. 100bp) (default = 0.12)
-e = <flag> greedily extend LCBs? experimental! (default = NO)
-n = <string>: alignment program (default: libMUSCLE)
-u = <flag>: output unaligned regions? .unaligned (default: NO)

-x = <flag>: enable filtering of SNPs located in PhiPack identified regions of recombination? (default: NO)

-h = <flag>: (h)elp: print this message and exit
-p = <int>: number of threads to use? (default= 1)
-P = <int>: max partition size? limits memory usage (default= 15000000)
-v = <flag>: (v)erbose output? (default = NO)
-V = <flag>: output (V)ersion and exit

Installation from source¶

http://github.com/marbl/harvest

git clone https://github.com/marbl/parsnp.git parsnp_src
cd parsnp_src

cd muscle
./autogen.sh
./configure --prefix=`pwd` CXXFLAGS=’-fopenmp’
make install

cd ..
./autogen.sh
./configure
make install

export PARSNPDIR=/path/to/parsnp/install

cd muscle
./autogen.sh
./configure --prefix=/usr/local/
sudo make install

cd ..
./autogen.sh
./configure --prefix=/usr/local/
sudo make install

Tutorial¶

Parsnp download instructions

To further demonstrate the functionality of Parsnp we have prepared two small tutorial datasets. The first dataset is a MERS coronavirus outbreak dataset involving 49 isolates. The second dataset is a selected set of 31 Streptococcus pneumoniae genomes. Both of these datasets should run on modestly equipped laptops in a few minutes.

1. 49 MERS Coronavirus genomes

   • Download genomes:

     • gzipped tarball

   • Run parsnp with default parameters

     ./parsnp -g ./ref/EMC_2012.gbk -d ./mers49 -c
     

   • Command-line output

   
   • Visualize with Gingr GGR
   

   • Configure parameters

     • 95% of the reference is covered by the alignment. This is <100% mainly due to a 1kbp unaligned region from 26kbp to 27kbp.

     • To force alignment across large collinear regions, use the -C maximum distance between two collinear MUMs:

       ./parsnp -g ./ref/EMC_2012.gbk -d ./mers49 -C 1000 -c
       

   • Visualize again with Gingr GGR

     • By adjusting the -C parameter, this region is no longer unaligned, boosting the reference coverage to 97%.

     

   • Zoom in with Gingr for nucleotide view of region

     • On closer inspection, a large stretch of N’s in Jeddah isolate C7569 was the culprit

     

   • Inspect Output:

     • Multiple alignment: XMFA
     • SNPs: VCF
     • Phylogeny: Newick

2. 31 Streptococcus pneumoniae genomes

   • Download genomes:

     • gzipped tarball

   • Run parsnp

     ./parsnp -r ./strep31/NC_011900.fna -d ./strep31 -p <num threads>
     

   • Command-line output:

     

   • Force inclusion of all genomes (-c)

     ./parsnp -r ./strep31/NC_011900.fna -d ./strep31 -p <num threads> -c
     

• Command-line output:

  

• Visualize with Gingr GGR

  

• Enable recombination detection/filter (-x)

  ./parsnp -r ./strep31/NC_011900.fna -d ./strep31 -p <num threads> -c -x
  

• Re-visualize with Gingr GGR

  • Bootstrap values have improved after running recombination filter; columns with filtered SNPs are displayed in image:

  

• Inspect Output:

  • Multiple alignment: XMFA
  • SNPs: VCF
  • Phylogeny: Newick

FAQs¶

17. What is Parsnp ?

1. Parsnp takes both draft and finished genomes of closely related strains as input, performs conservative core genome alignment and as output returns multi-alignments (XMFA), variants (VCF), core genome phylogeny (Newick) and Gingr input format (GGR).

17. Why should I use Parsnp?

1. The main advantages of Parsnp over alternative approaches is robust filtration of variant (SNP) calls, multiple alignments as output and superior speed. Parsnp can align 200-300 bacterial strains in <30 minutes on a 16-core server and ~1000 in a couple of hours.

17. Why should I *not* use Parsnp?

1. If you are interested in pan genome/whole genome alignment, existing tools for the job that perform well include Mauve, Mugsy, among others. In addition, Parsnp is tailored for intraspecific genome analysis (outbreak analysis of a pathogen, etc). One main limitation of Parsnp is that it cannot handle subsets (core genome only) and is not as sensitive as existing methods.

17. How can I visualize the results?

1. Gingr (http://github.com/marbl/gingr) can open Parsnp output and provide an interactive display of multi-alignments, variants and the phylogenetic tree estimated from the core genome alignment.

17. What % of genome X is aligned? What is the core genome alignment size?

1. Within the log output, there are coverage values listed that individually indicate the percentage of a given genome that is included in the core genome alignment. Note, this includes the Muscle aligned-regions plus the maximal unique matches (MUMs). The core genome alignment size can then be calculated by multiplying the coverage value, for a given genome, by its length.

17. Only a small percentage (<40%) of the reference genome covered by the alignments, huh?

1. Parsnp is a conservative core genome alignment method that necessarily requires that all genomes are present in each aligned regions. The focus is on aligning 1000s of closely related bacterial strains quickly while maintaining sensitivity comparable to existing WGA methods. In additon, the core genome has been shown to contain as few as 30-40% of the gene content (even in very closely-related clades) due to reductive genome evolution and/or a large accessory genome (with plenty of IS/phage elements). However, for increased sensitivity w.r.t aligned regions, and alignments containing subsets, both Mugsy and Mauve are terrific tools for the job.

17. I am not sure which reference genome to choose, help!

1. Since the core necessarily includes all genomes, the choice of reference does not matter (with a couple important exceptions, continue reading). Feel free to use the parameter ‘-r !’ to randomly select a reference if you are feeling indecisive or if all of the genomes contained in the genome directory are of similar quality. Typically, finished/closed genomes are used a the reference strain to ensure they are high-quality and do not contain assembly artifacts or contaminant.

17. How are the genomes selected from the input directory? Not all of the genomes I wanted included are in the tree!

1. By default, parsnp calculates the MUMi distance between the reference and each of the genomes in the genome directory. All genomes with MUMi distance <= 0.01 are included, all others are discarded. To force all genomes present in the genome dir to be included simply include ‘-c’ as a command-line parameter.

17. Parsnp fails to report a SNP position output by method X, why is this?

1. The goal of parsnp is to capture all informative signals found in the core genome of the specified clade of interest. Any SNPs in regions not shared by all genomes will not reported. Additionally, any SNP found in a likely poorly aligned region would also be discarded. Finally, parsnp does not perform LCB extension and therefore may miss SNPs appearing at the end of locally conserved blocks or clusters.

Requirements¶

Tutorial¶

Running Gingr¶

Mac OS X¶

• Gingr.app can be moved to the Applications folder if desired
• Double-click Gingr.app to run
• Depending on your security settings, there may be an error that Gingr is not from the Mac App Store or is from an unidentified developer. To run it anyway:

• Right click on Gingr.app
• Select “Open” from the menu
• Click the “Open” button at the next prompt

Linux¶

• From the desktop

• Click on the “gingr” binary

• From a terminal

• Navigate to the folder with the “gingr” binary
• Run ”./gingr”

Browsing a Gingr file¶

• Download Gingr input file
• Open in Gingr (File->Open)

• The phylogeny appears on the left. Hover over a clade to highlight and outline the corresponding tracks to the right.

• Click to zoom in on the clade

• The multiple alignment appears on the right, shown as a SNP heatmap when zoomed out. To see the full alignment, zoom in with the mouse wheel or by selecting a region in the ruler.

Importing other files¶

• Create a new workspace (File->New)

• Download the data files
  • Alignment: xmfa
  • Reference: fasta
  • Annotations: genbank
  • Phylogeny: newick
• Open the XMFA alignment (File->Open). Since XMFA files can be accompanied by reference files, the Open dialog will appear. Choose the Fasta file as the reference in this window.

• The preview panes allow you to ensure that the header for the reference is the same as the first sequence in the XMFA. This allows sequences between LCBs to be shown and allows annotations to be added later.

• The track highlighted in blue (“england.gbk.fna.srt”) is the current reference for variants. Select a new reference by right-clicking on a track.

• Next, import the phylogenetic tree (File->Open)

• Reroot the tree at the midpoint (Tree->Reroot at midpoint)

• The tree will now be balanced at the center of the longest path

• Finally, import the annotations (File->Open)

• The workspace can be saved to share or return to later (File->Save)

File formats¶

The flowchart below describes the various file formats that can be imported or exported to/from Gingr (or the harvesttools command line utility).

• Alignments

  • Core only: The Gingr file format stores core alignments, or alignments that involve all genomes. When alignments are loaded, blocks that are not core will be discarded.
  • When importing MAF alignments, the first sequence of the first core block is used as the reference. Each reference contig will be padded with Ns up to its first LCB and between subsequent LCBs. If alignment blocks overlap in reference coordinate space, the block seen earlier in the file will be kept; the later one will be ignored.
  • Multi-fasta: This format does not store rearrangement information, so the alignment is treated as a single LCB.

• References

  • XMFA files can be accompanied by Fasta reference files to provide sequence between LCBs and to allow Genbank annotations (which must have matching GI numbers) to be loaded later. Genbank files that contain sequence can also be used as references.
  • Multi-fasta alignments will use the first sequence as the reference. Genbank annotations can be loaded later if the GIs match.

• Variants

  • VCF files must be imported with a Fasta reference. The only fields imported are sequence identifier (CHROM), position (POS), reference allele (REF), alternate alleles (ALT), quality (QUAL), filters (FILTER, including ##FILTER specifications in the header), and genotype (GT); all other information is ignored. Additionally, Since VCF does not store complete alignment information, any insertions larger than one base will be replaced by an LCB boundary when importing. If a genotype is diploid or polyploid, only the first haplotype is used in the multi-alignment (the others are ignored). Symbolic alleles and breakends are currently unsupported and will also be ignored. When writing to VCF, only the imported fields will be populated. Indel output is also currently unimplemented, so indels will be skipped when writing.
  • The multi-fasta SNP output is the same format as multi-fasta alignments, but only contains columns with unfiltered (“PASS”) variants (like a Mauve SNP file). This is useful for generating phylogenetic trees, but does not contain positional information or rearrangements.

Quickstart¶

harvesttools -g <reference_genbank_file1> -r <reference fasta file> -x <XMFA file> -o hvt.ggr

harvesttools -i input.ggr -X output.xmfa

harvesttools -i input.ggr -S output.snps

Installation from source¶

http://github.com/marbl/harvest

 git clone https://github.com/marbl/harvest-tools.git hvt_src
 cd hvt_src

./bootsrap.sh
./configure [--prefix=...] [--with-protobuf=...]
make
[sudo] make install