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

Question

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.

Answer 1

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).

Answer 2

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.

Answer 3

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.

Answer 4

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.

Answer 5

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