Easy debugging of C/C++/CUDA python extensions

Writing an extension called by python (in C, C++ or CUDA)? Not working? Typical.

When doing the same from R it’s pretty easy to debug, just run with R -d <debugger name> e.g. R -d valgrind or R -d gdb you get into the debugger, continue, then run interactively as usual. (For a more complex example using both at once see this blog post).

Doing this from python seems trickier to me. I started off following this guide: https://johnfoster.pge.utexas.edu/blog/posts/debugging-cc%2B%2B-libraries-called-by-python/
But I’m not really an ipython user, and prefer gdb over lldb (due to familiarity). I think this is a good way to do this if you need an interactive python session, but really this is overcomplicated for my typical use case.

Now all I do is run (say I’m debugging CUDA code):
cuda-gdb python

Then in the debugger (an example with pp-sketchlib):

set args ~/installs/pp-sketchlib/pp_sketch-runner.py --sketch --rfile rfile_test.txt --ref-db Pf6 --cpus 8 --use-gpu

Works just fine! (as long as you remember to compile with the debug flags on)

Things I have learnt about using R (as a python/C++ programmer)

  • For loops are ok. I thought these were impossibly slow in R, and you always had to vectorise code – but not so. A common error is to try and append to arrays. It always best to preallocate an array/list with result <- rep(NA, length) or similar, and then write to it by index.
  • strings aren’t the easiest thing to deal with, but stringr might help, and there is a regular expression package built-in. Use paste0 to concatenate.
  • That being said, if you can easily vectorise code you should do so.
  • It’s really easy to set a function as a variable and pass these around, bind its arguments (a partial in python) etc. Use this feature!
  • A matrix type has to contain the same data type, so if you have different data types and convert to a matrix, they will all be converted to a compatible type silently.
  • Use ‘%*%’ for matrix multiplication.
  • Use ‘% in %’ for sets.
  • Use array[array$column == value, ] or similar, for selecting values from a data.frame.
  • Use a single ‘&’ or ‘|’ when combining conditions in the above. Double ‘&&’ will short circuit.
  • data.frame columns must have the same type. If you want to mix types you’ll likely need a list.
  • Don’t use apply, as the first step is to convert to a matrix, making all your types the same. Use vapply, which defines the expected return type
  • R will use functions which are partial matches to the name you called without commenting on this behaviour (!?!). Add ‘options(warnPartialMatchAttr=TRUE,
    warnPartialMatchArgs=TRUE)’ to ~/.Rprofile to turn this strange default off.
  • A numeric (float) and an integer are different types. Use 1L to make an integer ‘1’.
  • furrr is a nice library for parallelisation.
  • In general, tidyverse packages offer good alternatives to many data science functions in base R.
  • Extremely basic OOP is available using ‘S3’ objects (pretty much inheritance, and overriding of some typical functions such as print, summary based on type). ‘S4’ is to be avoided, apparently. ‘R6’ gives you more typical features.
  • Use the devtools package for your development, it automates most building and testing of the code.
  • RStudio has a nice built-in profiler.
  • R has a good FFI with C and can automatically sort out compiling for you. C++ is also possible with RCpp, but I believe a little more involved.

Paper summary – PopPUNK for bacterial epidemiology

A paper describing our recent method for bacterial epidemiology PopPUNK has just been published in Genome Research, which you can read here:

You can install our software by running
conda install poppunk
and that full details and documentation can be found at https://poppunk.readthedocs.io

In this blog post I will attempt to describe some of our key features and findings in a shorter format. Broadly, I think there are three main parts:

  1. Core and accessory distances can be estimated using k-mer distances.
  2. In many species, finding clusters in core-accessory space gives good quality clusters of strains. The core and accessory distances from 1) provide further resolution within clusters of strains.
  3. These clusters can be expanded online as new samples are added, and their names stay consistent between studies.

I’ve also noted some of the work we added in our revision, for those that might have seen the first version as a pre-print on bioRxiv. We added more direct comparisons with phylogenies and cg/wgMLST schemes, showing that PopPUNK was preferable to wg/cgMLST, while still fulfilling the criteria desirable for an epidemiological typing system laid out by Nadon et al.

Core and accessory distances can be estimated using k-mer distances

