Disclaimer: I banged this out fast from existing scripts to help some folks, but haven’t tested it yet. Will do that shortly, in the meantime, be careful!
I use phylogenetic placement, namely the program pplacer, in a lot of my publications. It is also a core part of of the paprica metabolic inference pipeline. As a result I field a lot questions from people trying to integrate pplacer into their own workflows. Although the Matsen group has done an excellent job with documentation for pplacer, guppy, and taxtastic, the three programs you need to work with to do phylogenetic placement from start to finish (see also EPA), there is still a steep learning curve for new users. In the hope of bringing the angle of that curve down a notch or two, and updating my previous posts on the subject (here and here), here is a complete, start to finish example of phylogenetic placement, using 16S rRNA gene sequences corresponding to the new tree of life recently published by Hug et al. To follow along with the tutorial start by downloading the sequences here.
You can use any number of alignment and tree building programs to create a reference tree for phylogenetic placement. I strongly recommend using RAxML and Infernal. After a lot of experimentation this combination seems to be produce the most correct topologies and best supported trees. You should be aware that no 16S rRNA gene tree (or any other tree) is absolutely “correct” for domain-level let alone life-level analyses, but nothing in life is perfect. While you’re installing software I also recommend the excellent utility Seqmagick. Finally, you will need a covariance model of the 16S rRNA gene to feed into Infernal. You can find that at the Rfam database here.
The workflow will follow these steps:
- Create an alignment of the reference sequences with Infernal
- Create a phylogenetic tree of the alignment
- Create a reference package from the alignment, tree, and stats file
- Proceed with the phylogenetic placement of your query reads
Create an alignment of the reference sequences
The very first thing that you need to do is clean your sequence names of any wonky punctuation. This is something that trips up almost everyone. You should really have no punctuation in the names beyond “_”, and no spaces!tr " -" "_" < hug_tol.fasta | tr -d "%\,;():=.\\*\"\'" > hug_tol.clean.fasta
Next create an alignment from the cleaned file. I always like to degap first, although it shouldn’t matter.## Degap seqmagick mogrify --ungap hug_tol.clean.fasta ## Align cmalign --dna -o hug_tol.clean.align.sto --outformat Pfam 16S_bacteria.cm hug_tol.clean.fasta ## Convert to fasta format seqmagick convert hug_tol.clean.align.sto hug_tol.clean.align.fasta
Build the reference tree
At this point you should have a nice clean DNA alignment in the fasta format. Now feed it to RAxML to build a tree. Depending on the size of the alignment this can take a little bit. I know you’ve read the entire RAxML manual so of course you are already aware that adding additional cpus won’t help…raxmlHPC-PTHREADS-AVX2 -T 8 -m GTRGAMMA -s hug_tol.clean.align.fasta -n ref.tre -f d -p 12345
I like having a sensibly rooted tree; it’s just more pleasing to look at. You can do this manually, or you can have RAxML try to root the tree for you.raxmlHPC-PTHREADS-AVX2 -T 2 -m GTRGAMMA -f I -t RAxML_bestTree.ref.tre -n root.ref.tre
Okay, now comes the tricky bit. Clearly you’d like to have some support values on your reference tree, but the Taxtastic program that we will use to build the reference tree won’t be able to read the RAxML stats file if it includes confidence values. The work around is to build a second tree with confidence values. You will feed this tree to Taxtastic with the stats file from the tree we already generated.## Generate confidence scores for tree raxmlHPC-PTHREADS-AVX2 -T 8 -m GTRGAMMA -f J -p 12345 -t RAxML_rootedTree.root.ref.tre -n conf.root.ref.tre -s hug_tol.clean.align.fasta
Now we can use the alignment, the rooted tree with confidence scores, and the stats file without confidence scores to create our reference package.taxit create -l 16S_rRNA -P hug_tol.refpkg --aln-fasta hug_tol.clean.align.fasta --tree-stats RAxML_info.ref.tre --tree-file RAxML_fastTreeSH_Support.conf.root.ref.tre
Align the query reads
At this point you have the reference package and you can proceed with analyzing some query reads! The first step is to align the query reads in exactly the same fashion as the reference sequences. This is important as the alignments will be merged later.## Clean the names tr " -" "_" < query.fasta | tr -d "%\,;():=.\\*\"\'" > query.clean.fasta ## Remove any gaps seqmagick mogrify --ungap query.clean.fast ## Align cmalign --dna -o query.clean.align.sto --outformat Pfam 16S_bacteria.cm query.clean.fasta
Now we use the esl-alimerge command, included with Infernal, to merge the query and reference alignments.## Merge alignments esl-alimerge --outformat pfam --dna -o query.hug_tol.clean.align.sto query.clean.align.sto hug_tol.refpkg/hug_tol.clean.align.sto ## Convert to fasta seqmagick convert query.hug_tol.clean.align.sto query.hug_tol.clean.align.fast
Now we’re on the home stretch, we can execute the phylogenetic placement itself! The flags are important here, so it’s worth checking the pplacer documentation to insure that your goals are consistent with mine (get a job, publish some papers?). You can probably accept most of the flags for the previous commands as is.pplacer -o query.hug_tol.clean.align.jplace -p --keep-at-most 20 -c hug_tol.refpkg query.hug_tol.clean.align.fasta
At this point you have a file named query.hug_tol.clean.align.jplace. You will need to use guppy to convert this json-format file to information that is readable by human. The two most useful guppy commands (in my experience) for a basic look at your data are:## Generate an easily parsed csv file of placements, with only a single placement reported for each ## query read. guppy to_csv --point-mass --pp -o query.hug_tol.clean.align.csv query.hug_tol.clean.align.jplace ## Generate a phyloxml tree with edges fattened according to the number of placements. guppy fat --node-numbers --point-mass --pp -o query.hug_tol.clean.align.phyloxml query.hug_tol.clean.align.jplace