Identifying Microbes in Space | Towards Data Science

living in the International Space Station for over a quarter of a century. As the old tomato at the back of your fridge will tell you, that’s plenty of time to start growing mould. What microbes are tough enough to survive on the International Space Station with stellar radiation and no gravity? We can find the answer using metagenomics.

Let’s walk through how to find out what’s living on the International Space Station. You take your unlabelled DNA sequence and try to find a very similar sequence from a known species, which is very similar to searching for similar documents or web pages. As such biologists have built on the ideas of locality sensitive hashing to identify species. Here I’ll go through the code, discuss the results (including a surprisingly tasty find), and describe the algorithm used to identify species. Once you understand these techniques you can go and find out what mould is living anywhere in the solar system!

In metagenomics, rather than studying the DNA of just one organism, you look at the DNA of all the different species in one place. This paper by Urbaniak et al. used metagenomics to study all the bacteria and fungi from 8 different locations around the International Space Station. We’ll focus on just one of their samples; they took a piece of cloth (sterile, scientific cloth) and wiped it on the dining table, extracted all the DNA from that cloth, and sequenced it. Because this was a metagenomic study, the sampling method captures DNA from anything that might be living on the table. Crucially, in metagenomics you do not need to know what’s there before you can look at it’s DNA.

Analysing the DNA

Urbaniak et al. made their sequencing data available from the Sequence Read Archive under a creative commons zero licence. I downloaded sample F4_S5_P from the European Nucleotide Archive because it’s easier to work with. If you’d like to choose a different sample you can see all the details on where they’re from around the station on NASA’s website.

# Download PMA treated dining table sample
wget -P data/rawreads -nc ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR171/065/SRR17140465/SRR17140465_2.fastq.gz
wget -P data/rawreads -nc ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR171/065/SRR17140465/SRR17140465_1.fastq.gz

In raw reads, some sequences may be full of errors, have adapter sequences, or be too short to be useful. So before we analyse them we can trim and filter the reads to get rid of any low quality reads or sections of reads. I’m using fastp because it’s quick and automatically detects adapters.

fastp \
    -i data/rawreads/${accession}_1.fastq.gz \
    -I data/rawreads/${accession}_2.fastq.gz \
    -o data/fastp/${accession}_1.fastq.gz \
    -O data/fastp/${accession}_2.fastq.gz \
    -h data/fastp/${accession}.html \
    -j data/fastp/${accession}.json \

Now that we have high quality DNA sequences from a dinner table in space, it’s time to find out what was living up there. kraken2 is a very fast and widely used metagenomics tool for identifying what species a DNA sequence is from. kraken2 is much faster than BLAST, another tool you may have heard of. BLAST is very good if you have 10s or even 100s of sequences to identify, but just for our one sample we have >75,000 pairs of reads! And this is a small read set, it’s not uncommon for modern sequencing runs to have millions of reads. With all these reads we need something very fast, like kraken2.

Why is kraken2 so much faster than BLAST? BLAST is doing alignment, lining up your sequence against possible matches and accounting for differences between them. That way, even if your sequence doesn’t exactly match another similar sequence, BLAST still knows exactly how they’re similar and how they’re different. This can be very useful information, but it is quite slow.

# BLAST alignment
CAAATAATTTGAAA
||  || |||||||
CACGTATTTTGAAA

In kraken2 we get a big speed boost by only looking for exact matches. We don’t try and match the whole sequence all at once, instead we break it up into shorter pieces of length $k$ called $k$-mers. $k$-mers don’t count as matching unless they are exactly the same, but if sequences are similar some of the multiple $k$-mers do match exactly and we get a hit.

# Query sequence
CAAATAATTTGAAA
# k-mers
CAAATAATT
 AAATAATTT
  AATAATTTG
   ATAATTTGA
    TAATTTGAA
     AATTTGAAA

# Similar sequence with 4/6 identical k-mers
CGAATAATTTGAAA
# k-mers
CGAATAATT           No exact match
 GAATAATTT          No exact match
  AATAATTTG             Exact match
   ATAATTTGA            Exact match
    TAATTTGAA           Exact match
     AATTTGAAA          Exact match

An important point to keep in mind is that kraken2 is comparing your sequences to a database. This means you can only identify the correct species, if it is in your database. If the species your DNA comes from is not in the database, you may get no answer or a false match against another closely related species. Ben Langmead, one of the developers of kraken2, makes several databases available on this website. I selected a database with fungi as well as bacteria, and to make things faster I chose one capped at 8GB, the PlusPF-8 database. Once I had downloaded and extracted the database I ran kraken2 for paired reads, with the --memory-mapping setting to make it easier on my laptop.

