Monthly Archives: October 2022

mmseqs2 logo

Splitting the Non-Redundant Database for MMSeqs2

We have a computational problem searching the Non-Redundant database, but we can solve that!

If you use MMSeqs2 to search the NR database, it needs about 1.75TB of RAM (that is approximate!). Our Flinders deepthought cluster has several nodes with 2.5TB RAM so at most you can run three concurrent jobs searching the database. However, database searches are completely independent: the scores you get for one pair of sequences are independent from all the other sequences in the database.

Note: That is true, but there is a caveat. The E-value of the search is dependent on the length of the query and the size of the database.

E value depends on the database size!

Therefore, the E-value that is reported is only approximate. There is a GitHub Issue to resolve this with MMSeqs2

Our solution is to split the database into smaller pieces, and then we can use those across more nodes of the cluster.

Download the NR database

Step 1, download the preformatted NR database using mmseqs2

mkdir --parents NR
mmseqs databases --threads 8 NR NR/NR (mktemp -d)

This will download the non-redundant database into the directory NR and the database will be called NR.

Split the database

Let’s split that database into many smaller chunks. From my tests it appears that 50 chunks requires about 50-60 GB RAM per compute, and 100 chunks requires about 25 GB RAM per compute.

CHUNKS=100
mkdir NR.split.$CHUNKS
cd NR.split.$CHUNKS
cut -f 1 ../NR/NR.index > ids.txt
split -n l/$CHUNKS --numeric-suffixes=1 --suffix-length=3 ids.txt
for N in $(seq 1 $CHUNKS); do
        echo $N;
        X="000$N";
        X="x${X:(-3)}";
        mkdir NR.$N
        mmseqs createsubdb $X ../NR/NR NR.$N/NR;
        mmseqs createsubdb $X ../NR/NR_h NR.$N/NR_h;
done
rm -f x???

(Hint: this works quite well as a slurm script!)

There are now 1 directory called NR.split.100 and within that there are 100 directories called NR.1, NR.2, NR.3, … NR.100.

Search against all the directories

We can now use a slurm array job to search against all of those:

#!/bin/bash
###############################################################
#                                                             #
# Search the NR using an array job.                           #
#                                                             #
# submit with this command:                                   #
#                                                             #
# sbatch --array=1-100:1 ./nr_chunks.slurm                    #
#                                                             #
# Written by Rob Edwards. Good luck!                          #
#                                                             #
###############################################################

#SBATCH --job-name=NRmmseqs
#SBATCH --time=5-0
#SBATCH --ntasks=1
#SBATCH --cpus-per-task=24
#SBATCH --mem=25G
#SBATCH -o mmseqsNRA_%a.out
#SBATCH -e mmseqsNRA_%a.err

## Memory requirement: for 100 chunks 25G (average from 72 searches: 16.3889, range: 12-20G)

set -euo pipefail

# change this to point to the output from `mmseqs createdb` on your data:
QUERY='mmseqs_formatted_data.db'

NR=$SLURM_ARRAY_TASK_ID
CHUNKS=100
mmseqs search --threads 24 $QUERY NR.split.$CHUNKS/NR.$NR/NR $QUERY.nr.$CHUNKS.$NR.db $(mktemp -d)
mmseqs convertalis --threads 24 $QUERY NR.split.$CHUNKS/NR.$NR/NR $QUERY.nr.$CHUNKS.$NR.db $QUERY.nr.$CHUNKS.$NR.m8

You can save this file as nr_chunks.slurm, edit the line with the QUERY location, and submit 100 MMSeqs2 jobs as an array job (using sbatch --array=1-100:1 ./nr_chunks.slurm). You now have 2,400 processors computing on your search, and each 24 processors is only consuming 25 GB RAM.

Note: I recommend that you adjust the number of threads, I am not sure if 24 is optimum!

Identifying Metagenomes from the SRA in the Cloud

PARTIE

We have several projects that look for all the metagenomes in the cloud, and we have several ways of searching the SRA. Here, we’ll search for all the WGS metagenomes in the SRA using a Google Big Table query.

Log into Google Console

You’ll need to log into Google console and access a project or create a new one.

Use SQL to find the metagenome/microbiome/metatranscriptome results

We use temporary tables to store the two main searches: what are amplicon projects and what are metagenome/microbiome/metatranscriptome projects, and then we find the projects that are metagenomes:

create temp table AMPLICON(acc STRING) as select acc as amplicon from `nih-sra-datastore.sra.metadata` where assay_type = 'AMPLICON' or libraryselection = 'PCR';
create temp table METAGENOMES(acc STRING) as select acc from `nih-sra-datastore.sra.metadata` where librarysource = "METAGENOMIC" or librarysource = 'METATRANSCRIPTOMIC' or organism like "%microbiom%" OR organism like "%metagenom%"  or organism like '%metatran%';
select acc from METAGENOMES where acc not in (select acc from AMPLICON);

Then save that as a JSON file to Google Drive.

Use jq to parse the JSON file

This is probably overkill because we only have one attribute in our data.

jq -r '.acc' bq-results-20221006-054328-1665035790273.json > SRA-metagenomes.txt

Find all the information about all the sequences

We can edit the above SQL to get all the information about all the metagenomes. Basically, we just change the second select statement.

create temp table AMPLICON(acc STRING) as select acc as amplicon from `nih-sra-datastore.sra.metadata` where assay_type = 'AMPLICON' or libraryselection = 'PCR';
select * from `nih-sra-datastore.sra.metadata` where acc not in (select acc from AMPLICON) and (librarysource = "METAGENOMIC" or librarysource = 'METATRANSCRIPTOMIC' or organism like "%microbiom%" OR organism like "%metagenom%");

Note: In this query, the parenthesis are important to make sure we do the and and or in the right place.

Then you can export the data as a JSON Newline file to Google Drive.

Current results

At the moment, this returns 642,842 runs from the SRA

Some things we can’t find

  • The old study_type field that we searched (using study_type = "Metagenomics") does not appear to have mapped to bigtable.
  • THe old scientific name that we searched (using sample.scientific_name like "%microbiom%" OR sample.scientific_name like "%metagenom%") does not appear to have mapped to bigtable.