Identifying Airborne Viruses | About Data Science

living on the International Space Station for more than a quarter of a century. As the old tomato in the back of your fridge will tell you, that's plenty of time for mold to start growing. What bacteria are hardy enough to survive on the International Space Station with star radiation and zero gravity? We can find the answer using metagenomics.

Let's explore how we can find out what lives on the International Space Station. You take your unlabeled DNA sequence and try to find the same sequence in a known species, which is exactly the same as searching for similar documents or web pages. As such biologists have built on spatially sensitive hashing ideas to identify species. Here I'll walk through the code, discuss the results (including some surprisingly sweet finds), and explain the algorithm used to identify species. Once you understand these techniques you can go find out what fungi live anywhere on the solar planet!

In metagenomics, rather than studying the DNA of just one organism, you look at the DNA of all organisms in one place. This paper by Urbaniak et al. used metagenomics to study all bacteria and fungi from 8 different locations around the International Space Station. We will focus on just one of their samples; they took a piece of cloth (a sterile, scientific cloth) and wiped it on the dining room 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. Most importantly, in metagenomics you don't need to know what's there before you look at DNA.

DNA analysis

Urbanak et al. have made their sequence data available in the Read Sequence Archive under a zero creation license. I downloaded the F4_S5_P sample from the European Nucleotide Archive because it is easy to work with. If you would like to choose a different sample you can see all the details of where they are from the station on the NASA 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, contain adapter sequences, or be too short to be used. So before we analyze them we can trim and filter the reads to remove any low quality reads or parts of the reads. I use it fastp because it is fast 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 the dinner table in space, it's time to figure out what was living up there. kraken2 is a fast and widely used metagenomics tool for identifying which DNA sequence comes from which species. kraken2 is something a lot faster than BLAST, another tool you may have heard of. BLAST is great if you have 10s or 100s of sequences to index, but for our single sample we have >75,000 pairs of reads! And this is a small set of reads, it's not uncommon for a sequence to run into millions of reads. With all these studies we need something quick, like kraken2.

why is that kraken2 much faster than BLAST? BLAST performs an alignment, comparing your sequences against possible matches and counting the differences between them. That way, even if your sequence doesn't exactly match other similar sequences, BLAST still knows exactly how it's similar and how it's different. This can be very useful information, but it is slow.

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

In kraken2 we get a big increasing speed by only looking for exact matches. We don't try to put the whole program together at once, instead we break it up into shorter pieces of length $k$ called $k$-ama. $k$-mers are not counted as identical unless they are exactly the same, but if the sequences are identical in some of the repeats $k$-mers are exactly the same 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

The important point to remember is that kraken2 compares your sequence with the database. This means that you can only identify the correct species, if in your database. If your DNA type is not in the database, you may not get an answer or a false match to another closely related type. Ben Langmead, one of the developers of kraken2make several information sites available on this website. I chose a database with fungi and viruses, and to do things quickly I chose one installed on 8GB, the PlusPF-8 database. Once I downloaded and extracted the database I used it kraken2 paired reading, ne --memory-mapping editing makes it easy on my laptop.

w get -P ~/databases 

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?