w get -P ~/databases https://genome-idx.s3.amazonaws.com/kraken/k2_pluspf_08_GB_20260226.tar.gz

tar -xvzf ~/databases/k2_pluspf_08_GB_20260226.tar.gz -C ~/databases

kraken2 \
    --db ~/databases/k2_pluspf_08_GB_20260226 \
    --paired \
    data/fastp/${accession}_1.fastq.gz \
    data/fastp/${accession}_2.fastq.gz \
    --threads 12 \
    --memory-mapping \
    --report data/kraken2/${accession}_kraken2_report.txt \
    --output data/kraken2/${accession}_kraken2_output.txt

What’s on the table?

The kraken2 report tells you how many reads mapped to each observed taxon, see kraken2 wiki for full details. The report tells you the number of reads mapping exactly to this taxon (column 2) and also the number of reads mapping to any of its descendant taxons (column 3). For example, the genus Acinetobacter is assigned 134 reads directly, but 238 reads if you include the reads assigned to the nine observed Acinetobacter species and two sub-genera. (I’ve added the column titles to help).

%       descendant# exact#  rank    txid    name
0.32	238     	134 	G   	469	     Acinetobacter
0.10	76      	76  	S   	40214	       Acinetobacter johnsonii
0.02	12      	0   	G1  	196816	       unclassified Acinetobacter
0.01	9       	9   	S   	2743575	         Acinetobacter sp. NEB 394
0.00	3       	3   	S   	2953738	         Acinetobacter sp. Z1
0.01	7       	0   	G1  	909768	       Acinetobacter calcoaceticus/baumannii complex
0.01	7       	7   	S   	470	         Acinetobacter baumannii
0.00	3       	3   	S   	1776742	       Acinetobacter vivianii
0.00	3       	3   	S   	40215	       Acinetobacter junii
0.00	1       	1   	S   	29430	       Acinetobacter haemolyticus
0.00	1       	1   	S   	106649	       Acinetobacter guillouiae
0.00	1       	1   	S   	28090	       Acinetobacter lwoffii

We can sort the table by number of reads and filter taxon rank to S to show the most common species found on dining table in the International Space Station. We can also plot this data using figtreemap, see the source code for exactly how. Lineage information was added using taxonkit to make it easier to colour the plot.

cat data/kraken2/SRR17140465_kraken2_report.txt | \
taxonkit lineage \
    --data-dir ~/databases/taxdump/ \
    -i 5 \
    > data/kraken2/SRR17140465_kraken2_report_lineage.txt

figtreemap.squarify_images.figtreemap(
    sizes,
    imgs,
    facecolor ="black",
    edgecolor ="white",
    label = [name.replace(" ", "\n") if size > .1 else "" for name, size in zip(names,sizes)],
    text_kwargs={
        "color":"white",
        "size": "xx-small",
    }
)
plt.show()

The figure above illustrates the species detected and how abundant their DNA was in the sample. The size of each species illustration is proportional to the number of reads from that species in the sample. Blue is used for bacteria and orange for fungi. The most common species detected are Pseudolactococcus raffinolactis (found in raw milk and used in food production), Cutibacterium acnes (found on human skin, related to acne), Ralstonia pickettii (which is an environmental bacteria which has caused some hospital outbreaks by contaminating medical equipment), and most interestingly Leuconostoc mesenteroides.

L. mesenteroides is used in fermentation, so I wanted to try and find out they had been fermenting. I used seqkit to extract the reads identified for this species and megahit to assemble them. The assembly wasn’t great, but it was good enough to BLAST the contigs and get more details about the strain of L. mesenteroides.

# Extract reads for Leuconostoc mesenteroides txid 1245
seqkit grep \
    -r \
    -n \
    -p "kraken:taxid\|1245" \
    --threads 4 \
    data/kraken2/${accession}_kraken2_classified_1.fastq > data/kraken2/${accession}_1245_1.fastq

# Assemble extracted reads
megahit \
    -1 data/kraken2/${accession}_1245_1.fastq \
    -2 data/kraken2/${accession}_1245_2.fastq \
    -o data/megahit/1245

Several contigs, including the longest, map perfectly to strain MSL129 which was isolated from kimchi! And we know that kimchi was eaten on the ISS because in 2008 Korea’s first astronaut was sent up with kimchi specially developed for space by the Korea Atomic Energy Research Institute! Before sending the kimchi to space all bacteria was killed by high dose radiation, described in this paper. However, the DNA we’ve detected here is from intact, presumably viable cells because the sample was treated with PMA. So somehow, this kimchi bacteria has made it to space and now that it’s there it is one of the most abundant species on the ISS dinner table!

