Tutorial
This guide is a brief workshop for gaining hands-on experience with Mycotools.
The following commands were executed against the reference MycotoolsDB in the
test/ directory of the Mycotools
repo. Full descriptions of each command are
available from its help menu (-h/--help) and in the usage
guide.
Note: Mycotools v2 consolidated the former standalone scripts under two commands —
mtdbfor the database andmycotoolsfor analysis. The old names (fa2tree,db2search, …) still forward for one transition release, but print a deprecation notice. If a command does not behave as expected, confirm the current invocation with-hor in the usage guide.
- A UNIX/Linux command-line environment. Mycotools is developed for Linux, though
it should be compatible with UNIX-based BASH environments, Linux virtual
machines, and potentially Windows Subsystem for Linux.
⚠️ Macs with the M* series CPUs do not have native package support for many of the dependencies. On an Intel Mac, switch to BASH as your default shell:chsh -s /bin/bash - Install Mycotools using miniconda3 as an environment manager
- Install FigTree for viewing phylogenetic trees
- Basic bash shell navigation knowledge
- Setup an NCBI account
- Setup a JGI account and review the use-restricted data terms and conditions
Install the following in your activated Mycotools environment (conda activate mycotools):
conda install clinker-py -c conda-forge -c bioconda
conda install orthofinder -c biocondaOne of the best habits while programming is to make clean, organized directory
hierarchies for projects. I recommend an overarching projects directory with
subdirectories titled <PROJECT_NAME>_YYYYmm, where YYYYmm is the year and
month. Navigate somewhere you want to start your directories (Desktop or Home)
and initialize the hierarchy:
mkdir <PROJECT_DIR>/mycotools_<YYYYmm>
cd <PROJECT_DIR>/mycotools_<YYYYmm>curl -o reference.mtdb https://raw.githubusercontent.com/xonq/mycotools/master/test/reference.mtdbThere are multiple ways to initialize a MycotoolsDB (MTDB): from a de novo
assimilation of publicly available data (mtdb update -i <INIT_DIR>); from a
spreadsheet of metadata referencing local genomes (mtdb update -i <INIT_DIR> --predb <PREDB.tsv>); or from a reference .mtdb file that contains publicly
available assembly accessions. We will implement a reference.
# access mycotools by making sure your environment is activated
conda activate mycotools
# initialize the mycotoolsdb according to a reference file
mtdb update -i <PROJECT_DIR>/mycotools_<YYYYmm>/ -r reference.mtdbYour own genomes are added by filling out a spreadsheet of metadata that
references them. Both an assembly and a gff3 annotation are required. In this
example, we will add a single genome acquired from JGI — Mycotools can
automatically assimilate all JGI genomes via mtdb update -i <INIT_DIR>, but
that takes quite some time.
# make a directory to run these commands
mkdir jgi_dwnld_<YYYYmm>
cd jgi_dwnld_<YYYYmm>
# download the assembly and gff of a MycoCosm accession
mycotools download jgi -i Ustbr1 -a -gThis outputs a predb file that is ready for assimilation into the database. To
add your own genomes, you would fill out one of these files manually — generate a
blank copy via mtdb predb > predb.tsv — then run the same commands we run here:
# curate the data via mtdb predb
mtdb p Ustbr1.predb.tsv
# add the curated data to the primary MTDB
mtdb u -a predb_<YYYYmmdd>/predb.mtdbNow we can check whether the file was added by querying its genome code:
mtdb ustbro1mtdb controls your interface to generated MTDBs. You can maintain different
primary MTDBs for different projects by initializing separate databases, and/or
keep one large primary database and interact with subsets of it.
To print the path of the primary database you are linked to:
mtdbThat's it! We can now interface with that primary data using some basic shell wizardry, e.g. to preview the contents of the primary MTDB:
cat $(mtdb)To interface with a subset of the primary MTDB, extract the genomes of interest:
# extract a database of Ustilago genomes
mtdb extract -l Ustilago > ust.mtdbLet's get into some computational biology! We will start by obtaining basic annotation statistics for what is in the database.
# return to the main project directory if you have not
cd <PROJECT_DIR>/mycotools_<YYYYmm>
# obtain annotation statistics from the primary MTDB
mycotools stats annotation $(mtdb) > annotation_stats.tsvYou can open that .tsv in whatever spreadsheet software or command-line text
editor you prefer. We can run a similar analysis for the assemblies:
mycotools stats assembly $(mtdb) > assembly_stats.tsvBoth commands can also be run against a single genome .gff3 or .fna
respectively, IF it has been curated into an MTDB.
It is often useful to group genes into homology groups for downstream analyses.
Single-copy orthologs (SCOs) are a particular type of homology group that is
useful for phylogenetic reconstruction, because SCO evolution is assumed to be
congruent with the evolution of the species. OrthoFinder is a robust automated
method of SCO determination, but we can more swiftly circumscribe homology groups
and identify putative SCOs using a faster algorithm, MMseqs cluster, implemented
in mycotools cluster db. Let's run a phylogenetic analysis on a subset of our
genomes.
# extract the lineages we want
mtdb e -l Ustilago >> tree.mtdb
mtdb e -l Puccinia >> tree.mtdb
mtdb e -l Trichoderma >> tree.mtdb
mtdb e -l Psilocybe >> tree.mtdb
mycotools cluster db -d tree.mtdbA phylogenomic tree is a phylogenetic analysis of a group of genomes that incorporates an arbitrary minimum number of genes. Mycotools implements the multigene partition method, which determines the best-fit evolutionary model for each SCO, concatenates an alignment of all SCOs in the dataset, and reconstructs a maximum-likelihood phylogenomic tree by applying the best-fit model to each individual gene partition. As a general rule-of-thumb, the more tractable genes we have, the better — we have quite a few, and it would be ideal to use them all!... but we don't have time, so let's work with a subset:
# make a folder for single copy genes
mkdir sco_<YYYYmm>
# copy the top three SCOs
for i in $(ls cluster_db_<YYYYmmdd>/single_copy_genes/ | head -3)
do cp cluster_db_<YYYYmmdd>/single_copy_genes/$i sco_<YYYYmm>/
done
# run the tree building pipeline
mycotools phylo tree -i sco_<YYYYmm>/ --partitionWhen complete, open concatenated.nex.contree — the consensus tree with 1000
ultrafast bootstrap replicates — from the resulting phylo_tree_<YYYYmmdd>
directory in FigTree.
You will note that the tips are labeled with the ome code, but we probably want to see the actual genus, species, and strain names. Let's convert the tips from genome codes to full names:
mycotools rename phylo_tree_<YYYYmmdd>/concatenated.nex.contree o \
> phylo_tree_<YYYYmmdd>/full_name.newickGo ahead and open this one in FigTree, and let's glance at how well supported the phylogeny is.
Horizontal transfer, gene duplication, convergence, and gene loss all lead to discordance between species and gene evolution. To study a gene's evolution, we therefore have to reconstruct a phylogeny of it — let's study a gene associated with nitrate assimilation in our dataset.
First, we need to identify homologs of this gene across our database via a BLAST
search of the protein sequence, which we obtain using mtdb accession fa.
# extract the protein accession of interest
mtdb accession fa -a ustbro1_1795 > ustbro1_1795.faa
# run a blast search on this gene against the primary MTDB
mycotools homology db -a blastp -q ustbro1_1795.faa -e 2With the BLAST results in hand, we can build a phylogeny from the outputted
.fasta of homologs. First we move the phylogenomic analysis to a different
directory, so the tree builder doesn't output to the same place because the date
is the same. This time we will use FastTree to expedite tree building at the cost
of some quality:
# move the phylogenomic directory
mv phylo_tree_<YYYYmmdd> phylogenomic_<YYYYmmdd>/
# run the single gene phylo
mycotools phylo tree -i homology_db_<YYYYmmdd>/fastas/ustbro1_1795.search.fa -fNow, we can view this tree in FigTree.
Oftentimes, phenotypes are derived from multiple genes that are proximal to one another in the linear genome. Studying the evolution of a locus from both a gene-by-gene and a locus-wide (shared synteny) perspective can therefore be important for understanding the evolution of a phenotype.
For this example, we will look at the nitrate assimilation gene cluster, which encodes three genes that import extracellular nitrate and reduce it to intracellular ammonium as a nitrogen source.
The Cluster Reconstruction and phylogenetic Analysis Pipeline (CRAP) facilitates swiftly analyzing the evolution of a gene cluster. We implement it by first extracting a representative nitrate assimilation locus, then inputting it into the pipeline:
# extract a locus of interest, and store in a file
mtdb accession locus -a ustbro1_1795 -p 1 > nitrate_cluster.txt
# run the CRAP analysis
mycotools phylo crap -q nitrate_cluster.txt -s blastpOnce complete, the output directory phylo_crap_<YYYYmmdd> contains .svg files
of the nitrate assimilation gene phylogenies, with synteny diagrams of the
individual genes extracted applied to each tip.
We can also make more visually appealing synteny diagrams that depict percent identity through Clinker. CRAP outputs the most similar loci to the query sequence; we invoke this option by having CRAP rerun off the initial results, referencing the initial output directory. This is the general way of resuming Mycotools commands.
# rerun a Mycotools command referencing a previous output directory
mycotools phylo crap -q nitrate_cluster.txt -s blastp --loci -o phylo_crap_<YYYYmmdd>There should now be a loci/ subdirectory containing the most similar loci, and
we want to generate GenBanks of these loci to run Clinker. Let's make a directory
in our main project folder:
mkdir clinker_<YYYYmmdd>Then generate GenBanks of each locus file using some basic BASH scripting:
for i in phylo_crap_<YYYYmmdd>/loci/*txt; do
o=$(basename ${i} .txt)
mtdb accession gbk -i ${i} > clinker_<YYYYmmdd>/${o}.gbk
doneWe can now run Clinker on this dataset:
cd clinker_<YYYYmmdd>
clinker ./ -p nitrate_assimilation.htmlWe previously circumscribed homologous gene groups using MMseqs2 cluster, which
is a fairly crude method. OrthoFinder is a much more robust method for
circumscribing genes into orthologous gene groups (orthogroups) that are
approximately set with respect to the most recent common ancestral genome of the
dataset. OrthoFinder is not built into Mycotools yet, though we can implement it —
and other external software — by acquiring the necessary files.
# make your OrthoFinder analysis directory
mkdir orthofinder_<YYYYmmdd>/
cd orthofinder_<YYYYmmdd>/
# extract the proteomes you need for an analysis, referencing a MycotoolsDB we
# extracted previously
mtdb files -p -d ../tree.mtdbThis generates a folder of proteome fastas, faa/. Note that it contains
"symlinked" files, which are essentially empty links to the original file;
editing them therefore edits the files in the database. To make hard copies,
append the --hard flag to the mtdb files command.
orthofinder -f faa/If your dataset is small enough, I recommend implementing OrthoFinder over
mycotools cluster db.
BioPython: Chapman, B. & Chang, J. Biopython: Python Tools for Computational Biology. SIGBIO Newsl. 20, 15–19 (2000).
BLAST: Ye, J., McGinnis, S. & Madden, T. L. BLAST: improvements for better sequence analysis. Nucleic Acids Research 34, W6–W9 (2006).
Clinker: clinker & clustermap.js: automatic generation of gene cluster comparison figures | Bioinformatics | Oxford Academic. https://academic.oup.com/bioinformatics/article/37/16/2473/6103786.
ClipKIT: Steenwyk, J. L., Iii, T. J. B., Li, Y., Shen, X.-X. & Rokas, A. ClipKIT: A multiple sequence alignment trimming software for accurate phylogenomic inference. PLOS Biology 18, e3001007 (2020).
CRAP: Slot, J. C. & Rokas, A. Horizontal Transfer of a Large and Highly Toxic Secondary Metabolic Gene Cluster between Fungi. Current Biology 21, 134–139 (2011).
FastTree: Price, M. N., Dehal, P. S. & Arkin, A. P. FastTree 2 – Approximately Maximum-Likelihood Trees for Large Alignments. PLoS ONE 5, e9490 (2010).
IQ-TREE: Minh, B. Q. et al. IQ-TREE 2: New Models and Efficient Methods for Phylogenetic Inference in the Genomic Era. Molecular Biology and Evolution 37, 1530–1534 (2020).
MAFFT: Katoh, K., Misawa, K., Kuma, K. & Miyata, T. MAFFT: a novel method for rapid multiple sequence alignment based on fast Fourier transform. Nucleic Acids Res 30, 3059–3066 (2002).
MMseqs2: Steinegger, M. & Söding, J. MMseqs2 enables sensitive protein sequence searching for the analysis of massive data sets. Nat Biotechnol 35, 1026–1028 (2017).
OrthoFinder: Emms, D. M. & Kelly, S. OrthoFinder: phylogenetic orthology inference for comparative genomics. Genome Biology 20, 238 (2019).