---
tags: stamps2025
---

[toc]

# Metagenomics intro II: bioinformatics techniques for measuring genomes in metagenomes (STAMPS 2025)

Titus Brown
July 17, 2025

## What do genomes look like?

Run:

head /opt/shared-2/genomes/GCF_037236875.1_genomic.fna


This is [Methanobrevibacter smithii](https://www.ncbi.nlm.nih.gov/search/all/?term=GCF_037236875.1), an archaeal species that is (probably) present in our pig gut metagenome, SRR11125891.

It's in [FASTA format](https://en.wikipedia.org/wiki/FASTA_format).

Count number of characters:

wc /opt/shared-2/genomes/GCF_037236875.1_genomic.fna

There are 1,826,901 characters in the file, or roughly 1.8 million base pairs.

Examine size:

ls -lh /opt/shared-2/genomes/GCF_037236875.1_genomic.fna
-rw-rw-r-- 1 stamps stamps 1.8M Jul 17 06:09 /opt/shared-2/genomes/GCF_037236875.1_genomic.fna

It's about 1.8 MB in size. This is pretty typical (albeit a bit small) for a bacterial or archaeal genome.

This archaeal genome has been assembled in chunks, or _contigs_ - how many contigs?

grep ^'>' /opt/shared-2/genomes/GCF_037236875.1_genomic.fna | wc -l

Approximately 30. Why??

## What do metagenomes look like?

Run:

gunzip -c /opt/shared-2/metagenomes/SRR11125891_1.fastq.gz | head


Basically: many, many reads in [FASTQ format](https://en.wikipedia.org/wiki/FASTQ_format) :).

To count lines:

gunzip -c /opt/shared-2/metagenomes/SRR11125891_1.fastq.gz | wc -l

This will tell you that there are 14732348 lines in that file. How many records are there? How many bases are there? How many bases are there in this metagenome overall, roughly?

What does this show?

ls -lh /opt/shared-2/metagenomes/SRR11125891_*

tl;dr genomes are ~small, metagenomes are BIG. (This is actually a pretty small one.)

## Reference based approaches: mapping and k-mers

"Analyzing" the presence of a genome is a common (perhaps the most common?) way to embark on metagenome analysis. Once you know some or all of a particular genome is present, you can then "transfer" information from the genome and genome annotation onto the metagenome.

Other techniques include looking for some or all known genes, analyzing marker genes for taxonomy, assembling into contigs and analyzing those, and, err... other things?

## Introducing mapping

The reads in a metagenome are drawn ~randomly from a weighted mixture (see: dart board analogy). There are several approaches to analyzing the reads and linking them to known information. One of the most common approaches is **mapping**.

In **mapping**, you try to "place" each read on a (at least one) reference genome. Mapping is the basis of many detailed analysis tools, including parts of anvi'o (Saturday!), inStrain, and the BioBakery software (MetaPhlan, StrainPhlan, Humann).

### Let's run some mapping!

Prepare work folder and install some software:

mkdir -p ~/day3-mapping
cd ~/day3-mapping

conda create -n mapping -y minimap2 samtools

Stickies up when conda finishes. (& then wait.)

Stickies down.

Then:

conda activate mapping

minimap2 -ax sr /opt/shared-2/genomes/GCF_037236875.1_genomic.fna \
    /opt/shared-2/metagenomes/SRR11125891_1.fastq.gz \
    /opt/shared-2/metagenomes/SRR11125891_2.fastq.gz \
    | samtools view -b > GCF_037236875.mapping.bam

This will take about 30 seconds - mappers are *fast*.

Do some reorganizing and indexing of data:

samtools sort GCF_037236875.mapping.bam -o GCF_037236875.mapping.sorted.bam
samtools index GCF_037236875.mapping.sorted.bam

And now let's look at the actual placement of reads!

samtools tview GCF_037236875.mapping.sorted.bam /opt/shared-2/genomes/GCF_037236875.1_genomic.fna

This opens up a terminal viewer. Commands:
* 'q' to quit
* '.' to toggle between consensus view and read mapped view
* left/right to move left and right in the contig
* up/down to look at more reads
* 'r' to toggle read name
* 'g' will let you select a new position; ESC to get out of that window.

What things do you notice?

### Summary statistics

We can also get some stats on things like coverage of each contig:

samtools coverage GCF_037236875.mapping.sorted.bam

The interesting columns are:
* numreads
* coverage - what percentage of bases are "hit" by at least one read?
* meandepth - on average, how many reads cover the bases?

You can also get read mapping percentage:

samtools flagstat GCF_037236875.mapping.sorted.bam

note that only 224689 reads mapped, out of 7366697 total!!

So, this is mapping. It underlies many analysis types. Iva will talk a bit about functional annotation of genomes on Saturday; please ask her about how it relates to mapping, too!

### Discussion points:

Things to discuss (titus will edit some answers in hackmd as we go ;):

#### Q: how many known genomes are there, roughly?

Somewhere north of 2m. There's actually about ~2m isolate genomes that aren't yet in the databases; see AllTheBacteria project.

#### Q: how do you identify what genome to map to?

This is a hard question! There are a variety of approaches!

You can map to many of them at once (centrifuge).

You can identify marker genes in your metagenome, and then identify which species are there, and map to those species genomes.

Another approach is to use k-mers to identify likely species. This technique is used by Kraken, Kaiju, sourmash, and sylph, among others. This is discussed next!

16S! You can probably use 16S to do this! But I'm not sure anyone does.

#### Q: what happens if the reference genome is incomplete or mismatched?

Great question! Mapping is _entirely_ reference based. If the reference genome doesn't match to the reads, then the reference genome will be ignored for that particular analysis.

#### Q: what happens if the metagenome is "missing" some or all of the reference genome for a species?

(Why might this happen?)

This happens a lot, and we'll talk about some specific reasons why later. Basically, it'll be the basis for the whole afternoon :sweat_smile:

#### Q: what's the best mapping software to use?

Start with minimap2 or bowtie2 for short reads.

Not sure for long reads.

It (mostly) is not critically important. The parameters for mapping, as well as filtering afterwards, are more important!

But, realistically, the reference genomes you use are MUCH more important than the mapping software you use (in the sense that many commonly used mappers work great already, but you're very much at the mercy of your reference genomes, whatever mapper you use!)

## Introducing k-mers

### k-mers are fixed-length (k) subsequences of DNA (or protein, or...)

Consider the 150 bp read:

CCCTCAGGCAGTCCGAAAGCTAATTTATACAACAAACGTTATAGAGGGATTTAACCGTCAACTCCGTAAAGTCACAAAGTCAAAATCTGTATTTCCAACAGATGACAGCCTGTTTAAAATGCTGTATCTGGCGATGATAGACATCACGAAA

or the first 80 characters or so of the genome above:

GTCACTTTCAAAAACTTGTGTGATTTGAATTTTCACTAGTTTTTGAGTCACGAATTTTTTTCGTTCTATTTTTTCTTTAT

Run a sliding window of length 21 across the read:

CCCTCAGGCAGTCCGAAAGCTAATTTATACAACAAACGTTATAGA...
CCCTCAGGCAGTCCGAAAGCT
 CCTCAGGCAGTCCGAAAGCTA
  CTCAGGCAGTCCGAAAGCTAA
   ...

and collect all subsequences.

Do the same for the contig sequence.

Now, ask: are there any shared k-mers (exact matches)? If so, that indicates that (maybe) these sequences are related.

That is really the central point of long k-mers in metagenomics: they let you infer potential overlap between different sequences, based on exact matches.

Titus to fill in below as we talk:

#### How do you choose a k-mer size?

Hint: there are 4**k distinct k-mers.

Think: k=21, k=31, k=51... Why these numbers?

#### Why is 'k' often odd?

#### How many k-mers are there in a (non-repetitive) genome of 1 million base pairs?

Does it depend on k?

#### What happens with k-mers when sequence has lots of errors?

#### Can we say anything about the likely sensitivity and specificity of k-mers?

If you see a match between two sequences, how likely is it to be a false positive (i.e. those two sequences are not real?)

#### Where do k-mers come from?

#### In what way are k-mers computationally efficient?

### Examining k-mer statistics with sourmash

sourmash is a bioinformatics "multitool" for playing with k-mers. Let's try it out!

#### Sketching some sequencse

mkdir -p ~/day3-kmers
cd ~/day3-kmers

conda activate sourmash

Now, convert the metagenome and the genome into k-mers:

sourmash sketch dna -p k=21,k=31,k=51 \
    /opt/shared-2/genomes/GCF_037236875.1_genomic.fna \
    -o GCF_037236875.1.sig.zip \
    --name GCF_037236875.1
    
sourmash scripts singlesketch -p k=21,k=31,k=51,abund \
    /opt/shared-2/metagenomes/SRR11125891_*.fastq.gz \
    -o SRR11125891.sig.zip \
    --name SRR11125891

This produces two files - sourmash "signatures" - that contain the k-mers in the genome sequence and read sequences, respectively. "Signatures" contain both the k-mers and some extra information like a name.

Discussion points:
* you have to pick one or more k-mer sizes (well, the default is k=31);
* by default, sourmash just tracks presence/absence; for metagenomes, however, you want to track the *number of times* you see a k-mer, as that can offer valuable information about genome coverage.
* by default, sourmash combines all the sequences in all the files you give it into one signature. You can change this, of course.
* you can use `sketch dna` or `scripts singlesketch` equally. The latter is just much faster for reads. The output is the same.

Link: [sourmash sketch documentation](https://sourmash.readthedocs.io/en/latest/sourmash-sketch.html)

#### Ask about the overlap!

Now let's ask a basic question: how much of this genome is in this metagenome??

sourmash scripts manysearch GCF_037236875.1.sig.zip \
    SRR11125891.sig.zip \
    -o manysearch.csv

Link: [branchwater plugin documentation for 'manysearch'](https://github.com/sourmash-bio/sourmash_plugin_branchwater/tree/main/doc#Running-manysearch)

You will see:

query             p_genome avg_abund   p_metag   metagenome name
--------          -------- ---------   -------   ---------------
GCF_037236875.1     92.6%    14.4       2.8%     SRR11125891

Here,
* p_genome is the fraction of the genome's k-mers present at least once in the metagenome (presence/absence);
* avg_abund is the average number of times each k-mer shows up in the metagenome;
* p_metag is the fraction of the _metagenome_ k-mers that are present in the genome. This loosely approximates how much of the metagenome data will map to the genome.

Discussion points:
* what k-mer size is being used? We have three to choose from!
* what happens if you change the k-mer size to 21 or 51 with `-k 21` or `-k 51`?
* how do these numbers compare with 'samtools coverage', above?

#### A weighted Venn diagram

Ugly, but maybe helps build intuition?

sourmash scripts weighted_venn SRR11125891.sig.zip GCF_037236875.1.sig.zip -o venn.png

#### Some intuition around k-mers, and also abundance

Remember the samtools tview output above? Let's look at it through the lens of k-mers...

Consider (same views, just toggled b/t showing sequence and not):



* What will you see if you count k-mers here?
* Where will k-mers not match the reference?

## Discussion: tradeoffs between approaches

We have:
* short reads
* long reads
* genomes

Meditate on the following questions:

* It is *easy* to match reads to genomes with mapping, and also with k-mers. Why?
* It is *hard* to match reads to reads with mapping. Why?
* Why doesn't everyone use k-mers and not mapping?
* Why doesn't everyone use mapping and not k-mers?

Think about them for 2 minutes. Then chat with a seat mate near you for 5 minutes. Then share with the class!

## My favorite k-mer metrics!!

You can do several things with k-mers and overlaps; here a suggestive Venn diagram showing intersections/overlap between genomes A, B, and C.


Here is some math for comparing genomes based on presence/absence of k-mers, and the sourmash commands that produce numbers based on that comparison:


(credit: Cassandra Olivas, UC Davis)

Try comparing our genome and our metagenome:

sourmash sig overlap GCF_037236875.1.sig.zip \
    SRR11125891.sig.zip -k 31 

you should get:

--- Similarity measures ---
jaccard similarity:          0.00678
first contained in second:   0.92605 (cANI: 0.99752)
second contained in first:   0.00678 (cANI: 0.85123)
average containment ANI:     0.92438


--- Hash overlap summary ---
number of hashes in first:   1839
number of hashes in second:  251092

number of hashes in common:  1703
only in first:               136
only in second:              249389
total (union):   

Can you estimate these numbers between sets of reads?

## sourmash and k-mers

sourmash offers a lot of different functionality that is almost all based around using a small subset of k-mers to estimate more complete statistics. This makes it quite fast in many circumstances. See the bibliography for some pointers here, and I'm happy to chat more about how and why it works.

Here's a venn diagram view that might help motivate the secret sauce:


## K-mer metrics and relationship to mapping and ANI

Sequence is sequence, and it turns out you can connect one type of measurement to another just fine - at least within a particular taxonomic/evolutionary range. For DNA k-mers in the 21-51 base size the taxonomic range is genus, species, and strain.

Some things we've seen or can infer:
* k-mers are sensitive, specific, and fast;
* they're not very good if you have lots of mismatches (a distant reference genome, or high strain variation) or errors (Nanopore sequencing); mapping is better here!
* there is a point at which sequences can still map but k-mers cannot be matched.
* However, mapping based approaches struggle to distinguish betwen strains.

We'll discuss all of this later, too! For now, just know that you can estimate a lot of mapping-based statistics such as coverage/depth, fraction detected, mapping fraction, and Average Nucleotide Identity quickly and easily with k-mers.

### Generally speaking, k-mer matches <-> mapping

Consider for this gut metagenome: red dots are 31-mer matches, blue dots are mapping matches.


## Annotated bibliography

A detailed sourmash tutorial about what _exactly_ it is doing under the hood: [link](https://sourmash.readthedocs.io/en/latest/kmers-and-minhash.html)

[Large-scale sequence comparisons with sourmash](https://f1000research.com/articles/8-1006). Pierce-Ward et al., 2019. The first sourmash paper.

[Mash: fast genome and metagenome distance estimation using MinHash](https://genomebiology.biomedcentral.com/articles/10.1186/s13059-016-0997-x), Ondov et al., 2016. The OG "do-many-genome-comparisons-with-kmers".

MetaPalette: a k-mer Painting Approach for Metagenomic Taxonomic Profiling and Quantification of Novel Strain Variation, Koslicki and Falush, 2016. [link](https://journals.asm.org/doi/10.1128/msystems.00020-16) Why k-mers work well for bacteria in the first place, and which k-mer sizes to use; one of the two inspirations for sourmash.

[Deriving confidence intervals for mutation rates across a wide range of evolutionary distances using FracMinHash](https://pubmed.ncbi.nlm.nih.gov/37344105/). Rahman Hera et al., 2023. The math behind k-mers and mapping/ANI.

[Informed and automated k-mer size selection for genome assembly](https://academic.oup.com/bioinformatics/article/30/1/31/235479), Chikhi and Medvedev, 2013.

[Merqury: reference-free quality, completeness, and phasing assessment for genome assemblies](https://genomebiology.biomedcentral.com/articles/10.1186/s13059-020-02134-9), Rhie et al., 2020. Thoughts on k-mer size!