In gtonkinhill/panstripe: Pangenome Plots
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>",
fig.path = "inst/vignette-supp/",
echo=TRUE,
warning=FALSE,
message=FALSE,
tidy=TRUE
)
Installation
panstripe is currently available on GitHub. It can be installed with devtools
install.packages('remotes')
remotes::install_github("gtonkinhill/panstripe")
If you would like to also build the vignette with your installation run:
remotes::install_github("gtonkinhill/panstripe", build_vignettes = TRUE)
Quick Start
Panstripe takes as input a phylogeny and gene presence/absence matrix in Rtab format as is produced by Panaroo, Roary, PIRATE and other prokaryotic pangenome tools.
library(panstripe)
library(ape)
library(patchwork)
set.seed(1234)
###
# NOTE: here we load example files from the panstripe package. You should replace these variables with the relevant paths to the files you are using.
phylo.file.name <- system.file("extdata", "tree.newick", package = "panstripe")
rtab.file.name <- system.file("extdata", "gene_presence_absence.Rtab", package = "panstripe")
###
# Load files
pa <- read_rtab(rtab.file.name)
tree <- read.tree(phylo.file.name)
# Run panstripe
fit <- panstripe(pa, tree)
fit$summary
# Plot results
plot_pangenome_params(fit)
plot_pangenome_cumulative(fit)
A significant p-value for the tip term indicates that there is a different rate of gene exchange at the tips of the phylogeny compared with the internal branches. This is usually driven by annotation errors or highly mobile elements that do not persist long enough to be observed in multiple genomes.
A significant p-value for the core term indicates that there is a significant association between the core genome branch length and the number of gene exchange events.
The depth term is less interesting but indicates when there is a difference in our ability to detect older gene exchange events.
Citation
To cite panstripe please use: Tonkin-Hill, G. et al. Robust analysis of prokaryotic pangenome gene gain and loss rates with Panstripe. bioRxiv (2022) doi:10.1101/2022.04.23.489244
Comparing Pangenomes
As panstripe uses a GLM framework it is straightforward to compare the slope terms between datasets. The compare_pangenomes function considers the interaction term between the data sets and the core, tip and depth terms.
IMPORTANT: It is important that the phylogenies for the pangenomes being compared are on the same scale. This can be achieved most easily by using time scaled phylogenies. Alternatively, it is possible to build a single phylogeny and partition it into clades or use SNP scaled phylogenies.
Here, we simulate two pangenomes with different gene exhange rates, fit the model using panstripe and run the compare_pangenomes function.
# Simulate a fast gene gain/loss rate with error
sim_fast <- simulate_pan(rate = 1e-3, ngenomes = 200)
sim_slow <- simulate_pan(rate = 5e-4, ngenomes = 200)
# Run panstripe
fit_fast <- panstripe(sim_fast$pa, sim_fast$tree)
fit_slow <- panstripe(sim_slow$pa, sim_slow$tree)
# Compare the fits
result <- compare_pangenomes(fit_fast, fit_slow)
result$summary
A significant p-value for the tip term indicates that the two pangenomes differ in the rates of gene presence and absence assigned to the tips of the phylogeny. This is usually driven by either differences in the annotation error rates between data sets or differences in the gain and loss of highly mobile elements that do not persist long enough to be observed in multiple genomes.
A significant p-value for the core term indicates that two data sets have different rates of gene gain and loss.
The depth term is less interesting but indicates when there is a difference in our ability to detect older gene exchange events between the two pangenomes.
The dispersion parameter indicates if there is a significant difference in the dispersion of the two pangenomes. This suggests that the relationship between the rate of gene exchange and the size of each event differs in the two pangenomes. Here, the p-value is obtained using a Likelihood Ratio Test.
## Plot the results
plot_pangenome_params(list(fast=fit_fast, slow=fit_slow), legend = FALSE) +
plot_pangenome_cumulative(list(fast=fit_fast, slow=fit_slow)) +
plot_layout(nrow = 1)
Open vs Closed
The definition of what constitutes an open or closed pangenome is somewhat ambiguous. We prefer to consider whether there is evidence for a temporal signal in the pattern of gene gain and loss. This prevents annotation errors leading to misleading results.
After fitting a panstripe model the significance of the temporal signal can be assessed by considering the coefficient and p-value of the 'core' term in the model. The uncertainty of this estimate can also be investigated by looking at the bootstrap confidence intervals of the core term.
Let's simulate a 'closed' pangenome
sim_closed <- simulate_pan(rate = 0, ngenomes = 100)
fit_closed <- panstripe(sim_closed$pa, sim_closed$tree)
fit_closed$summary
The p-value indicates that there is not a significant association between core branch lengths and gene gain/loss. This is typical of species that undergo very little to no recombination such as M. tuberculosis.
Rate vs Size
While comparing the core parameter of the model identifies differences in the association between branch lengths and gene gain and loss it does not indicate whether this is driven by higher rates of recombination or simply larger recombination events involving more genes.
The panstripe model allows these two scenarios to be investigated by allowing the dispersion parameter in the model to be different for each pangenome.
Here, we simulate two data sets with the same recombination rate but where each recombination event differs in the number of genes involved.
sim_large <- simulate_pan(rate = 1e-3, ngenomes = 100, mean_trans_size = 4)
sim_small <- simulate_pan(rate = 1e-3, ngenomes = 100, mean_trans_size = 3)
fit_large <- panstripe(sim_large$pa, sim_large$tree)
fit_small <- panstripe(sim_small$pa, sim_small$tree)
# Compare the fits
result <- compare_pangenomes(fit_large, fit_small)
result$summary
Output
The panstripe function generates a list with the following attributes
summary
A table indicating a subset of the inferred parameters of the GLM. The p-values and bootstrap confidence intervals can be used to determine whether each term in the model is significantly associated with gene gain and loss.
-
core indicates whether the branch lengths in the phylogeny are associated with gene gain and loss.
-
tip indicates associations with genes observed to occur on the tips of the phylogeny. These are usually driven by a combination of annotation errors and depending upon the temporal sampling density also highly mobile elements that are not observed in multiple genomes.
-
depth indicates whether the rate of gene gain and loss changes significantly with the depth of a branch. Typically, our ability to detect gene exchange events reduces for older ancestral branches.
-
p the inferred index parameter of the underlying Tweedie distribution used in the GLM
-
phi the inferred dispersion parameter of the GLM
model
The output of fitting the GLM model. This object can be used to predict how many gene gains and losses we would expect to observe given the parameters of a particular branch.
data
A table with the data used to fit the model. The acc column indicates the inferred number of gene gain/loss events for a branch; the core column is the branch length taken from the phylogeny; the istip column indicates whether the branch occurs at the tip of the phylogeny and the depth column indicates the distance from the root node to the branch.
ci_samples
A 'boot' object, generated by the boot package. Can be used to investigate the uncertainty in the parameter estimates.
tree/pa
The original data provided to the panstripe function.
Alternative models
The default Panstripe model assumes a Compound Poisson (Tweedie) distribution as implemented in the Tweedie R package. In some cases there may be insufficient data or the data may not reliably fit the Tweedie distribution. This is often the case when there is nearly no gene exchange events inferred to have occurred at the internal branches of the phylogeny.
To help account for these issues and to add flexibility to the package it is also possible to fit alternative distributions in place of the Tweedie model. The most common alternative that we suggest is to use a Gaussian distribution. This is usually more robust and is less likely to run into convergence issues in the model fit. In our tests both of these distributions generally give very similar results.
Panstripe can be run using an alternative distribution as
fit_gaussian <- panstripe(pa, tree, family='gaussian')
Plots
Panstripe includes a number of useful plotting functions to help with interpretation of the output of panaroo.
Pangenome fit
The inferred parameters of the Panstripe model fit can be plotted as
plot_pangenome_params(fit_fast)
We also suggest plotting the cumulative number of gene gain and loss events against the core branch length. As Panstripe models each branch individually and accounts for the depth of the branch it is not possible to add the Panstripe model fit to this plot.
plot_pangenome_cumulative(fit_fast)
These two plots can easily be combined using the patchwork R package as
plot_pangenome_params(fit_fast) +
plot_pangenome_cumulative(fit_fast) +
plot_layout(nrow = 1)
The functions can also take a named list as input allowing for easy comparisons between data sets.
plot_pangenome_params(list(fast=fit_fast, slow=fit_slow), legend = FALSE) +
plot_pangenome_cumulative(list(fast=fit_fast, slow=fit_slow)) +
plot_layout(nrow = 1)
A plot of the residuals of the regression can also be generated using
plot_residuals(fit_fast)
Tree with presence/absence
A plot of the phylogeny and the corresponding gene presence/absence matrix can be made using the plot_tree_pa function. This function also takes an optional vector of gene names to include in the plot.
#look at only those genes that vary
variable_genes <- colnames(pa)[apply(pa, 2, sd)>0]
plot_tree_pa(tree = tree, pa = pa, genes = variable_genes, label_genes = FALSE, cols = 'black')
Inferred ancestral states
The plot_gain_loss function allows for the visualisation of the fitted gene gain and loss events on the given phylogeny. The enrichment for events at the tips of a tree is often driven by a combination of highly mobile elements and annotation errors.
plot_gain_loss(fit)
tSNE
The tSNE dimension reduction technique can be used to investigate evidence for clusters within the pangenome.
plot_tsne(pa)
The Mandrake method can also be used as an alternative to tSNE.
Accumulation curves
While we do not recommend the use of accumulation curves as they do not account for population structure, sampling bias or annotation errors we have included a function to plot them to make it easier for users to compare methods.
plot_acc(list(fast=sim_fast$pa, slow=sim_slow$pa))
Etymology
The name panstripe is an adaptation of 'pinstripe', the name of a long-nosed potoroo who was a villain in the playstation game Crash Bandicoot. The name was chosen as the program relies on the output of panaroo (named after the potoroo).
gtonkinhill/panstripe documentation built on Feb. 27, 2025, 9:01 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
knitr::opts_chunk$set( collapse = TRUE, comment = "#>", fig.path = "inst/vignette-supp/", echo=TRUE, warning=FALSE, message=FALSE, tidy=TRUE )
Installation
panstripe is currently available on GitHub. It can be installed with devtools
install.packages('remotes') remotes::install_github("gtonkinhill/panstripe")
If you would like to also build the vignette with your installation run:
remotes::install_github("gtonkinhill/panstripe", build_vignettes = TRUE)
Quick Start
Panstripe takes as input a phylogeny and gene presence/absence matrix in Rtab format as is produced by Panaroo, Roary, PIRATE and other prokaryotic pangenome tools.
library(panstripe) library(ape) library(patchwork) set.seed(1234) ### # NOTE: here we load example files from the panstripe package. You should replace these variables with the relevant paths to the files you are using. phylo.file.name <- system.file("extdata", "tree.newick", package = "panstripe") rtab.file.name <- system.file("extdata", "gene_presence_absence.Rtab", package = "panstripe") ### # Load files pa <- read_rtab(rtab.file.name) tree <- read.tree(phylo.file.name) # Run panstripe fit <- panstripe(pa, tree) fit$summary # Plot results plot_pangenome_params(fit) plot_pangenome_cumulative(fit)
A significant p-value for the tip term indicates that there is a different rate of gene exchange at the tips of the phylogeny compared with the internal branches. This is usually driven by annotation errors or highly mobile elements that do not persist long enough to be observed in multiple genomes.
A significant p-value for the core term indicates that there is a significant association between the core genome branch length and the number of gene exchange events.
The depth term is less interesting but indicates when there is a difference in our ability to detect older gene exchange events.
Citation
To cite panstripe please use: Tonkin-Hill, G. et al. Robust analysis of prokaryotic pangenome gene gain and loss rates with Panstripe. bioRxiv (2022) doi:10.1101/2022.04.23.489244
Comparing Pangenomes
As panstripe uses a GLM framework it is straightforward to compare the slope terms between datasets. The compare_pangenomes function considers the interaction term between the data sets and the core, tip and depth terms.
IMPORTANT: It is important that the phylogenies for the pangenomes being compared are on the same scale. This can be achieved most easily by using time scaled phylogenies. Alternatively, it is possible to build a single phylogeny and partition it into clades or use SNP scaled phylogenies.
Here, we simulate two pangenomes with different gene exhange rates, fit the model using panstripe and run the compare_pangenomes function.
# Simulate a fast gene gain/loss rate with error sim_fast <- simulate_pan(rate = 1e-3, ngenomes = 200) sim_slow <- simulate_pan(rate = 5e-4, ngenomes = 200) # Run panstripe fit_fast <- panstripe(sim_fast$pa, sim_fast$tree) fit_slow <- panstripe(sim_slow$pa, sim_slow$tree) # Compare the fits result <- compare_pangenomes(fit_fast, fit_slow) result$summary
A significant p-value for the tip term indicates that the two pangenomes differ in the rates of gene presence and absence assigned to the tips of the phylogeny. This is usually driven by either differences in the annotation error rates between data sets or differences in the gain and loss of highly mobile elements that do not persist long enough to be observed in multiple genomes.
A significant p-value for the core term indicates that two data sets have different rates of gene gain and loss.
The depth term is less interesting but indicates when there is a difference in our ability to detect older gene exchange events between the two pangenomes.
The dispersion parameter indicates if there is a significant difference in the dispersion of the two pangenomes. This suggests that the relationship between the rate of gene exchange and the size of each event differs in the two pangenomes. Here, the p-value is obtained using a Likelihood Ratio Test.
## Plot the results plot_pangenome_params(list(fast=fit_fast, slow=fit_slow), legend = FALSE) + plot_pangenome_cumulative(list(fast=fit_fast, slow=fit_slow)) + plot_layout(nrow = 1)
Open vs Closed
The definition of what constitutes an open or closed pangenome is somewhat ambiguous. We prefer to consider whether there is evidence for a temporal signal in the pattern of gene gain and loss. This prevents annotation errors leading to misleading results.
After fitting a panstripe model the significance of the temporal signal can be assessed by considering the coefficient and p-value of the 'core' term in the model. The uncertainty of this estimate can also be investigated by looking at the bootstrap confidence intervals of the core term.
Let's simulate a 'closed' pangenome
sim_closed <- simulate_pan(rate = 0, ngenomes = 100) fit_closed <- panstripe(sim_closed$pa, sim_closed$tree) fit_closed$summary
The p-value indicates that there is not a significant association between core branch lengths and gene gain/loss. This is typical of species that undergo very little to no recombination such as M. tuberculosis.
Rate vs Size
While comparing the core parameter of the model identifies differences in the association between branch lengths and gene gain and loss it does not indicate whether this is driven by higher rates of recombination or simply larger recombination events involving more genes.
The panstripe model allows these two scenarios to be investigated by allowing the dispersion parameter in the model to be different for each pangenome.
Here, we simulate two data sets with the same recombination rate but where each recombination event differs in the number of genes involved.
sim_large <- simulate_pan(rate = 1e-3, ngenomes = 100, mean_trans_size = 4) sim_small <- simulate_pan(rate = 1e-3, ngenomes = 100, mean_trans_size = 3) fit_large <- panstripe(sim_large$pa, sim_large$tree) fit_small <- panstripe(sim_small$pa, sim_small$tree) # Compare the fits result <- compare_pangenomes(fit_large, fit_small) result$summary
Output
The panstripe function generates a list with the following attributes
summary
A table indicating a subset of the inferred parameters of the GLM. The p-values and bootstrap confidence intervals can be used to determine whether each term in the model is significantly associated with gene gain and loss.
-
core indicates whether the branch lengths in the phylogeny are associated with gene gain and loss.
-
tip indicates associations with genes observed to occur on the tips of the phylogeny. These are usually driven by a combination of annotation errors and depending upon the temporal sampling density also highly mobile elements that are not observed in multiple genomes.
-
depth indicates whether the rate of gene gain and loss changes significantly with the depth of a branch. Typically, our ability to detect gene exchange events reduces for older ancestral branches.
-
p the inferred index parameter of the underlying Tweedie distribution used in the GLM
-
phi the inferred dispersion parameter of the GLM
model
The output of fitting the GLM model. This object can be used to predict how many gene gains and losses we would expect to observe given the parameters of a particular branch.
data
A table with the data used to fit the model. The acc column indicates the inferred number of gene gain/loss events for a branch; the core column is the branch length taken from the phylogeny; the istip column indicates whether the branch occurs at the tip of the phylogeny and the depth column indicates the distance from the root node to the branch.
ci_samples
A 'boot' object, generated by the boot package. Can be used to investigate the uncertainty in the parameter estimates.
tree/pa
The original data provided to the panstripe function.
Alternative models
The default Panstripe model assumes a Compound Poisson (Tweedie) distribution as implemented in the Tweedie R package. In some cases there may be insufficient data or the data may not reliably fit the Tweedie distribution. This is often the case when there is nearly no gene exchange events inferred to have occurred at the internal branches of the phylogeny.
To help account for these issues and to add flexibility to the package it is also possible to fit alternative distributions in place of the Tweedie model. The most common alternative that we suggest is to use a Gaussian distribution. This is usually more robust and is less likely to run into convergence issues in the model fit. In our tests both of these distributions generally give very similar results.
Panstripe can be run using an alternative distribution as
fit_gaussian <- panstripe(pa, tree, family='gaussian')
Plots
Panstripe includes a number of useful plotting functions to help with interpretation of the output of panaroo.
Pangenome fit
The inferred parameters of the Panstripe model fit can be plotted as
plot_pangenome_params(fit_fast)
We also suggest plotting the cumulative number of gene gain and loss events against the core branch length. As Panstripe models each branch individually and accounts for the depth of the branch it is not possible to add the Panstripe model fit to this plot.
plot_pangenome_cumulative(fit_fast)
These two plots can easily be combined using the patchwork R package as
plot_pangenome_params(fit_fast) + plot_pangenome_cumulative(fit_fast) + plot_layout(nrow = 1)
The functions can also take a named list as input allowing for easy comparisons between data sets.
plot_pangenome_params(list(fast=fit_fast, slow=fit_slow), legend = FALSE) + plot_pangenome_cumulative(list(fast=fit_fast, slow=fit_slow)) + plot_layout(nrow = 1)
A plot of the residuals of the regression can also be generated using
plot_residuals(fit_fast)
Tree with presence/absence
A plot of the phylogeny and the corresponding gene presence/absence matrix can be made using the plot_tree_pa function. This function also takes an optional vector of gene names to include in the plot.
#look at only those genes that vary variable_genes <- colnames(pa)[apply(pa, 2, sd)>0] plot_tree_pa(tree = tree, pa = pa, genes = variable_genes, label_genes = FALSE, cols = 'black')
Inferred ancestral states
The plot_gain_loss function allows for the visualisation of the fitted gene gain and loss events on the given phylogeny. The enrichment for events at the tips of a tree is often driven by a combination of highly mobile elements and annotation errors.
plot_gain_loss(fit)
tSNE
The tSNE dimension reduction technique can be used to investigate evidence for clusters within the pangenome.
plot_tsne(pa)
The Mandrake method can also be used as an alternative to tSNE.
Accumulation curves
While we do not recommend the use of accumulation curves as they do not account for population structure, sampling bias or annotation errors we have included a function to plot them to make it easier for users to compare methods.
plot_acc(list(fast=sim_fast$pa, slow=sim_slow$pa))
Etymology
The name panstripe is an adaptation of 'pinstripe', the name of a long-nosed potoroo who was a villain in the playstation game Crash Bandicoot. The name was chosen as the program relies on the output of panaroo (named after the potoroo).
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.