---
tags: stamps2025
---

# Metagenomics intro IV: comprehensive metagenome content analysis and the challenge of bacterial strains (STAMPS 2025)

Titus Brown
July 17, 2025

## Identifying _all_ overlapping genomes for a metagenome

conda activate sourmash

mkdir -p ~/day3
cd ~/day3

sourmash prefetch -k 51 --scaled 10_000 \
    /opt/shared-2/sourmash-mgx/SRR11125891.sig \
    /opt/shared-2/sourmash-db/gtdb-rs226.k51.rocksdb \
    -o SRR11125891.prefetch.csv

This calculates _all_ overlaps between the metagenome and the database, and will... take a few minutes.

How many matches are there?

Why, and what's going on? :sob:

## The challenge of bacterial strains!

Let's take a step back, and look at one particular species - the one we chose to map to, which in NCBI land is `Methanobrevibacter smithii`, but in GTDB-land is known as `s__Methanocatella`.

Let's pull out all the genome names that GTDB says belong to `s__Methanocatella`:

sourmash tax grep s__Methanocatella \
    -t /opt/shared-2/sourmash-db/gtdb-rs226.lineages.csv \
    -o methano.tax.csv

How many do we get? 349!

Now let's grab these 349 genome sketches:

sourmash sig cat --picklist methano.tax.csv:ident:identprefix \
    /opt/shared-2/sourmash-db/gtdb-rs226.k31.sig.zip \
    -o methano.k31.sig.zip

and compare them:

sourmash compare methano.k31.sig.zip -o methano.cmp
sourmash plot methano.cmp

Check out `methano.cmp.matrix.png` - what do we see?

Let's ask how many of these match to our metagenome:

sourmash scripts manysearch \
    methano.k31.sig.zip \
    /opt/shared-2/sourmash-mgx/SRR11125891.sig \
    -o SRR11125891.methano.manysearch.csv

and plot that presence diagram...

sourmash scripts presence_filter \
    SRR11125891.methano.manysearch.csv \
    -o SRR11125891.methano.presence.png \
    -N 0 \
    --detection-column-name containment 

You will see this - what does it mean??


## Identifying a _usefully small_ set of overlapping genomes for a metagenome

What 'gather' does is choose a minimal set of genomes necessary to represent the known content. It does that by (essentially) taking the genome with the largest fraction detected, and then removing anything that overlaps with it: a "greedy" approach of taking the best match, first.

Let's run it again:

conda activate sourmash

sourmash gather -k 51 --scaled 10_000 \
    /opt/shared-2/sourmash-mgx/SRR11125891.sig \
    /opt/shared-2/sourmash-db/entire-2025-07-11.k51.rocksdb \
    --threshold-bp=0 \
    -o SRR11125891.t0.gather.csv

Next, to illustrate the method of taxonomy construction, we're going to run `sourmash tax annotate`, which annotates the CSV file with taxonomic lineages:

sourmash tax annotate -g SRR11125891.t0.gather.csv \
     -t /opt/shared-2/sourmash-db/entire-2025-07-11.lineages.sqldb

Check out our new file `SRR11125891.t0.gather.with-lineages.csv` - how does this differ from `SRR11125891.t0.gather.csv`? Hint - look at the last column...

And finally we're going to run `taxburst` on it:

conda activate taxburst
taxburst -F tax_annotate \
    SRR11125891.t0.gather.with-lineages.csv \
    -o SRR11125891.t0.taxburst.html

:::danger
If `conda activate taxburst` doesn't work, you may need to run:

conda create -n taxburst -y python=3.12

followed by

conda activate taxburst
pip install -U taxburst

to install taxburst.

Then run the taxburst command above.
:::

### Looking at the taxburst view

Open the taxburst file.
What we did was bring in all the genome level information -
this is a pretty busy diagram, what with ~1500 genomes on it!

Try decreasing the max depth here ;).

## What are other ways to deal with massive strain overlap?

A common way to deal with genomes that have lots of shared sequence is to dereplicate them to "species representatives". This is an approach that GTDB uses to produce it's "reps" database, and it's also something that MetaPhlan and others do. I'm not sure I understand the consequences of doing things this way, however.

But it's still an evolving field.

## Annotated bibliography

Mash: fast genome and metagenome distance estimation using MinHash, Ondov et al. [link](https://genomebiology.biomedcentral.com/articles/10.1186/s13059-016-0997-x) The OG "let's use clever techniques to compare genomes using a subset of k-mers" paper that inspired sourmash :).

Lightweight compositional analysis of metagenomes with FracMinHash and minimum metagenome covers, Irber et al., 2023. [link](https://www.biorxiv.org/content/10.1101/2022.01.11.475838v2). This is the core sourmash paper that explains how sourmash works and what 'sourmash gather' does.