A Dockstore WDL Workflow that uses UShER to place your consensus.fasta sequences on the SARS-CoV-2 phylogeny
UShER (Ultrafast Sample Placement on Existing Trees) uses maximum parsimony to rapidly place new samples onto an existing phylogenetic tree (1). We built this UShER workflow so that Public Health Labs with SARS-CoV-2 consensus fasta sequences can more easily identify the most similar viral genomes available in public databases. The output files could potentially be used to facilitate not only outbreak investigations, but also genomic surveillance. By modifying the input files, the same approach could be easily applied to monitor other pathogens.
You must have already aligned the fastq file reads to create a bam file for each sample, and then used that bam file to create a so-called "consensus fasta" file for each sample (e.g. using IVAR, freebayes, etc.). Concatenate these fasta files into a single fasta file, and make sure that each sequence has a unique fasta header. When this Workflow is running on Terra the text strings in the fasta file header rows is converted into Sample IDs in the output tables and data structures.
There are four suggested ways to run the publicTreeSamplePlacement workflow. In order of increasing complexity:
Depending on the number of samples (sequences) in your input fasta file a run will sometimes fail with an "Out of Memory" error. In our experience if this happens it is usually during the extractSubtrees task. In the runtime block of this task, one of the parameters, named 'maxRetries' has been set to '5'. When you launch the analysis from the Workflow tab in Terra, you have the option of selecting the option to "Retry with more memory", and if you click in this checkbox then Terra will attempt to rerun the workflow using more RAM.
In the current configuration, which uses Terra's filepath specification to provide the input parameters, all of your output data files from each WDL task should be located in Terra's Execution Directories (as shown in this example):
The new phylogenetic tree, named 'user_seqs.pb' is located in the Execution Directory of the usher task. By itself this file, which contains all of your sample sequences placed on the Global SARS-CoV-2 phylogeny, may not seem very useful (the phylogenetic tree it stores is extremely large with over a million sequences). However, one of the other WDL tasks named 'taxodium' uses the matUtils tools to reformat the user_seqs.pb file into one named 'user_seqs.taxodium.pb', which is in the taxodium Execution Directory. After downloading and saving this file to your localhost you can upload it into the taxonium visualization website located at https://cov2tree.org/. Just click on the 'Choose File' button and it will open a file selector widget and let you upload the new tree containing your sequences placed on the public BigTree. Taxonium lets you zoom in and zoom out of of locations in the tree, in real time (conceptually similar to the way that Google Maps lets you rapidly zoom in and out).
The UShER output files that often get used the most for visualization are located in the extractSubtrees Execution Directory on Terra. There is a TSV table, named subtree-assignments.tsv, produced by the matUtils tool, which lists all of your samples and the name of each sample's subtree (containing the 500 nearest neighbours in the larger tree), as well as related metadata. These subtree.json files are located in a subdirectory whose name begins 'glob-xxxxxx' (as shown in this screenshot from Terra):
Inside the glob-xxxx directory will be a collection of JSON documents that contain the "sub-trees" with your sample sequences, that have been extracted by matUtils from the user_seqs.pb tree:
After downloading those JSON renditions of the subtrees extracted from the huge user_seqs.pb tree, you can upload individual examples to Nextstrain's Auspice Tree Viewer, located here: https://auspice.us/ In Auspice there are many different options that you can set that lets you view and label the subtree using different facets, as shown in this example:
This is the only required file that the user must supply to run the publicTreeSamplePlacement Workflow. Each sequence must have a unique fasta header preceding it, where the first character on the header line is '>'. The file can contain from one SARS-CoV-2 consensus fasta sequence, up to thousands of sequences.
These five files are required by UShER and the other tools for the various steps. One of them is static, two of them are more stable, and two of them change every night.
The maftt multiple sequence aligner is used at this step. It is relatively rapid for such short genomes. This step makes sure that all of the sequences are the same length as the SARS-CoV-2 RefSeq, in preparation for the next step.
This step uses a program named faToVcf to convert all of your aligned, input fasta sequences into a much more compact, text-based data structure, using the common VCF format
UShER ingests both the VCF file with all of the annotated variants in your sample sequences, and the SARS-CoV-2 phylogenetic tree built using sequences available in public databases. The resulting user_seqs.pb file contains the information for a Newick tree, that has been converted into Google's protobuf file format.
matUtils integrates annotated genomic features in a GTF file from NCBI, and adds the predicted amino acid changes resulting from the nucleotide sequence variants found in your samples.