The importance of accessory genome evolution and divergence has been increasingly recognised over the past few years. To analyse the accessory genome, one typically attempts to find clusters of orthologous genes (COGs) using roary, panX or another similar method. These methods compare all annotated genes to all others, which results in a number of comparisons which increases with the square of the number of sequences. Though efficiencies in these pieces of software keep this computation possible, for larger populations this takes a significant amount of time, especially if reruns are needed due to new samples or poorly chosen clustering thresholds.

For some downstream purposes just extracting the core and accessory distances between pairs of samples is sufficient, as information on individual COGs and annotations is not needed. We wanted to use a k-mer based approach to do this, so that we:

  • wouldn’t be reliant on gene annotations, which differ in quality between species and sample collection, and can be hard to standardise across labs and studies.
  • do not require an alignment step, which takes a long time, and again can be hard to standardise.
  • could take advantage of recent developments in efficient k-mer based software. Specifically, we used mash, which rapidly calculates distances between sequences by taking the lowest scoring subset (a sketch) of k-mers from a sequence assembly.

Noting that longer k-mers are more likely to mismatch between samples due to SNP in a shared (core) region, but that k-mers of all sizes are equally likely to mismatch between samples due to a missing accessory element (longer than the k-mer length), we were able to formulate a relation between mash distances at various k-mer lengths and core and accessory distances. Ultimately, this allows us to calculate core and accessory distances between all sequences in a population tens or hundreds of times faster than from clustering and aligning genes. In a population of 128 Listeria monocytogenes PopPUNK took about ten minutes, whereas a run of roary alone took 31 hrs. We also compared our results to this method in simulations (figure 2 in the paper) and in ten varied species (figure 3 in the paper) and found our faster estimates to be consistent with both our simulations and the real data.

We can then plot the core and accessory distances for all pairwise comparisons of samples, adding density countours where many points overlap. Here is the L. monocytogenes example:

L. monocytogenes distance distribution
Core and accessory distances between every pair of 128 L. monocytogenes samples, with density contours

This distribution is useful for a number of things, particularly clustering – the focus of the rest of the paper – but can also tell us about overall core-accessory evolution, and can be used to pick out samples which have unexpected divergence in either core or accessory content (see supplementary figure 11 for a detailed example of this).

Using core-accessory distances to define sequence clusters (strains)

A good, widely used method to define clusters of closely related sequences in a population is hierBAPS, or the recent upgrade fastbaps (both fit the same model, but the newer version is significantly faster at doing so and can also use a phylogeny to constrain the possible clusters). While this approach has many nice features such as being able to cluster recursively and being able to extract likelihoods for fits and assignments, the following limitations make it challenging to directly apply in all the places where subtyping of a population is useful for epidemiology:

  • Requires a good quality alignment as input (time-consuming).
  • Is difficult to update (every new genome requires a complete refit).
  • Difficult to keep cluster names consistent (different studies and runs will have different cluster IDs, and possibly assignments).
  • Usually forms a bin cluster of outlier sequences (putting all low frequency clusters together, rather than separating them).
  • Can be slow to fit (though perhaps no longer, with the development of fastbaps).

These drawbacks make a species-level definition of subtypes potentially challenging. Additionally, for some species (of particular interest to us was the Streptococcus genus) the solution found by optimising the BAPS likelihood does not provide great quality clusters across the tree, I think due to unmodelled recombination events and many small clusters.

This is what we set out to try and improve with PopPUNK, hoping that our fast estimation of core and accessory distances could be used for this purpose.

By finding clusters that are clearly separated in core-accessory space (using one of two standard machine-learning methods) we are able to determine a cutoff for which distances are within the same strain. Applying this to the same distances as above:

GMM fit to L. monocytogenes
Four component Gaussian mixture model fitted to pairwise distances in L. monocytogenes

The light-blue cluster closest to the origin is the within-strain cluster – distances in this cluster represent comparisons between samples in the same strain. We can then draw links between any pair of samples less than this distance apart. Linked samples form the clusters of strains in a network, the connected components (see figure 5A,B). In the network samples A and B may be greater than the cutoff distance apart, if both are close enough to a third sample C they will be in the same cluster. For most of the species we applied our method to this approach gave good clusters very quickly (table 1, figure 4). For two Streptococcal species where extensive recombination blurred the separation between the components, we needed to apply a final step to adjust the position of the boundary. Optimising properties of the network, avoiding clusters which are straggly and linked by only a few samples connected to many things, and reducing the overall number of links, then gave good results in these species without further extensive computation.

