Recipes
Here we introduce our cookbook which is made of recipes. Each recipe is a customized piece of code to complete different advanced tasks. By using these recipes (or creating new ones) the user can take full profit of functionalities provided by pagoo and easily interact with other R packages to perform a variety of analyses including phylogenetics, pangenome-wide association studies, sequence comparisons, ecological measures, prepare publication-quality figures, among others. If you want to standardize any analysis from your pangenome data that is not covered in this tutorial, we can help you and write a recipe. Contact us!
Start by loading the pagoo object, as in previous tutorials:
library(pagoo) # Load package
toy_rds <- system.file('extdata', 'campylobacter.RDS', package = 'pagoo')
p <- load_pangenomeRDS(toy_rds)
Annotate clusters
A common practice to improve the annotation quality of the pangenome is to retrieve a single representative sequence of each cluster and to ‘blast’ it against a given database. The best hits annotation is then transferred to the clusters. The EggNOG database is a very comprehensive and well curated resource which includes annotation from many other databases (e.g. KEGG, COG, CAZy, etc), and the eggnog-mapper tool allows to easily transfer their annotation to any given query sequence.
This recipe shows how to write a multi-fasta file containing one representative sequence of each cluster to annotate with the EggNOG-mapper tool, and then it shows how to add this new information to the pangenome object.
The first step is to retrieve a single representative sequence from each cluster, then to translate them, and finally to write them to a multi-fasta file:
p$sequences %>% # Get the sequences
lapply("[[", 1) %>% # Select the first one as representative.
unlist() %>% # Unlist and..
DNAStringSet() %>% # ..transform to DNAStringSet.
translate(if.fuzzy.codon = "solve") %>% # Translate.
writeXStringSet(filepath = "representatives.fasta") # Write a fasta file.
Now the user has to provide the representatives.fasta as input file to the eggnog-mapper. This can be done through the web server, or through the command line if it is installed. The command (on your terminal) would look somewhat like this:
emapper.py -i representative.fasta -o representative_emapper
Then we have to read the *.emapper.annotations file into the R session, and feed the pangenome object with its information:
library(magrittr)
# Read the annotations file
emap <- read.csv("representatives.emapper.annotations",
sep = "\t",
comment.char = "#",
header = FALSE,
na.strings = "-")
# Set column names
colnames(emap) <- c("query", "seed_ortholog", "evalue", "score", "eggNOG_OGs",
"max_annot_lvl", "COG_category", "Description", "Preferred_name",
"GOs", "EC", "KEGG_ko", "KEGG_Pathway", "KEGG_Module", "KE",
"GG_Reaction", "KEGG_rclass", "BRITE", "CAZy",
"BiGG_Reaction", "PFAMs")
# Lets take only some of them which I usually find useful. Subset the data.frame:
cluster_meta <- emap[, c("query", "COG_category", "KEGG_ko", "CAZy")]
# Clean and parse the fields before feeding it to the pangenome
cluster_meta$COG_category <- cluster_meta$COG_category %>% strsplit("")
cluster_meta$KEGG_ko <- clusert_meta$KEGG_ko %>%
gsub("ko:", "", .) %>%
strsplit(",")
cluster_meta$CAZy <- cluster_meta$CAZy %>% strsplit("")
# Add the metadata to the pangenome clusters
p$add_metadata("cluster", cluster_meta)
# Now the object contains the new information in the clusters field:
p$clusters
Neutral genes
The core genome is composed by those genes that are present in every or almost every genome in the sample. As these genes are present in all genomes, they can be used to extract important biological information like phylogenetic relationships, study recombination or selective pressures. For doing these, core genes need to be aligned. In the next few sections we explain how to do this to perform some of the above mentioned downstream analyses.
To reveal the vertical evolutionary history of a bacterial population we should take into account the presence of horizontal acquisition of genetic material, like by means of recombination, that is maintained by the action of natural selection. Hence, we can use the Tajima’s D test of neutrality as implemented in pegas ( Paradis, 2010) to identify those core genes that are not subjected to strong selective pressures and are likely evolving neutrally.
# Load required packages
library(magrittr)
library(DECIPHER)
library(pegas)
library(ape)
p$core_level <- 100 # Set core_level to 100% to avoid
# AlignTranslation() errors.
tajimaD <- p$core_seqs_4_phylo() %>% # Core genome sequences
lapply(DECIPHER::AlignTranslation) %>% # Align translation
lapply(ape::as.DNAbin) %>% # Transform class to DNAbin
lapply(pegas::tajima.test) %>% # Compute Tajima's test
sapply('[[', 'D') # Get Tajima's 'D' statistic from each
# Which are neutral?
which(tajimaD <= 0.2 & tajimaD >= -0.2)
Quick dirty phylogeny
The best practices on how to build a phylogeny from core genomes is a topic of debate. In this first example we provide a one-liner to align individual core genes (like in the previous case), produce a concatenated core genome alignment, calculate a core genome phylogeny using the Neighbor-Joining method as implemented by phangorn (Schliep, 2011), and visualize it using host metadata to color tree tips using ggtree (Yu et al., 2017). Phylogenetic tree is assigned to the phy variable, and plotted as side effect by using magrittr’s %T>% operator.
# Load required packages
library(magrittr)
library(DECIPHER)
library(Biostrings)
library(phangorn)
library(ggtree)
phy <- p$core_seqs_4_phylo() %>% # Core genome sequences
lapply(DECIPHER::AlignSeqs) %>% # Align
do.call(Biostrings::xscat, .) %>% # Concatenate alignments
setNames(p$organisms$org) %>% # Set sequence names
as('matrix') %>% # Transform to matrix
phangorn::phyDat(type = 'DNA') %>% # Transform to phangorn's phyDat
phangorn::dist.ml() %>% # Compute distance
phangorn::NJ() %T>% { # Compute NJ, and assign "phy"
{
ggtree::ggtree(.) %<+% # Create ggtree
as.data.frame(p$organisms) + # Get organisms metadata
ggtree::geom_tippoint(aes(colour = host)) + # Add coloured tip points
scale_color_brewer(palette = 'Set1') # Set color palette
} %>%
print()
}
Maximum Likelihood phylogeny
The following method is similar to the previous, but in this case we are optimizing the topology and branch lengths by maximum likelihood method implemented in phangorn package. Although this is a fully working example, its purpose is just to illustrate the idea. You should consider tuning the parameters to better fit your dataset (e.g.: just 4 discrete gamma distributions for a whole coregenome alignment may be too low).
# Load required packages
library(magrittr)
library(DECIPHER)
library(Biostrings)
library(phangorn)
library(ggtree)
phy <- p$core_seqs_4_phylo() %>% # Core genome sequences
lapply(DECIPHER::AlignSeqs) %>% # Align
do.call(Biostrings::xscat, .) %>% # Concatenate alignments
setNames(p$organisms$org) %>% # Set sequence names
as('matrix') %>% # Transform to matrix
phangorn::phyDat(type = 'DNA') %T>% # Transform to phangorn's phyDat
assign('dat', ., .GlobalEnv) %>% # Assign to "dat" in .GlobalEnv
phangorn::dist.ml() %>% # Compute distance
phangorn::NJ() %>% # Compute NJ (initial tree)
phangorn::pml(data = dat, k = 4) %>% # Compute likelihood with 4 discrete
# gamma distributions.
phangorn::optim.pml(rearrangement = "stochastic", # Optimize likelihood with
# stochastic rearrangements,
optGamma = TRUE, # optimize gamma rate parameter,
optInv = TRUE, # optimize prop of variable size,
model ="GTR") %>% # and use "GTR" model.
magrittr::extract2("tree") %T>% { # Extract the tree only, and pass
{ # it to ggtree.
ggtree::ggtree(.) %<+% # Create ggtree
as.data.frame(p$organisms) + # Get organisms metadata
ggtree::geom_tippoint(aes(colour = host)) + # Add coloured tip points
scale_color_brewer(palette = 'Set1') # Set color palette
} %>%
print()
}
Population structure
Identifying population structure from genomic information is a common problem in microbial ecology. This aims to identify discrete sub-populations within a more heterogeneous population, which is helpful to detect associations with particular phenotypes, geographic origin, host-association, etc. There are various methods for doing this, but probably the most used in microbial genomics is hierBAPS (Cheng et al., 2013) which has been re-implemented in R as rhierBAPS (Tonkin-Hill et al., 2018). In the following recipe we will combine previous examples to 1) align core clusters; 2) compute Tajima’s D statistic over aligned core clusters to identify the ones that are likely to be evolving neutrally; 3) concatenate selected neutral clusters; 4) Run the hierBAPS algorithm; 5) extract lineage information and add it to the pagoo object as organism’s metadata; and 6) compute and plot a tree with lineage information as colour tips.
library(magrittr)
library(DECIPHER)
library(rhierbaps)
library(ape)
library(phangorn)
# 0. Always use core_level at 100% when using DECIPHER::AlignTranslation()
p$core_level <- 100
# 1. Align translation of core genes
ali <- p$core_seqs_4_phylo() %>% # Core genome sequences
lapply(DECIPHER::AlignTranslation) # Align translation
# 2. Identify neutral core clusters
tajD <- ali %>%
lapply(ape::as.DNAbin) %>% # Transform class to DNAbin
lapply(pegas::tajima.test) %>% # Compute Tajima's test
sapply('[[', 'D') # Subset D statistic
neutral <- which(tajD <= 2 & tajD >= -2)
# 3. Concatenate neutral core clusters
concat_neu <- ali[neutral] %>% # Select neutral clusters
do.call(Biostrings::xscat, .) %>% # Concatenate alignments
setNames(p$organisms$org) %>% # Set sequence names
as('matrix') %>% # Transform to matrix
tolower() # Translate to lower case
# 4. Compute structure
rhb <- hierBAPS(snp.matrix = concat_neu, # Input matrix alignment
n.pops = 10, # Max number of subpopulations
max.depth = 1, # Max depth for hierarchical clustering
n.extra.rounds = 5) # Extra rounds to ensure convergence
# 5. Add lineage as metadata to organisms in pagoo object
res <- rhb$partition.df
lin <- data.frame(org = as.character(res[, 1]),
lineage = as.factor(res[, 2]))
p$add_metadata(map = 'org', data = lin)
# 6. Compute tree and plot it with lineage information
concat_neu %>%
phangorn::phyDat(type = 'DNA') %>% # Transform to phangorn's phyDat
phangorn::dist.ml() %>% # Compute distance
phangorn::NJ() %>% # Compute NJ
ggtree::ggtree() %<+% # Create ggtree
as.data.frame(p$organisms) + # Get organisms metadata
ggtree::geom_tippoint(aes(colour = lineage)) # Colour tips with lineage info
Pairwise distance vs pairwise nucleotide diversity
The following recipe computes accessory genome jaccard’s distance and nucleotide diversity of synonym sites for each pair of genomes, and then plot it with a linear regression curve. This method provides insights of accessory genome adaptative evolution. If both magnitudes are not correlated, it could mean that there are some selective pressures shaping clade’s evolution. It first sets the core_level threshold to 100%, then creates a matrix to store pairwise distances and synonym nucleotide diversities. Accessory genome Jaccard distance between each pair of organisms is computed by a pagoo method which is basically a wrapper of vegan::vegdist function. Then there is a long code block with the method to align each core gene and select their polymorphic synonym sites. With this input, pairwise nucleotide diversity is computed with pegas::nuc.div function. At the end, a correlation between both magnitudes is plotted.
library(magrittr)
library(IRanges)
library(Biostrings)
library(DECIPHER)
library(ape)
library(pegas)
# Set core level to 100%. This recipe only works if this is set to 100%.
p$core_level <- 100
# Create pairs matrix
pairs <- data.frame(t(combn(nrow(p$organisms), 2)))
colnames(pairs) <- c('org1', 'org2')
# Compute paired jaccard similarity, transform to matrix
jaccard_sim <- as.matrix(p$dist(method = "jaccard", binary = TRUE))
# Fill results matrix
pairs$jaccard_sim <- apply(pairs, 1, function(i){
ii <- i[1]
jj <- i[2]
jaccard_sim[ii, jj]
})
# Return only synonymous polymorphic sites.
# First, it removes non-synonymous codons, and then retains only
# polymorphyc sites.
syn_poly_sites <- p$core_seqs_4_phylo() %>%
lapply(function(x){
lns <- elementNROWS(x) # Align translatation filtering
tali <- x[which(lns != 0)] %>% # truncated codons and returning
Biostrings::subseq(1L, lns %/% 3 * 3) %>% # both DNA and AA alignments.
DECIPHER::AlignTranslation(type = "both")
syno <- tali[[2]] %>% # Identify non-synonymous
Biostrings::consensusMatrix() %>% # codons.
magrittr::equals(0) %>%
magrittr::not() %>%
colSums() %>%
magrittr::equals(1) %>%
which()
neut <- tali[[1]] %>% # Remove non-synonymous codons.
lapply(function(x){
IRanges::successiveViews(
x, rep.int(3L, length(x) %/% 3L))
}) %>%
lapply('[', syno) %>%
lapply(unlist) %>%
Biostrings::DNAStringSet()
poly <- neut %>% # Identify polymorphic sites.
Biostrings::consensusMatrix() %>%
magrittr::equals(0) %>%
magrittr::not() %>%
colSums() %>%
magrittr::is_greater_than(1) %>%
which()
lapply(neut, '[', poly) %>% # Retain only polymorphic
Biostrings::DNAStringSet() # sites
}) %>%
do.call(Biostrings::xscat, .) %>% # Concatenate.
setNames(p$organisms$org) %>% # Set names.
ape::as.DNAbin() # Convert to DNAbin class.
# Compute paired nucleotide diversity, and fill results matrix
pairs$nuc_div <- apply(pairs, 1, function(i){
ii <- i[1]
jj <- i[2]
pair <- c(syn_poly_sites[ii], syn_poly_sites[jj])
pegas::nuc.div(pair)
})
# Plot correlation with R-base graphics
plot(jaccard_sim ~ nuc_div, pairs)
abline(lm(jaccard_sim ~ nuc_div, pairs))