I kraken2 the report tells you how many reads were mapped to each taxon monitored, see the kraken2 wiki for full details. The report tells you the number of mapping readings that directly correspond to this taxon (column 2) and the number of mapping readings for any of its parent taxa (column 3). For example, the genus Acinetobacteria assigned 134 reads directly, but 238 reads when you include the assigned reads of the nine noted. Acinetobacteria species and two subspecies. (I added column headings 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 readings and sort tax rank to S to show the most common species found on the International Space Station dining table. And we can edit this data using figtreemapsee the source code for how. Genealogy information was added using taxonkit to make it easier to color the structure.

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 shows the type of species found and how abundant their DNA was in the sample. The size of the representation for each species is proportional to the number of reads of that species in the sample. Blue is used for bacteria and orange for fungi. The most common types found are: Pseudolactococcus raffinolactis (found in raw milk and used in food production), Cutibacterium acnes (found in human skin, associated with acne), Ralstonia pickettii (which is the natural bacteria that caused the hospital outbreak by contaminating medical equipment), and very interesting Leuconostoc mesenteroides.

L. mesenteroides it is used in fermentation, so I wanted to try to find out how long it has been fermenting. I used seqkit to extract identified readings of this type and megahit to combine them. The assembly wasn't perfect, but it was good enough to ANALYZE the contigs and get more information about the species 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 a very long, complete map of the MSL129 repressor isolated from kimchi! And we know that kimchi was eaten on the ISS because in 2008 the first Korean astronaut was sent with kimchi specially developed in space by the Korea Atomic Energy Research Institute! Before sending kimchi into space all bacteria are killed by high-dose radiation, described in this paper. However, the DNA we obtained here is from intact cells, which may be active because the sample was treated with PMA. So somehow, this kimchi bacteria has made it to space and since it's there it's one of the most common varieties on the ISS dinner table!

This is a technical explanation of how kraken2 gives a quick taxonomy sequence. Kraken2 it takes your query sequence and splits it in half $k$-mers, it only stores a subset of those $k$-mers, converts the remainder $k$-mers to a number called a hash, then look up the hash in the database to see what type it comes from. There are very clever tricks used to make this very fast, both ways $k$-mers are used and how the database is structured.

Overlapping $k$-mers have lots of junk, around $k$-mers have multiple times in the same sequence. “Reducers” are a way to avoid interactions and duplicates that are almost identical $k$-ama. By taking the reducer of the same set $k$-mers you can find the following single to represent multiples $k$-ama. For minimizer, each $k$-mer is divided into small $k$-mers (we use $l$ in this small length and call them $l$– to avoid confusion). You only last $l$-mer i.e. first in alphabets, i.e. diminutive of $k$-um. The same $k$-mers will give you the same reducer so you don't need to save each one $k$-um.

Let's illustrate this with an example $k=16$ again $l=10$. Our sequence is 25 bases long and generates 10 variants $k$-ama. But those $k$-mers break down into only 3 different dimers.

# 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 it actually ignores other positions in the minimizer sequence; counter intuitively this improves the learning curve. $s$ alternating positions at the end of the reducer are ignored, so if $s=3$ AGTGTGGAGA is effectively AGTG-GAA. Ignoring these bases means that the difference in those 3 positions does not break the match. Therefore, kraken2 it can still recognize species despite genetic mutations and sequencing errors.

The next step is to “hash” these reducers. A hash is a random number but the same sequence will always give you the same number. When you create a database this hash is stored in the order of the types from which it appears. Kraken2 it uses a clever technique called collective hash code to make the database smaller and faster to search. Kraken2 it splits the hash into two parts (eg 123456 becomes 123 and 456) and saves only the first part in the database. Which sounds crazy, but the second part is coded as a backup. So in this example 123 will be stored in row 456 of the database. It's not exactly 456, but it looks for an empty line to include a record from 456, so 123 will be stored on line 456 or so. This compact hash code makes the database smaller because we don't store all the hashes, and it's faster to search because we know where to start looking.

Once you have each taxonomic assignment $k$-In your study, how do you combine this to give a taxa to all the studies? After all, each study has more $k$-people who may disagree. The readings are assigned to whichever tax has the most $k$-mer itself or its ancestors. And so it is kraken2 can assign a taxonomy in sequence at interstellar speed 🚀.

The main problem we solved here is taking a DNA sequence of unknown origin and quickly comparing it to a large number of other labeled sequences that are identical and incomplete. This is almost exactly the same as searching for nearby web pages or finding cheats. So it is not surprising to see techniques like minhashing and location sensitive hashing appear in these bioinformatics tools like kraken2 and sourmash.

Traditional tools for comparing biological sequences such as BLAST were mainly focused on the evolutionary perspective but as the amount of sequence data available continues to increase bioinformatics experts have had to turn to faster methods that do not focus too much on the biological origin of the sequences.

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