Using recursive CTEs for calculating taxonomic lineages

This topic falls into the category of "things that I didn't know that I didn't know." Many of my projects involve the manipulation of taxonomies, particularly the NCBI taxonomy, but I didn't understand the efficiency with which it was possible to calculate a lineage in SQL. I recently revisited the issue, and here's what I learned.

The NCBI taxonomy is provided as a collection of tables that can be downloaded in its entirety by FTP. There are two tables in the taxonomy database that we care about for the purposes of this post. The first is nodes. The table from the NCBI taxonomy contains some additional fields that aren't important in this context, but at its core, it looks something like this:

create table nodes (
       tax_id integer,
       parent integer,
       rank text
);

The second table that we will use later is called names, with a schema that (in its simplified form) looks like this:

create table names (
       tax_id integer,
       tax_name text,
       is_primary boolean
);

Multiple synonyms are allowed for each tax_id; the column names.is_primary indicates which is the definitive name for the organism (this isn't present in the table downloaded from NCBI but is derived from other attributes).

The taxonomic hierarchy is defined in nodes using an adjacency list model. It turns out that relational databases supporting recursive queries can easily traverse the graph defined by the model for a given node. This is well-described elsewhere (the linked post compares the adjacency model with other representations of trees in sql).

The query that we will use is called a recursive CTE (common table expression):

WITH RECURSIVE a AS (
 SELECT tax_id, parent_id, rank
  FROM nodes
  WHERE tax_id = '1280'
UNION ALL
 SELECT p.tax_id, p.parent_id, p.rank
  FROM a JOIN nodes p ON a.parent_id = p.tax_id
)
SELECT * FROM a;

And voila: we have traversed the tree from leaf to root, and in doing so, the lineage is defined. I'm using postgresql for this example, but sqlite supports recursive CTEs as well, and an identical query can be used with both databases.

At this point, we can easily add taxonomic names as well:

WITH RECURSIVE a AS (
 SELECT tax_id, parent_id, rank
  FROM nodes
  WHERE tax_id = '1280'
UNION ALL
 SELECT p.tax_id, p.parent_id, p.rank
  FROM a JOIN nodes p ON a.parent_id = p.tax_id
)
SELECT a.tax_id, a.rank, names.tax_name
  FROM a JOIN names USING(tax_id)
  WHERE is_primary;

It's worth pausing to point out an important feature of the representation of the adjacency list model of the nodes table downloaded from NCBI: the root is identified as the node for which tax_id = 1, but the value of parent_id for the root node is also 1! The recursive CTE as written expects a value of NULL in the parent_id field of the root node; with a value of 1, it will never terminate, and your laptop will get quite hot! So do yourself a favor and execute the following before trying this as home:

UPDATE nodes set parent_id = NULL where tax_id = 1;

As I mentioned, all of this is pretty well documented (as I learned once I bothered to look). But what about a query for multiple tax_ids at once? Is it as easy as including more than one in the query?

WITH RECURSIVE a AS (
 SELECT tax_id, parent_id, rank
  FROM nodes
  WHERE tax_id in ('562', '1280')
UNION ALL
 SELECT p.tax_id, p.parent_id, p.rank
  FROM a JOIN nodes p ON a.parent_id = p.tax_id
)
SELECT a.tax_id, a.rank, names.tax_name
  FROM a JOIN names USING(tax_id)
  WHERE is_primary;

Hmm, no, that doesn't work so well. Even if the ordering was right, we would find it difficult to keep track of the individual lineages. So maybe we can keep track?

WITH RECURSIVE a AS (
 SELECT tax_id as tid, 1 as ord, tax_id, parent_id, rank
  FROM nodes
  WHERE tax_id in ('562', '1280')
UNION ALL
 SELECT a.tid, a.ord + 1, p.tax_id, p.parent_id, p.rank
  FROM a JOIN nodes p ON a.parent_id = p.tax_id
)
SELECT a.tid, a.ord, a.tax_id, a.rank, names.tax_name
  FROM a JOIN names USING(tax_id)
  WHERE is_primary
  ORDER BY tid, ord;

Ok, that's better! If we like, we can group and aggregate the lineages:

(The results are shown using the equivalent of psql -x; the first line just provides a header row.)

SELECT 'value' as key;
WITH RECURSIVE a AS (
 SELECT tax_id as tid, 1 as ord, tax_id, parent_id, rank
  FROM nodes
  WHERE tax_id in ('562', '1280')
UNION ALL
 SELECT a.tid, a.ord + 1, p.tax_id, p.parent_id, p.rank
  FROM a JOIN nodes p ON a.parent_id = p.tax_id
)
SELECT a.tid, array_agg(a.tax_id), array_agg(a.rank), array_agg(names.tax_name)
  FROM a JOIN names USING(tax_id)
  WHERE is_primary
  GROUP BY tid;

Here's another approach: accumulate values in an array in the recursive expression, and then use the first value in the array to identify the first tax_id in each lineage.

WITH RECURSIVE a AS (
 SELECT tax_id as tid,
	parent_id as pid,
	rank,
	tax_id,
	ARRAY[tax_id] as lineage
  FROM nodes
  WHERE tax_id = '1280'
UNION ALL
  SELECT p.tax_id,
	 p.parent_id,
	 p.rank,
	 lineage[1],
	 lineage || ARRAY[p.tax_id]
  FROM a JOIN nodes p ON a.pid = p.tax_id
)
SELECT * FROM a;

We only actually care about the final lineage when it terminates at the root node.

WITH RECURSIVE a AS (
 SELECT tax_id as tid,
	parent_id as pid,
	rank,
	tax_id,
	ARRAY[tax_id] as lineage
  FROM nodes
  WHERE tax_id = '1280'
UNION ALL
  SELECT p.tax_id,
	 p.parent_id,
	 p.rank,
	 lineage[1],
	 lineage || ARRAY[p.tax_id]
  FROM a JOIN nodes p ON a.pid = p.tax_id
)
SELECT tax_id, lineage FROM a
WHERE a.rank = 'root';

With some embellishment, we can see how this strategy can be used to retrieve lineages plus additional annotation for multiple tax_ids (at arbitrary ranks) as once.

SELECT 'value' as key;
WITH RECURSIVE a AS (
 SELECT tax_id as tid,
	parent_id as pid,
	rank,
	tax_id,
	ARRAY[tax_id] as lineage,
	ARRAY[rank] as ranks
  FROM nodes
  WHERE tax_id in ('1279', '1280', '562')
UNION ALL
  SELECT p.tax_id,
	 p.parent_id,
	 p.rank,
	 lineage[1],
	 lineage || ARRAY[p.tax_id],
	 ranks || ARRAY[p.rank]
  FROM a JOIN nodes p ON a.pid = p.tax_id
)
SELECT tax_id, tax_name, nodes.rank, lineage, ranks
FROM a
JOIN nodes USING(tax_id)
JOIN names USING(tax_id)
WHERE a.rank = 'root'
AND names.is_primary;

How well does this scale? Both of the approaches executed in less than one second for 1000 tax_id's, but took more like 40s for 2000, so not fantastically well. But I expect the some optimization is possible.

Unfortunately, in the absence of an ARRAY data type, this approach does not work for sqlite, but you could do something similar by concatenating strings (though I doubt that this would be very efficient).