Adding new isolates in to a fitted cluster model – quickly and consistently

We found some very useful advantages to representing the clusters as a network, which solve many of the above issues:

  • Outlier sequences or small groups are correctly placed in their own clusters, rather than being grouped in a bin.
  • New isolates can be added in by calculating distances to samples already in the network (avoiding calculating everything again).
  • Removing fully connected sets of samples (cliques) removes redundancy, and further improves speed/memory/storage.
  • Cluster names are by size, and can be constrained to be the same between studies, or as a database is added to.

This means that you can download a PopPUNK database (usually 10-100Mb) and run using --assign-samples with new assemblies. This will cluster new samples within the context of an existing population, without having to redo/care about the model fit. The databases can be expanded without having to refit the model, or worry about cluster names changing (which is one of the nice features of MLST). We tested this with an emerging E. coli strain not seen in the database at the first time point in a longitudinal series, and PopPUNK was able to track its emergence and expansion (see figure 5D,E).

Comparison with gene-by-gene methods

For the second version of this paper we were asked to add in a more explicit comparison with gene-by-gene methods and phylogenies. My understanding of how MLST and cgMLST/wgMLST schemes are applied in epidemiology is:

  1. Download the MLST allele database, or upload your sample to an online database.
  2. BLAST (or another alignment/matching tool) is used to find the typing genes in the uploaded sample.
  3. These genes are compared to a list of previously seen sequences for that gene (alleles). If they exactly match they are assigned the same identifier, otherwise a new one.
  4. A sequence type (ST) is assigned to each unique combination of alleles in the typed genes.

In step 3, counting any number of changes within a gene as a single change loses some resolution, but has the advantage that it does not overcount recombination events. With a good choice of genes making up the scheme, MLST schemes have been shown to capture population structure very well. It is faster than alignment and modelling with hierBAPS, a single sample can easily be added, and with centralised databases it can also deal with keeping names of clusters (STs) consistent.

However, some drawbacks are:

  • MLST lacks resolution, and throws data from whole genome sequencing away. Extended schemes (using core genes or all genes) improve this, but still ignore intergenic regions and group any number of changes together as a single allele difference.
  • Defining allele databases is difficult and requires continual curation for people to add new alleles, and avoid duplicate genes. There are some impressive efforts to do this (e.g. enterobase, BIGSdb), but these only cover some species.
  • It is unclear how to form clusters from allele assignment. How many changes of allele should be in the same cluster? Common interpretations of a single change is far too specific for extended schemes, and seems to be a convenience-based choice rather than a biological one.

In practice, I found that downloading a cgMLST scheme and applying it to my own data was quite challenging due to how the gene database needed to be formatted, and to make sure all the dependencies worked (thanks to João Carriço and Mickael Silva for helping me with this, and for their chewBBACA software which made the comparison possible). MLST methods and databases have been around for longer, and so this was easier to work with. Defining and maintaining a new scheme for a species which doesn’t have one yet seems like it would be a significant undertaking, though I didn’t try this myself.

To directly compare PopPUNK and these methods, we performed MLST and cgMLST assignment on two different species with good typing schemes (L. monocytogenes and E. coli). We then calculated pairwise distances in terms of the number of allele changes, which gave a gene distance matrix rather than a core and accessory distance matrix. By using these with the PopPUNK network we could find how many allele changes to connect to form similar clusters, and how good projections at various cutoffs are (see supplementary tables 6 and 7). We could make the clusters similar between PopPUNK and (cg)MLST, but only by manually testing many values of the cutoff for number of allele changes.

I don’t have lots of experience using gene-by-gene methods or analysing surveillance datasets, but from these tests I ended up concluding that PopPUNK has the following advantages over gene-by-gene methods:

  • The clusters are likely to be real biological entities based on their separation in core-accessory space. With cgMLST/MLST it is unclear how many allele changes should be chosen as a cutoff, but PopPUNK optimises this.
  • The calculation of core distances allow trees to be drawn within clusters (e.g. with the --microreact option), giving further resolution and relationships within clusters.
  • PopPUNK is faster to run.
  • You don’t need curated database of genes (difficult to format/curate when it exists, more difficult to create when it doesn’t).

So we ended up concluding that PopPUNK also retains the advantages of gene-by-gene approaches, and meets the criteria of Nadon et al for a genomic surveillance scheme.

Challenging species – low diversity

