diff --git a/docs/conf.py b/docs/conf.py index d1d2ffea3..db9ede381 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -70,8 +70,6 @@ def prose_list(items): 'faq/seq_traits.md', 'faq/translate_ref.md', 'faq/vcf_input.md', - 'tutorials/tb_tutorial.md', - 'tutorials/zika_tutorial.md', 'usage/augur_snakemake.md', ] diff --git a/docs/index.rst b/docs/index.rst index ead0a788c..43700b677 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -7,9 +7,9 @@ Augur: A bioinformatics toolkit for phylogenetic analysis .. note:: The documentation you are viewing is Augur's reference guide, which means it is information-oriented and targeted at users who just need info about how Augur works. - + * If you have a question about how to achieve a specific goal with Augur, check out our :doc:`Augur-focused How-to Guides section ` in the main Nextstrain documentation. - * If you want to learn the basics of how to use Augur from scratch, check out our :doc:`Zika tutorial ` in the main Nextstrain documentation. + * If you want to learn the basics of how to use Augur from scratch, check out our :doc:`Zika tutorial ` in the main Nextstrain documentation. * If you want to understand how Augur fits together with Auspice to visualize results, check out our :doc:`Data Formats section ` in the main Nextstrain documentation. diff --git a/docs/tutorials/tb_tutorial.md b/docs/tutorials/tb_tutorial.md deleted file mode 100644 index c40b2965e..000000000 --- a/docs/tutorials/tb_tutorial.md +++ /dev/null @@ -1,547 +0,0 @@ -# MTb tutorials using vcf data - -This tutorial explains how to create a Nextstrain build for Tuberculosis sequences. -However, much of it will be applicable to any run where you are starting with [VCF](https://en.wikipedia.org/wiki/Variant_Call_Format) files rather than [FASTA](https://en.wikipedia.org/wiki/FASTA_format) files. - -As in the Zika fasta-input [tutorial](zika_tutorial), we'll build up a Snakefile step-by-step for each step of the analysis. - -## Setup - -To run this tutorial you'll need to [install augur](../guides/install/augur_install.md) and [install Snakemake](https://snakemake.readthedocs.io/en/stable/getting_started/installation.html). - -## Build steps -Nextstrain builds typically require the following steps: - -* Preparing data -* Constructing a phylogeny -* Annotating the phylogeny -* Exporting the results - -However, the commands that make up each step can vary depending on your pathogen and your analysis. -Here, we'll follow these steps: - - 1. Prepare pathogen sequences and metadata - 1. Filter the sequences (remove unwanted sequences and/or sample sequences) - 2. Mask the sequences (exclude regions of the sequence that are unreliable) - - 2. Construct a phylogeny - 1. Construct an initial tree (to get topology) - 2. Convert this into a time-resolved tree (to get accurate branch lengths) - - 3. Annotate the phylogeny - 1. Infer ancestral sequences - 2. Translate genes and identify amino-acid changes - 3. Reconstruct ancestral states (like location or host) - 4. Identify clades on the tree - 5. Identify drug resistance mutations - - 4. Export the final results into auspice-readable format - -## Download Data - -The data in this tutorial is public and is a subset of the data from Lee et al.'s 2015 paper [Population genomics of *Mycobacterium tuberculosis* in the Inuit](http://www.pnas.org/content/112/44/13609). -As location was anonymized in the paper, location data provided here was randomly chosen from the region for illustrative purposes. - -First, download the Tuberculosis (TB) build which includes example data and a pathogen build script. -Then enter the directory you just cloned. - -```bash -git clone https://github.com/nextstrain/tb.git -cd tb -conda activate nextstrain -``` - -## Prepare the Sequences - -A Nextstrain build with VCF file input starts with: - -* A VCF file containing all the sequences you want to include (variable sites only) -* A FASTA file of the reference sequence to which your VCF was mapped -* A tab-delimited metadata file - -There are other files you will need if you want to perform certain steps, like masking. -If you are also working with TB sequences, you may be able to use the files provided here on your own data. -Otherwise, you'll need to provide files specific to your pathogen. - -Here, our input file is compressed with gzip - you can see it ends with `.vcf.gz`. -However, `augur` can take gzipped or un-gzipped VCF files. -It can also produce either gzipped or un-gzipped VCF files as output (detected from the file ending you provide). -Here, we'll usually keep our output VCF files gzipped, by giving our output files endings like `.vcf.gz`, but you can specify `.vcf` instead. - -All the data you need to make the TB build are in the `data` and `config` folders. - -### Filter the Sequences - -Sometimes you may want to exclude certain sequences from analysis. -You may also wish to downsample your data based on certain criteria. `filter` lets you do this. - -For this example, we'll just exclude sequences in the file `dropped_strains.txt`. - -We'll need to specify these starting files at the top of our Snakefile: -```bash -seq_file = "data/lee_2015.vcf.gz", -meta_file = "data/meta.tsv", -exclude_file = "config/dropped_strains.txt" -``` - -And we'll add this as our first rule: -```bash -rule filter: - input: - seq = seq_file, - meta = meta_file, - exclude = exclude_file - output: - "results/filtered.vcf.gz" - shell: - """ - augur filter --sequences {input.seq} \ - --metadata {input.meta} \ - --exclude {input.exclude} \ - --output {output} - """ -``` - -Now run filter. If you are using the Snakefile included with the TB tutorial, you can run: -```bash -snakemake --cores 1 filter -``` - -If you have created your own Snakefile, you'll need to specify its name. For example, if it is called `TB_snakefile`, you would run: -```bash -snakemake --cores 1 -s TB_snakefile filter -``` - -### Mask the Sequences - -There may be regions in your pathogen sequences that are unreliable. -For example, areas that are hard to map because of repeat regions. -Often, these are excluded from analysis so that incorrect calls in these areas don't influence the results. -The areas to be masked are specified in a BED-format file. -This is a standard, tab-delimited format with five columns: Chrom, ChomStart, ChromEnd, locus tag, and Comment. -You can open up `config/Locus_to_exclude_Mtb.bed` in the TB tutorial to see the file format. - -The first, fourth, and fifth columns (Chrom, locus tag, and Comment) can be blank or contain anything - they will be ignored. -All sites between each ChromStart and ChromEnd will be removed from the analysis. - -We'll need to add this BED-format file to the top of the Snakefile (below the files already there): -```bash -mask_file = "config/Locus_to_exclude_Mtb.bed" -``` - -Now we can add the `mask` rule: -```bash -rule mask: - input: - seq = rules.filter.output, - mask = mask_file - output: - "results/masked.vcf.gz" - shell: - """ - augur mask --sequences {input.seq} \ - --mask {input.mask} \ - --output {output} - """ -``` - -## Construct the Phylogeny - -Now our sequences are ready to start analysis. - -With VCF files, we'll do this in two steps that are slightly different from FASTA-input. -1. First, we'll use only the variable sites to construct a tree quickly. This will give us the topology, but the branch lengths will be incorrect. -2. Next, we'll consider the entire sequence to correct our branch lengths. -At the same time, the sample date information will be used to create a time-resolved tree. - -### Get the Topology - -You can use different tree-building programs to build your initial tree, and specify some parameters. -Here, we'll use IQTree. -We specify it here with the argument `--method`, but it's also the default. - -In `tree`, we pass in the VCF file and the reference it was mapped to. -We also pass in a list of sites that we'd like to exclude from building the topology (optional). -These are sites associated with drug-resistance mutations that can influence the topology. -We exclude them here, but they'll be allowed to influence branch length and be included in ancestral sequence reconstruction later. - -We must add the reference sequence our VCF file was mapped to, and our list of sites to exclude from tree-building to the top of the Snakefile: -```bash -ref_file = "data/ref.fasta" -sites_file = "config/drm_sites.txt" -``` - -And add the `tree` rule to the Snakefile: -```bash -rule tree: - input: - aln = rules.mask.output, - ref = ref_file, - sites = sites_file - output: - "results/tree_raw.nwk" - params: - method = 'iqtree' - shell: - """ - augur tree --alignment {input.aln} \ - --vcf-reference {input.ref} \ - --method {params.method} \ - --exclude-sites {input.sites} \ - --output {output} - """ -``` - -### Fix Branch Lengths & Get a Time-Resolved Tree - -Now we'll use the topology from `tree`, but get more accurate branch lengths and a time-resolved tree. -This adjusts branch lengths in the tree to position tips by their sample date and infer the most likely time of their ancestors, using [TreeTime](https://github.com/neherlab/treetime). -There are many options that can be specified here in `refine` to help you get a good tree. - -`refine` will produce as output: -* another tree (newick format) -* a JSON format file with the inferred dates and mutations on each node/branch - -```bash -rule refine: - input: - tree = rules.tree.output, - aln = rules.mask.output, - metadata = meta_file, - ref = ref_file - output: - tree = "results/tree.nwk", - node_data = "results/branch_lengths.json", - params: - root = 'min_dev', - coal = 'opt' - shell: - """ - augur refine --tree {input.tree} \ - --alignment {input.aln} \ - --vcf-reference {input.ref} \ - --metadata {input.metadata} \ - --timetree \ - --root {params.root} \ - --coalescent {params.coal} \ - --output-tree {output.tree} \ - --output-node-data {output.node_data} - """ -``` - -In addition to assigning times to internal nodes, the `refine` command filters tips that are likely outliers. -Branch lengths in the resulting Newick tree measure adjusted nucleotide divergence. -All other data inferred by TreeTime is stored by strain or internal node name in the JSON file. - -## Annotate the Phylogeny - -Now that we have an accurate tree and some information about the ancestral sequences, we can annotate some interesting data onto our phylogeny. -TreeTime can infer ancestral sequences and ancestral traits from an existing phylogenetic tree and metadata to annotate each tip of the tree. - -### Infer Ancestral Sequences - -We can reconstruct the ancestral sequences for the internal nodes on our phylogeny and identify any nucleotide mutations on the branches leading to any node in the tree. - -For VCF runs, `ancestral` will produce another VCF that contains the reconstructed sequence of all the internal nodes and the sequences from the tip nodes, as well as a JSON-format file that contains nucleotide mutation information for each node. - -```bash -rule ancestral: - input: - tree = rules.refine.output.tree, - alignment = rules.mask.output, - ref = ref_file - output: - nt_data = "results/nt_muts.json", - vcf_out = "results/nt_muts.vcf" - params: - inference = "joint" - shell: - """ - augur ancestral --tree {input.tree} \ - --alignment {input.alignment} \ - --vcf-reference {input.ref} \ - --inference {params.inference} \ - --output-node-data {output.nt_data} \ - --output-vcf {output.vcf_out} - """ -``` - -### Identify Amino-Acid Mutations - -With `translate` we can identify amino acid mutations from the nucleotide mutations and a GFF file with gene coordinate annotations. -The resulting JSON file contains amino acid mutations indexed by strain or internal node name and by gene name. -`translate` will also produce a VCF-style file with the amino acid changes for each gene and each sequence, and FASTA file with the translated 'reference' genes which the VCF-style file 'maps' to. - -Because of the number of genes in TB, we will only translate genes associated with drug resistance to save time. -We can pass in a list of genes to translate using `--genes`. -Note that the `--reference-sequence` option is how you pass in the GFF file with the gene coordinates. - -We'll need to add the GFF file with the gene annotations and the file with a list of genes to translate to the list of files at the top of the Snakefile: -```bash -generef_file = "config/Mtb_H37Rv_NCBI_Annot.gff", -genes_file = "config/genes.txt" -``` - -```bash -rule translate: - input: - tree = rules.refine.output.tree, - ref = ref_file, - gene_ref = generef_file, - vcf = rules.ancestral.output.vcf_out, - genes = genes_file - output: - aa_data = "results/aa_muts.json", - vcf_out = "results/translations.vcf", - vcf_ref = "results/translations_reference.fasta" - shell: - """ - augur translate --tree {input.tree} \ - --vcf-reference {input.ref} \ - --ancestral-sequences {input.vcf} \ - --genes {input.genes} \ - --reference-sequence {input.gene_ref} \ - --output-node-data {output.aa_data} \ - --alignment-output {output.vcf_out} \ - --vcf-reference-output {output.vcf_ref} - """ -``` - -### Reconstruct Ancestral States - -`traits` can reconstruct the probable ancestral state of traits like location and host (or others). -This is done by specifying a column or columns in the meta-data file. - -`--confidence` will give confidence estimates for the reconstructed states. -The output will be a JSON file with the state (and confidence, if specified) information for each node. - -```bash -rule traits: - input: - tree = rules.refine.output.tree, - meta = meta_file - output: - "results/traits.json" - params: - traits = 'location' - shell: - """ - augur traits --tree {input.tree} \ - --metadata {input.meta} \ - --columns {params.traits} \ - --output-node-data {output} - """ -``` - -### Identify Specified Clades - -In the [original paper](http://www.pnas.org/content/112/44/13609), the authors identified 'sublineages' within the dataset. -We can add these to our dataset as 'clades' by defining the sublineages with amino-acid or nucleotide mutations specific to that sublineage, given here in the file `config/clades.tsv`. -Open it up in a text editor to have a look at the format. - -The `clades.tsv` file must be tab-delimited with four columns: clade, gene, site, and alt. -The 'clade' column gives the name of the clade being defined - you can have more than one row per clade - it will only be defined from the branch where all criteria are met. -The 'gene' and 'site' columns specify the gene (or `nuc` for nucleotide) and location (by AA position in the gene, or nucleotide position in the genome) where the branch must have the 'alt' (4th column) value to be considered this clade. - -As clades, these sublineages will be labelled and we'll be able to color the tree by them. - -You can specify clades for your own data by first doing a run without clades, then mousing over branches where you'd like to start defining a clade to see what mutations are present. - -We'll need to add the file that defines the clades to the top of our Snakefile: -```bash -clades_file = "config/clades.tsv" -``` - -```bash -rule clades: - input: - tree = rules.refine.output.tree, - aa_muts = rules.translate.output.aa_data, - nuc_muts = rules.ancestral.output.nt_data, - clades = clades_file - output: - clade_data = "results/clades.json" - shell: - """ - augur clades --tree {input.tree} \ - --mutations {input.nuc_muts} {input.aa_muts} \ - --clades {input.clades} \ - --output-node-data {output.clade_data} - """ -``` - -### Identify Drug Resistance Mutations - -`sequence-traits` can identify any trait associated with particular nucleotide or amino-acid mutations, not just drug resistance mutations (DRMs). - -This dataset doesn't actually contain any drug resistance mutations, but identifying such mutations is often of interest to those working on tuberculosis. -Here, we'll run this step as an example, even though it won't add anything to the tree for this dataset. - -Open up the `config/DRMs-AAnuc.tsv` file to see the format of a file that specifies sequence traits. -It contains five columns: GENE, SITE, ALT, DISPLAY_NAME, and FEATURE. DISPLAY_NAME can be blank. - -For drug resistance, we list the gene, the AA position in the gene, the AA mutation that confers resistance (you can list a site multiple times if multiple bases give resistance), and the name of the drug this mutation gives resistance to: -```bash -GENE SITE ALT DISPLAY_NAME FEATURE -gyrB 461 N Fluoroquinolones -gyrB 499 D Fluoroquinolones -rpoB 432 E Rifampicin -rpoB 432 K Rifampicin -``` -We can leave DISPLAY_NAME blank, as auspice will by default display the gene, site, and original and alternative base. - -For mutations outside of protein-coding genes, we can specify their position using nucleotides: -```bash -GENE SITE ALT DISPLAY_NAME FEATURE -nuc 1472749 A rrs: C904A Streptomycin -nuc 1473246 G rrs: A1401G Amikacin Capreomycin Kanamycin -nuc 1673423 T fabG1: G-17T Isoniazid Ethionamide -nuc 1673425 T fabG1: C-15T Isoniazid Ethionamide -``` -In the literature, these mutations are still referred to by their position within non-protein-coding genes (`rrs`) or location near genes (`-17 fabG1`), not their nucleotide location. -We can ensure auspice displays the more useful common nomenclature by giving entries for the DISPLAY_NAME column. - -`sequence-traits` will return a value for each "feature" - for example, all the mutations on the tree that lead to resistance to Streptomycin. -It will also generate a count either of the total number of "features" each node has (ex: the total number of drugs a sequence is resistant to), or the total number or mutations specified in the file each node has (ex: the total number of DRMs a sequence has, even if some are for the same drug). -You can specify a name for this count using the `--label` argument (here: "Drug_Resistance"). -The `--count` argument value specifies whether to count the number of traits (ex: drugs resistant to) (use `traits`) or number of overall mutations (use `mutations`). - -We'll need to add the file that defines the sequence traits (DRMs) to the top of our Snakefile: -```bash -drms_file = "config/DRMs-AAnuc.tsv" -``` - -```bash -rule seqtraits: - input: - align = rules.ancestral.output.vcf_out, - ref = ref_file, - trans_align = rules.translate.output.vcf_out, - trans_ref = rules.translate.output.vcf_ref, - drms = drms_file - output: - drm_data = "results/drms.json" - params: - count = "traits", - label = "Drug_Resistance" - shell: - """ - augur sequence-traits \ - --ancestral-sequences {input.align} \ - --vcf-reference {input.ref} \ - --translations {input.trans_align} \ - --vcf-translate-reference {input.trans_ref} \ - --features {input.drms} \ - --count {params.count} \ - --label {params.label} \ - --output-node-data {output.drm_data} - """ -``` - -## Export the Results - -Finally, collect all node annotations and metadata and export it all in auspice’s JSON format. -The resulting tree and metadata JSON files are the inputs to the auspice visualization tool. - -The names of the output tree and meta data files are here specified by a rule called `all` at the beginning of our Snakefile. -It should be even before the list of files, and looks like this: -```bash -rule all: - input: - auspice_tree = "auspice/tb_tree.json", - auspice_meta = "auspice/tb_meta.json" -``` - -This rule tells Snakemake what the final output of our entire run should look like. -It will run all rules necessary to produce these files, so they should be the names of your final step. -If you have an "all" rule, you can run your entire analysis just by running `snakemake --cores 1` or `snakemake --cores 1 --snakefile Snakefile2` (if the name of your Snakefile is not 'Snakefile'). - -We'll need to add a few remaining files to our list of files at the start of our Snakefile: -```bash -colors_file = "config/color.tsv", -config_file = "config/config.json", -geo_info_file = "config/lat_longs.tsv" -``` - -The `color.tsv` file is optional, but allows us to specify our own colors for particular traits. -If you open it up, you can see that we choose our own colors for values in 'region', 'country', 'location' and 'clade_membership'. -If you don't supply a `color.tsv` file, auspice will choose colors for you. -This can be the simplest way to start - then you can add colors for any traits where you don't like what auspice has chosen. - -The `lat_longs.tsv` file contains the latitudes and longitudes for the geographic locations of your data, and may or may not be needed for your data. -Augur contains many latitudes and longitudes for countries and regions, but if you want to specify data at a different level (state, province, county, city), you can include your own file as well (it will be used in addition to the defaults, so country location can still be retrieved from the augur file, for example). -At the bottom of the `config/lat_longs.tsv` file in the TB tutorial, notice there are entries for 'location', listing each village. - -The `config.json` file should be familiar from the Zika Fasta-input [tutorial](zika_tutorial). -There are a couple of changes worth pointing out, though. - -Since all the samples come from the region of North America and the country of Canada, we don't include these anywhere in our data - all the samples would be the same. -Instead, we have 'location' as a `color_options` entry, and also as our `geo` (where the samples will be drawn on the map), and as a `filters` option. - -We also have a `color_options` entry for 'clade_membership', since we designated clades with the `clades` rule. -The trait is added to our tree as `clade_membership` which is why this is the name of the option and the `key` value, but we could set the `legendTitle` and `menuItem` to be anything we wish, if we wanted. - - -```bash -rule export: - input: - tree = rules.refine.output.tree, - metadata = meta_file, - branch_lengths = rules.refine.output.node_data, - traits = rules.traits.output, - nt_muts = rules.ancestral.output.nt_data, - aa_muts = rules.translate.output.aa_data, - drms = rules.seqtraits.output.drm_data, - color_defs = "config/colors.tsv", - config = "config/config.json", - geo_info = "config/lat_longs.tsv", - clades = rules.clades.output.clade_data - output: - auspice_json = "auspice/tb.json", - shell: - """ - augur export v2 \ - --tree {input.tree} \ - --metadata {input.metadata} \ - --node-data {input.branch_lengths} {input.traits} {input.drms} {input.aa_muts} {input.nt_muts} {input.clades} \ - --auspice-config {input.config} \ - --colors {input.color_defs} \ - --lat-longs {input.geo_info} \ - --output {output.auspice_json} \ - """ -``` - - -As mentioned previously, this dataset has no drug resistance, so it's not included in the `config.json` file to display, even though we ran the `sequence-traits` rule. -If you did have drug resistance information that you wanted to display, you would need to add it to the `config.json` file as `color_options`. - -First, you would want to add a color-by for the total number of drugs each node is resistant to. -Since we gave the label 'Drug_Resistance' when we ran the rule, this will be the name of the option, and the `key`, but we can make the `menuItem` and `legendTitle` different if we wish: -``` - "Drug_Resistance": { - "menuItem": "Drug_Resistance", - "legendTitle": "Drug Resistance", - "type": "discrete", - "key": "Drug_Resistance" - }, -``` -If you had given a different label when you ran the rule, you would change this entry to match. - -You would then need an option for each drug where you have resistance information (or each FEATURE where you have information). -For example, to show the mutations present that confer resistance to Streptomycin and Rifampicin: -``` - "Streptomycin": { - "menuItem": "Streptomycin", - "legendTitle": "Streptomycin Resistance", - "type": "discrete", - "key": "Streptomycin" - }, - "Rifampicin": { - "menuItem": "Rifampicin", - "legendTitle": "Rifampicin Resistance", - "type": "discrete", - "key": "Rifampicin" - }, -``` -You would need an entry for every FEATURE in your original file (though you could then remove any that had no information on the tree). - diff --git a/docs/tutorials/zika_tutorial.md b/docs/tutorials/zika_tutorial.md deleted file mode 100644 index f7bafb66b..000000000 --- a/docs/tutorials/zika_tutorial.md +++ /dev/null @@ -1,355 +0,0 @@ -# Zika tutorial -- fasta input - -The tool-chain [*augur*](https://github.com/nextstrain/augur) is the bioinformatics engine of nextstrain and produces the files that can be visualized in the webbrowser using [*auspice*](https://github.com/nextstrain/auspice). -Augur consists of a number of tools that allow the user to filter and align sequences, build trees, and integrate the phylogenetic analysis with meta data. -The different tools are meant to be composable and the output of one tool will serve as the input of other tools. -We will work off the tutorial for Zika virus on the [nextstrain web site](https://nextstrain.org/docs/getting-started/zika-tutorial) and the github repository [nextstrain/zika-tutorial](https://github.com/nextstrain/zika-tutorial). - -## Setup - -To run this tutorial you'll need to [install augur](../guides/install/augur_install.md) and [install Snakemake](https://snakemake.readthedocs.io/en/stable/getting_started/installation.html). - -## Augur commands - -As an example, we'll look that the `filter` command in greater detail. -This command allows you to selected various subsets of your input data for different types of analysis. -A simple example use of this command would be -```bash -augur filter --sequences data/sequences.fasta --metadata data/metadata.tsv --min-date 2012 --output filtered.fasta -``` -This command will select all sequences with collection date in 2012 or later. -The filter command has a large number of options that allow flexible filtering for many common situations. -One such use-case is the exclusion of sequences that are known to be outliers (e.g.~because of sequencing errors, cell-culture adaptation, ...). -These can be specified in a separate file: -``` -BRA/2016/FC_DQ75D1 -COL/FLR_00034/2015 -... -``` -To drop such strains, you can pass the name of this file to the augur filter command: -```bash -augur filter --sequences data/sequences.fasta \ - --metadata data/metadata.tsv \ - --min-date 2012 \ - --exclude config/dropped_strains.txt \ - --output filtered.fasta -``` -(To improve legibility, we have wrapped the command across multiple lines.) -If you run this command (you should be able to copy-paste this into your terminal), you should see that one of the sequences in the data set was dropped since its name was in the `drooped_strain.txt` file. - -Another common filtering operation is subsetting of data to a achieve a more even spatio-temporal distribution or cut-down data set size to more manageable numbers. -The filter command allows you to select a specific number of sequences from specific groups, for example one sequence per month from each country: -```bash -augur filter \ - --sequences data/sequences.fasta \ - --metadata data/metadata.tsv \ - --min-date 2012 \ - --exclude config/dropped_strains.txt \ - --group-by country year month \ - --sequences-per-group 1 \ - --output filtered.fasta -``` -This subsampling and filtering will reduce the number of sequences in this tutorial data set from 34 to 24. - -## Chaining of augur commands with snakemake. - -The output of one augur command serves as input to the next and hence these commands need to be executed in order. -This is a common pattern in bioinformatics and multiple tools -- so called workflow managers -- exist to facilitate this. -Within nextstrain, we have relied on [Snakemake](https://snakemake.readthedocs.io/en/stable/). - -Snakemake breaks a workflow into a set of rules that are specified in a file called `Snakefile`. -Each rule takes a number of input files, specifies a few parameters, and produces output files. -A simple rule would look like this: -```python -rule filter: - input: - sequences = "data/sequences.fasta", - metadata = "data/metadata.tsv" - output: - sequences = "results/filtered.fasta" - params: - min_date = 2012 - shell: - """ - augur filter \ - --sequences {input.sequences} \ - --metadata {input.metadata} \ - --output {output.sequences} \ - --min-date {params.min_date} - """ -``` -This rule would produce `results/filtered.fasta` from the input files `data/sequences.fasta` and `data/metadata.tsv` using the `augur filter` command. -Note that we explicitly specify what is an input and what is an output file. -To filter our data, we would now call snakemake as -```bash -snakemake --cores 1 results/filtered.fasta -``` -and snakemake will run the same command as specified above. - -So far, this is just a complicated reformulation of what we did above, but the benefit of workflow management becomes clear once we add more steps. -The next natural step in our phylogenetic pipeline is aligning the filtered sequences and we define a rule `align`. -```bash -rule align: - input: - sequences = rules.filter.output.sequences, - reference = "config/zika_outgroup.gb" - output: - alignment = "results/aligned.fasta" - shell: - """ - augur align \ - --sequences {input.sequences} \ - --reference-sequence {input.reference} \ - --output {output.alignment} - """ -``` -If you now want to generate the alignment, you can type -```bash -snakemake --cores 1 results/aligned.fasta -``` -and snakemake will - - * determine that `results/aligned.fasta` is an output of rule `align` - * check whether all required input files are in place and run the necessary rules if not - * run rule `align` and check whether the file `results/aligned.fasta` appeared. - -If you supply a reference sequence to `augur align`, augur will include that reference sequence in the alignment and strip all insertions relative to that reference. -This will guarantee a consistent reference coordinate system and genome annotation later on. - - -### Construct the Phylogeny - -Infer a phylogenetic tree from the multiple sequence alignment. -```bash -rule tree: - input: - alignment = rules.align.output.alignment - output: - tree = "results/tree_raw.nwk" - shell: - """ - augur tree \ - --alignment {input.alignment} \ - --output {output.tree} - """ -``` - -The resulting tree is stored in [Newick format](http://evolution.genetics.washington.edu/phylip/newicktree.html). -Branch lengths in this tree measure nucleotide divergence. - -### Get a Time-Resolved Tree -Most phylogenies you see on nextstrain are time-resolved, that is the branch lengths of the tree correspond to calendar time rather than evolutionary distance. -The default method to infer time-resolved phylogenies in nextstrain is [TreeTime](https://github.com/neherlab/treetime) and this analysis is an done using the augur command `refine` (this steps "refines" the existing tree...). -The corresponding rule for in the snakefile would look as follows: -```bash -rule refine: - input: - tree = rules.tree.output.tree, - alignment = rules.align.output, - metadata = "data/metadata.tsv" - output: - tree = "results/tree.nwk", - node_data = "results/branch_lengths.json" - shell: - """ - augur refine \ - --tree {input.tree} \ - --alignment {input.alignment} \ - --metadata {input.metadata} \ - --timetree \ - --output-tree {output.tree} \ - --output-node-data {output.node_data} - """ -``` -This command requires the output of the tree and align rules as input and will produce the file `results/tree.nwk` and `results/branch_lengths.json`. -The latter contains the inferred clock model and the inferred dates of all nodes, including internal nodes. -Each internal node has been given a unique name (e.g.~`NODE_0000001211`) which we will use in the following steps to attach more information to the tree. -The `refine` command has many different options that allow to specify clock rates, filter sequences that don't follow the molecular clock, calculate confidence intervals, etc. -We'll get into these details later. - -## Annotate the Phylogeny - -### Reconstruct Ancestral Traits - -TreeTime can also infer ancestral traits from an existing phylogenetic tree and metadata annotating each tip of the tree. -The following command infers the region and country of all internal nodes from the time tree and original strain metadata. -As with the `refine` command, the resulting JSON output is indexed by strain or internal node name. - -Specifying `--confidence` means that the confidence intervals for the reconstructed trait values will be estimated. - -```bash -rule traits: - input: - tree = rules.refine.output.tree, - metadata = "data/metadata.tsv" - output: - node_data = "results/traits.json", - params: - columns = "region country" - shell: - """ - augur traits \ - --tree {input.tree} \ - --metadata {input.metadata} \ - --output-node-data {output.node_data} \ - --columns {params.columns} \ - --confidence - """ -``` - -### Infer Ancestral Sequences - -Using `ancestral`, we can reconstruct what the ancesters of our samples' sequences must have looked like, and record the mutations that occurred on each branch. - -Next, infer the ancestral sequence of each internal node and identify any nucleotide mutations on the branches leading to any node in the tree. - -```bash -rule ancestral: - input: - tree = rules.refine.output.tree, - alignment = rules.align.output - output: - node_data = "results/nt_muts.json" - shell: - """ - augur ancestral \ - --tree {input.tree} \ - --alignment {input.alignment} \ - --output-node-data {output.node_data} - """ -``` - -### Identify Amino-Acid Mutations - -Identify amino acid mutations from the nucleotide mutations and a reference sequence with gene coordinate annotations. -The resulting JSON file contains amino acid mutations indexed by strain or internal node name and by gene name. -To export a FASTA file with the complete amino acid translations for each gene from each node’s sequence, specify the `--alignment-output` parameter in the form of `results/aligned_aa_%GENE.fasta`. - -The reference sequence needs to be in Genbank format (if working with Fasta input) or GFF format (if working with VCF input) - see the Prerequisites page for more information on how to find and download a reference sequence from Genbank. - -Open the Zika tutorial reference sequence `config/zika_outgroup.gb` in a text editor. -For translation, the important parts are the `source` and `CDS` sections that are part of the `Features` part of the file: - -```bash -[...] -FEATURES Location/Qualifiers - source 1..10769 - /collection_date="25-Oct-2013" - /country="French Polynesia" - /db_xref="taxon:64320" - /host="Homo sapiens" - /isolation_source="serum" - /mol_type="genomic RNA" - /organism="Zika virus" - /strain="PF13/251013-18" -[...] - CDS 91..456 - /product="capsid protein" - /gene="CA" - CDS 457..735 - /product="propeptide" - /gene="PRO" - CDS 736..960 - /product="membrane protein" - /gene="MP" - CDS 961..2472 - /product="envelope protein" - /gene="ENV" -[...] -``` - -The `source` section tells `augur` how long the whole genome is. -The `CDS` sections describe each gene (using `/gene=`) that should be translated by giving its name and its start and end location in the genome. - -In your downloaded file from Genbank, you might notice that it has only one or no `CDS` sections, and that instead the sections you would like to translate as genes are described as `mat_peptide`, like the below: - -``` - mat_peptide 91..465 - /product="capsid protein" - mat_peptide 466..735 - /product="propeptide" - mat_peptide 736..960 - /product="membrane protein" - mat_peptide 961..2475 - /product="envelope protein" -``` - -In this case, you will have to manually modify this file to work in `augur`. -Anything you wish to translate as a gene should be changed from `mat_peptide` to `CDS` as its section name. -You will also need to add `/gene=` and a gene name to the sections you'd like to translate. -If you have one really long `CDS` section that covers most of the genome, we'd recommend not including this to translate (don't provide a `/gene=`) if you can include smaller genes instead, as such a long polyprotein is not very informative. - -If you are working with VCF-input, the GFF file contains the same information, but `augur` will only translate those items coded as `gene`. - -```bash -rule translate: - input: - tree = rules.refine.output.tree, - node_data = rules.ancestral.output.node_data, - reference = "config/zika_outgroup.gb" - output: - node_data = "results/aa_muts.json" - shell: - """ - augur translate \ - --tree {input.tree} \ - --ancestral-sequences {input.node_data} \ - --reference-sequence {input.reference} \ - --output-node-data {output.node_data} \ - """ -``` - -## Export the Results - -Finally, collect all node annotations and metadata and export it all in auspice’s JSON format. -This command pulls together the metadata file and all of the output files generated by the previous steps, and combines it into two output files that contain the annotated tree (`"auspice/zika_tree.json"`) and associated metadata (`"auspice/zika_meta.json"`). - -There is one input file below that is not metadata or output from a previous step: the config file (`"config/auspice_config.json"`). -This file is very important - it species the page title, maintainer(s), what filters are present, what places should be mapped by, and what the user is able to color the tree by ("color-bys"). - -Open up `"config/auspice_config.json"` in a text editor. -The `title` and `maintiner` fields should be easy to spot. -You can modify them if you wish (and should, if you are using your own data!). - -The `filters` section specifies what traits will be available for users to filter the data by. -For example, if "country" is a filter, users will be able to only show data from the countries they select. - -The `color_options` section specifies what users will be allowed to color the data by. You will almost always have the first two sections: `"gt"` and `"num_date"`, unless you are using trees without sequences or trees that aren't time-resolved. -After that, you can include any trait in your data. -If you have additional traits in your metadata, you will need to add them here to have them show up in auspice! You can copy the entry for `"country"` and re-name it to match the appropriate column from your metadata file. - -The `geo` section informs auspice what traits should be used to draw the samples onto the map. Obviously, they must be location-related, like "country" and "region". -You can use other locations if that's more relevant to your data. We won't cover this right now, but can explain in detail as you come to this with your own data. - -Don't worry about the `defaults` section for the moment. - - -The resulting tree and metadata JSON files are the inputs to the auspice visualization tool. - -```bash -rule export: - input: - tree = rules.refine.output.tree, - metadata = "data/metadata.tsv", - branch_lengths = rules.refine.output.node_data, - traits = rules.traits.output.node_data, - nt_muts = rules.ancestral.output.node_data, - aa_muts = rules.translate.output.node_data, - auspice_config = "config/auspice_config.json" - output: - auspice_json = "auspice/zika.json", - shell: - """ - augur export v2 \ - --tree {input.tree} \ - --metadata {input.metadata} \ - --node-data {input.branch_lengths} {input.traits} {input.nt_muts} {input.aa_muts} \ - --colors {input.colors} \ - --lat-longs {input.lat_longs} \ - --auspice-config {input.auspice_config} \ - --output {output.auspice_json} - """ -``` - - - diff --git a/docs/usage/cli/filter.rst b/docs/usage/cli/filter.rst index 3ebf97072..36250aee5 100644 --- a/docs/usage/cli/filter.rst +++ b/docs/usage/cli/filter.rst @@ -11,11 +11,11 @@ augur filter :func: make_parser :prog: augur :path: filter - + How we subsample sequences in the zika-tutoral ============================================== -As an example, we'll look that the ``filter`` command in greater detail using material form the :doc:`zika tutorial `. +As an example, we'll look that the ``filter`` command in greater detail using material from the :doc:`Zika tutorial `. The filter command allows you to selected various subsets of your input data for different types of analysis. A simple example use of this command would be @@ -45,7 +45,7 @@ To drop such strains, you can pass the name of this file to the augur filter com --output filtered.fasta (To improve legibility, we have wrapped the command across multiple lines.) -If you run this command (you should be able to copy-paste this into your terminal) on the data provided in the :doc:`zika tutorial `, you should see that one of the sequences in the data set was dropped since its name was in the ``dropped_strains.txt`` file. +If you run this command (you should be able to copy-paste this into your terminal) on the data provided in the :doc:`Zika tutorial `, you should see that one of the sequences in the data set was dropped since its name was in the ``dropped_strains.txt`` file. Another common filtering operation is subsetting of data to a achieve a more even spatio-temporal distribution or to cut-down data set size to more manageable numbers. The filter command allows you to select a specific number of sequences from specific groups, for example one sequence per month from each country: