# How do I efficiently perform a metagenome screen of “all” species?

I’ve got an RNA-seq dataset with a large proportion of environmental RNA “contamination”. BLASTing random reads reveals that much of the data comes from bacterial, plant and viral RNA. My target organism only accounts for ~5% of the RNA-seq read data.

I would like to get a more or less (more, if possible) comprehensive overview of what species are found in my sample — bacteria, plants, animals (?), but also viruses. How can I perform this?

I have been unable to find a “standard” way of performing a metagenome screen. BLAST online services all seem to be severely rate limited (certainly unable to handle uploading ~80 M reads). Installing BLAST (or a similar tool) locally is of course not an issue, but I can’t find a comprehensive database that spans all phyla — the best I’ve been able to find are databases, such as NCBI-NR, that are restricted to single phyla or classes, or to bacteria.

Ideally I’d like a workflow that I feed with my RNA-seq data, and get out a list of species with % of reads mapped.

• NCBI's nr database isn't restricted to any species as far as I know. Why do you say it's limited to single phyla or classes or to bacteria? Wouldn't blasting against nr be good enough? And are they really that rate limited? I admit I haven't done this in a few years, but back in the day, I could use NCBI's command-line blast client which would use a remote nr DB and could happily handle several thousand sequences. Have you tried that? Perhaps by splitting your input reads? Aug 1 '17 at 16:46
• @terdon From what I’ve read, it’s only for procaryotes. But I may be completely wrong. Regarding the rate limit, the NCBI BLAST FAQ states that they don’t support batch usage, and that large jobs shouldn’t be run through the URL API. Aug 1 '17 at 16:50
• Nr covers a wide range of taxa as @terdon says. You can download it and build your own blast db if you have ~100gb disk space and a medium memory server. However, BLAST will be slow, I would recommend DIAMOND instead Aug 1 '17 at 19:16

As you point out, BLAST is not going to work well for reads from a typical RNASeq experiment.

From an RNASeq dataset, you might be able to get away with just a ribosomal screen using Silva, because the ribosome is the most abundantly transcribed gene. This will usually be the case even if polyA selection or ribosomal depletion was carried out on a sample.

However, bear in mind that there will be a lot of false results due to database incompleteness. I used Kraken when I did something like this a few years ago, and ended up with one sequence annotated as a fruit bat (where my subject was a water-living planarian flatworm). Centrifuge would be a better option now; the authors even claim that they can store the entirety of nr in a 70 GB index with Centrifuge (which would be better / more comprehensive than just a ribosomal screen).

I think you're just looking for nr. It is absolutely not limited to prokaryotes, far from it. It is, according to the blurb on the NCBI's blast page:

The nucleotide collection consists of GenBank+EMBL+DDBJ+PDB+RefSeq sequences, but excludes EST, STS, GSS, WGS, TSA, patent sequences as well as phase 0, 1, and 2 HTGS sequences. The database is non-redundant. Identical sequences have been merged into one entry, while preserving the accession, GI, title and taxonomy information for each entry.

The link you gave in your question doesn't describe the nr database, but the collection of RefSeq non-redundant sequences. The two are not the same thing at all. The nr databases contains sequences from all available species so it should be perfect for what you're after.

Now, I admit I haven't tried this with millions of redas, but I did used to run many blast searches with several thousand sequences using the command line blast tools and the remote nr database. I just checked the blast+ docs now and see no mention of a rate limit. Which isn't to say there isn't one, for sure, but since each query is run sequentially when you give a multifasta file, there seems to reason to limit since the servers should be able to queue the jobs as needed.

So, don't use the API or roll your own, just download NCBI's blast client and try that:

blastn -query query.fa -db nr -task blastn


If that complains, you'll have to download the nr DB (well, since you want nucleotide blast, you probably want its nucleotide cousin, nt) which you can find here: ftp://ftp.ncbi.nlm.nih.gov/blast/db/. Just get all the foo.nt* files.

• According to the BLAST API Usage guidelines, this isn’t really supported, even if there is no hard limit. But you seem to be right about NR, this seems to be what I’m after, after all. Aug 1 '17 at 17:58
• @KonradRudolph I'm hoping that refers to using the API in your own tools and not to using their own standalone blast+ program. I may very well be wrong, but I'm not sure if the instructions for the API apply to them as well. I would expect them to have built their tool in such a way as to not overload the servers. I figure it's at least worth a try. Aug 1 '17 at 18:27

You might want to look into minihashes as a more efficient way of clustering and comparing your sequences to known species and samples. In the mash paper the authors show how it would be feasible to create minihash sketches that you could search against from large metagenomics datasets (check the "Clustering massive metagenomic datasets" section in the results):

For a large-scale test, samples from the Human Microbiome Project [36] (HMP) and Metagenomics of the Human Intestinal Tract [37] (MetaHIT) were combined to create a ~10 TB 888-sample dataset. Importantly, the size of a Mash sketch is independent of the input size, requiring only 70 MB to store the combined sketches (s = 10,000, k = 21) for these datasets.

Moreover, the creators of sourmash (a tool using the same approach) provide precomputed sketches for the whole of NCBI's RefSeq database (includes bacterial, viral and fungal genomes). Creating a similar sketch for the whole of ENA/SRA (or a subset of it) could require some time but would allow for very fast and accurate searches.

One thing that hasn't been suggested so far would be phylogenetic placement. First option there would be placing against one larger tree roughly touching the domains you mentioned to get an overview. Second option would be having three trees (bacteria, plants, viruses) made up of relevant taxa for the type of sample you have.

This would give you a detailed and accurate view of your sample, but it does involve some work. Sadly the computational throughput for placement is still a bit behind BLAST, and not many complete pipelines exist, but in my opinion it is worth the extra mile (I might be biased. I'm working to resolve as many of these issues as possible).

If you want to look into it, I can reccomend this excellent paper by Erick Matsen. Since your samples sound vaguely soil-ish, I will also link this recent Nature paper using placement to analyze rainforest soil samples.

I have had good success using the MEGAN software to detect and visualize contamination. My input data was WGS not RNA-seq But I would check it out. It can create a very nice weight phylogenetic tree for all the detected organisms.

https://ab.inf.uni-tuebingen.de/software/megan6

• I've downloaded MEGAN but my (admittedly cursory) reading of the documentation made it seem as if the software only visualised results computed by other tools. Aug 3 '17 at 22:32