The main place we found PopPUNK’s clusters to be worse than those from RhierBAPS was for populations with limited genetic diversity, for example within an identified strain. The calculation of core and accessory distances will in theory work to any resolution (but one may need to increase the sketch size to the genome length divided by the variants per genome). But if there is no clear within-strain versus between-strain separation in the distances and instead just a cloud of points, the spatial clustering methods are not likely to converge on a good solution. Network-based model refinement is needed in this case, though it is likely to split the strain into many substrains.

One example of this was Neisseria gonorrhoeae, which is essentially a strain of Neisseria meningitidis (which did work using default settings). Using refinement of core distances we did get a reasonable fit, and were able to use this to find accessory elements moving within strains (we also looked at this within a well studied strain of S. pneumoniae). In Mycobacterium tuberculosis diversity was even more limited, so while the core distance based phylogeny PopPUNK produced was consistent with the lineages estimated by the first level of hierBAPS, PopPUNK’s clustering split the population into many more substrains (comparable to spoligotyping).

See the supplementary text S1 and figure S11 for a full discussion of this point.

Recommendations for running PopPUNK

Some tips:

  • The most important part of the model fit is to uniquely select the cluster closest to the origin.
  • If you have high accessory distance points check for contamination in assemblies and remove them.
  • Run using one of the extra output options (e.g. --microreact) and you’ll get a lot more information out, and can make interactive visualisations.

See the quickstart guide for a walkthrough, the tutorial for all details on a more complex example and the troubleshooting for common issues.


PopPUNK was the result of a collaboration between many people, but I’d particularly like to thank Nick Croucher who jointly worked on the method, code and paper with me.

conda build: libarchive: cannot open shared object file: No such file or directory

I was getting the following error, attempting to run conda-build on a package, using a conda env:

Traceback (most recent call last):
File "/nfs/users/nfs_j/jl11/pathogen_nfs/large_installations/miniconda3/envs/conda_py36/bin/conda-build", line 7, in <module>
from conda_build.cli.main_build import main
File "/nfs/users/nfs_j/jl11/pathogen_nfs/large_installations/miniconda3/envs/conda_py36/lib/python3.6/site-packages/conda_build/cli/main_build.py", line 18, in <module>
import conda_build.api as api
File "/nfs/users/nfs_j/jl11/pathogen_nfs/large_installations/miniconda3/envs/conda_py36/lib/python3.6/site-packages/conda_build/api.py", line 22, in <module>
from conda_build.config import Config, get_or_merge_config, DEFAULT_PREFIX_LENGTH as _prefix_length
File "/nfs/users/nfs_j/jl11/pathogen_nfs/large_installations/miniconda3/envs/conda_py36/lib/python3.6/site-packages/conda_build/config.py", line 17, in <module>
from .variants import get_default_variant
File "/nfs/users/nfs_j/jl11/pathogen_nfs/large_installations/miniconda3/envs/conda_py36/lib/python3.6/site-packages/conda_build/variants.py", line 15, in <module>
from conda_build.utils import ensure_list, trim_empty_keys, get_logger
File "/nfs/users/nfs_j/jl11/pathogen_nfs/large_installations/miniconda3/envs/conda_py36/lib/python3.6/site-packages/conda_build/utils.py", line 10, in <module>
import libarchive
File "/nfs/users/nfs_j/jl11/pathogen_nfs/large_installations/miniconda3/envs/conda_py36/lib/python3.6/site-packages/libarchive/__init__.py", line 1, in <module>
from .entry import ArchiveEntry
File "/nfs/users/nfs_j/jl11/pathogen_nfs/large_installations/miniconda3/envs/conda_py36/lib/python3.6/site-packages/libarchive/entry.py", line 6, in <module>
from . import ffi
File "/nfs/users/nfs_j/jl11/pathogen_nfs/large_installations/miniconda3/envs/conda_py36/lib/python3.6/site-packages/libarchive/ffi.py", line 27, in <module>
libarchive = ctypes.cdll.LoadLibrary(libarchive_path)
File "/nfs/users/nfs_j/jl11/pathogen_nfs/large_installations/miniconda3/envs/conda_py36/lib/python3.6/ctypes/__init__.py", line 426, in LoadLibrary
return self._dlltype(name)
File "/nfs/users/nfs_j/jl11/pathogen_nfs/large_installations/miniconda3/envs/conda_py36/lib/python3.6/ctypes/__init__.py", line 348, in __init__
self._handle = _dlopen(self._name, mode)
OSError: /nfs/users/nfs_j/jl11/pathogen_nfs/large_installations/miniconda3/envs/conda_py36/lib/libarchive: cannot open shared object file: No such file or directory