This is a technical description of how kraken2 rapidly assigns taxonomy to sequences. Kraken2 takes your query sequence and splits it into $k$-mers, keeps only a subset of those $k$-mers, converts those remaining $k$-mers to a number called a hash, and then looks the hash up in a database to see which species it is from. There are some very clever tricks being used to make this very fast, both in how the $k$-mers are used and how the database is built.

Overlapping $k$-mers have a lot of redundancy, adjacent $k$-mers have mostly the same sequence. “Minimizers” are a technique to avoid dealing with multiple almost identical $k$-mers. By taking the minimizer of a set of similar $k$-mers you can get a single subsequence to represent multiple the $k$-mers. To get the minimizer, each $k$-mer is broken into smaller $k$-mers (we use $l$ for this smaller length and call them $l$-mers to avoid confusion). You only keep the $l$-mer that is alphabetically first, that is the minimizer of the $k$-mer. Similar $k$-mers will give you the same minimizer so you don’t have to save each $k$-mer.

Let illustrate this with an example using $k=16$ and $l=10$. Our sequence is 25 bases long and produces 10 different $k$-mers. But those $k$-mers collapse down to only 3 different minimizers.

# Sequence
CCTGAGTGTGGAGACCCAAGTGAGA

# kmers, k = 16
CCTGAGTGTGGAGACC
 CTGAGTGTGGAGACCC
  TGAGTGTGGAGACCCA
   GAGTGTGGAGACCCAA
    AGTGTGGAGACCCAAG
     GTGTGGAGACCCAAGT
      TGTGGAGACCCAAGTG
       GTGGAGACCCAAGTGA
        TGGAGACCCAAGTGAG
         GGAGACCCAAGTGAGA

k-mer	            Minimizer
CCTGAGTGTGGAGACC	AGTGTGGAGA    1
CTGAGTGTGGAGACCC	AGTGTGGAGA
TGAGTGTGGAGACCCA	AGTGTGGAGA
GAGTGTGGAGACCCAA	AGTGTGGAGA
AGTGTGGAGACCCAAG	AGTGTGGAGA
GTGTGGAGACCCAAGT	AGACCCAAGT    2
TGTGGAGACCCAAGTG	AGACCCAAGT
GTGGAGACCCAAGTGA	ACCCAAGTGA    3
TGGAGACCCAAGTGAG	ACCCAAGTGA
GGAGACCCAAGTGAGA	ACCCAAGTGA

kraken2 actually ignores some positions in the minimizer sequence; counter intuitively this improves read assignment. $s$ alternating positions at the end of the minimizer are ignored, so if $s=3$ AGTGTGGAGA is effectively AGTG-G-A-A. Ignoring these bases means that differences at those 3 positions don’t break the match. Therefore, kraken2 can still identify species despite mutations and sequencing errors.

The next step is to “hash” these minimizers. A hash is effectively a random number but the same sequence will always give you the same number. When building the database this hash is then stored in a row with the species it is from. Kraken2 uses a clever approach called a compact hash code to make the database smaller and faster to search. Kraken2 splits the hash into two parts (e.g. 123456 becomes 123 and 456) and only saves the first part in the database. Which sounds like madness, but the second part is encoded as the position in the database. So in this example 123 would be saved in row 456 of the database. It’s not exactly 456 though, it checks for an empty row to put the record in starting at 456, so strictly 123 would be saved in row 456 or greater. This compact hash code makes the database smaller because we’re not saving all of the hash, and faster to search because we know where to start looking.

Once you have the taxonomic assignments for each $k$-mer in your read, how do you combine these to assign a taxa to the whole read? After all, each read has multiple $k$-mers which could disagree with each other. The read is assigned to whichever taxa has the most $k$-mer hits itself or its ancestors. And that’s how kraken2 can assign taxonomy to a sequence with interstellar speed 🚀.

The core problem that we have solved here is taking a DNA sequence of unknown origin and rapidly comparing it to a large number of other labelled sequences that are imperfect matches. This is almost exactly the same problem as looking for near identical webpages or detecting plagiarism. So it’s no surprise to see techniques like minhashing and locality sensitive hashing popping up in these bioinformatics tools like kraken2 and sourmash.

Traditional tools for biological sequence comparison like BLAST were more rooted in evolutionary theory but as the amount of sequence data available keeps increasing bioinformaticians have had to turn to faster methods which place less emphasis on the biological origins of the sequence.

By reducing your problem to its core question, you can see where it overlaps with problems in different fields and explore their solutions.