The /lib directory does contain libarchive, both as a dynamic (.so) and static (.a) library. There turned out to be two relevant environment variables:


A workaround is then to run

export LIBARCHIVE=/nfs/users/nfs_j/jl11/pathogen_nfs/large_installations/miniconda3/envs/conda_py36/lib/libarchive.so

There is probably a proper reason why this has happened and a permanent solution to the issue, but this works for now.

Firth regression in python

Marco Galardini and I have recently reimplemented the bacterial GWAS software SEER in python. As part of this I rewrote my C++ code for Firth regression in python. Firth regression gives better estimates when data in logistic regression is separable or close to separable (when a chi-squared contingency table has small entries).

I found that although there is an R implementation logistf I couldn’t find an equivalent in another language, or python’s statsmodels. Here is a gist with my python functions and a skeleton of how to use them and calculate p-values, in case anyone would like to use this in future without having to write the optimiser themselves.

#!/usr/bin/env python
'''Python implementation of Firth regression by John Lees
See https://www.ncbi.nlm.nih.gov/pubmed/12758140'''
def firth_likelihood(beta, logit):
return (logit.loglike(beta) + 0.5*np.log(np.linalg.det(logit.hessian(beta))))
# Do firth regression
# Note information = -hessian, for some reason available but not implemented in statsmodels
def fit_firth(y, X, start_vec=None, step_limit=1000, convergence_limit=0.0001):
logit_model = smf.Logit(y, X)
if start_vec is None:
start_vec = np.zeros(X.shape[1])
beta_iterations = []
for i in range(0, step_limit):
pi = logit_model.predict(beta_iterations[i])
W = np.diagflat(np.multiply(pi, 1pi))
var_covar_mat = np.linalg.pinv(logit_model.hessian(beta_iterations[i]))
# build hat matrix
rootW = np.sqrt(W)
H = np.dot(np.transpose(X), np.transpose(rootW))
H = np.matmul(var_covar_mat, H)
H = np.matmul(np.dot(rootW, X), H)
# penalised score
U = np.matmul(np.transpose(X), y pi + np.multiply(np.diagonal(H), 0.5 pi))
new_beta = beta_iterations[i] + np.matmul(var_covar_mat, U)
# step halving
j = 0
while firth_likelihood(new_beta, logit_model) > firth_likelihood(beta_iterations[i], logit_model):
new_beta = beta_iterations[i] + 0.5*(new_beta beta_iterations[i])
j = j + 1
if (j > step_limit):
sys.stderr.write('Firth regression failed\n')
return None
if i > 0 and (np.linalg.norm(beta_iterations[i] beta_iterations[i1]) < convergence_limit):
return_fit = None
if np.linalg.norm(beta_iterations[i] beta_iterations[i1]) >= convergence_limit:
sys.stderr.write('Firth regression failed\n')
# Calculate stats
fitll = firth_likelihood(beta_iterations[1], logit_model)
intercept = beta_iterations[1][0]
beta = beta_iterations[1][1:].tolist()
bse = np.sqrt(np.diagonal(np.linalg.pinv(logit_model.hessian(beta_iterations[1]))))
return_fit = intercept, beta, bse, fitll
return return_fit
if __name__ == "__main__":
import sys
import warnings
import math
import statsmodels
import numpy as np
from scipy import stats
import statsmodels.api as smf
# create X and y here. Make sure X has an intercept term (column of ones)
# …
# How to call and calculate p-values
(intercept, beta, bse, fitll) = fit_firth(y, X)
beta = [intercept] + beta
# Wald test
waldp = []
for beta_val, bse_val in zip(beta, bse):
waldp.append(2 * (1 stats.norm.cdf(abs(beta_val/bse_val))))
lrtp = []
for beta_idx, (beta_val, bse_val) in enumerate(zip(beta, bse)):
null_X = np.delete(X, beta_idx, axis=1)
(null_intercept, null_beta, null_bse, null_fitll) = fit_firth(y, null_X)
lrstat = 2*(null_fitll fitll)
lrt_pvalue = 1
if lrstat > 0: # non-convergence
lrt_pvalue = stats.chi2.sf(lrstat, 1)