<p>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.</p>

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

<p>If you want to look into it, I can reccomend <a href="https://academic.oup.com/sysbio/article/64/1/e26/2847641" rel="nofollow noreferrer">this excellent paper by Erick Matsen</a>. Since your samples sound vaguely soil-ish, I will also link <a href="https://www.nature.com/articles/s41559-017-0091" rel="nofollow noreferrer">this recent Nature paper</a> using placement to analyze rainforest soil samples.</p>

CC BY-SA 3.0
2017-09-15T14:10:56.140Z
2017-09-15T14:10:56.140Z
Pierre Barbera
<p>I think you're just looking for <code>nr</code>. It is absolutely not limited to prokaryotes, far from it. It is, according to the blurb on the NCBI's blast page:</p>

<blockquote>
  <p>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.</p>
</blockquote>

<p>The link you gave in your question doesn't describe the <code>nr</code> 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. </p>

<p>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 <a href="https://www.ncbi.nlm.nih.gov/books/NBK279690/" rel="nofollow noreferrer">blast+ docs</a> 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. </p>

<p>So, don't use the API or roll your own, just <a href="https://blast.ncbi.nlm.nih.gov/Blast.cgi?CMD=Web&amp;PAGE_TYPE=BlastDocs&amp;DOC_TYPE=Download" rel="nofollow noreferrer">download NCBI's blast client</a> and try that:</p>

<pre><code>blastn -query query.fa -db nr -task blastn 
</code></pre>

<p>If that complains, you'll have to download the <code>nr</code> DB (well, since you want nucleotide blast, you probably want its nucleotide cousin, <code>nt</code>) which you can find here: <a href="ftp://ftp.ncbi.nlm.nih.gov/blast/db/" rel="nofollow noreferrer">ftp://ftp.ncbi.nlm.nih.gov/blast/db/</a>. Just get all the <code>foo.nt*</code> files.</p>

CC BY-SA 3.0
2017-08-01T17:58:48.210Z
According to the [BLAST API Usage guidelines](https://blast.ncbi.nlm.nih.gov/Blast.cgi?CMD=Web&PAGE_TYPE=BlastDocs&DOC_TYPE=DeveloperInfo), 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.
Konrad Rudolph
CC BY-SA 3.0
2017-08-01T18:27:24.883Z
@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.
terdon
CC BY-SA 3.0
2017-08-01T17:36:43.390Z
2017-08-01T17:36:43.390Z
terdon
<p>You might want to look into <a href="https://en.wikipedia.org/wiki/MinHash" rel="nofollow noreferrer">minihashes</a> as a more efficient way of clustering and comparing your sequences to known species and samples. <a href="https://genomebiology.biomedcentral.com/articles/10.1186/s13059-016-0997-x" rel="nofollow noreferrer">In the mash paper</a> 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):</p>

<blockquote>
  <p>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.</p>
</blockquote>

<p>Moreover, the creators of <a href="https://sourmash.readthedocs.io/en/latest/" rel="nofollow noreferrer">sourmash</a> (a tool using the same approach) provide <a href="https://sourmash.readthedocs.io/en/latest/databases.html" rel="nofollow noreferrer">precomputed sketches for the whole of NCBI's RefSeq database</a> (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.</p>

CC BY-SA 3.0
2017-08-03T21:38:09.977Z
2017-08-03T21:38:09.977Z
mgalardini
<p>As you point out, BLAST is not going to work well for reads from a typical RNASeq experiment.</p>

<p>From an RNASeq dataset, you might be able to get away with just a ribosomal screen using <a href="https://www.arb-silva.de/projects/ssu-ref-nr/" rel="noreferrer">Silva</a>, 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.</p>

<p>However, bear in mind that there will be a lot of false results due to database incompleteness. I used <a href="http://ccb.jhu.edu/software/kraken/" rel="noreferrer">Kraken</a> 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). <a href="https://github.com/infphilo/centrifuge/" rel="noreferrer">Centrifuge</a> would be a better option now; the authors even claim that <a href="http://genome.cshlp.org/content/26/12/1721.full" rel="noreferrer">they can store the entirety of nr in a 70 GB index</a> with Centrifuge (which would be better / more comprehensive than just a ribosomal screen).</p>

CC BY-SA 3.0
2017-08-01T21:54:51.507Z
2017-08-01T21:54:51.507Z
gringer
<p>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.</p>

<p><a href="https://ab.inf.uni-tuebingen.de/software/megan6" rel="nofollow noreferrer">https://ab.inf.uni-tuebingen.de/software/megan6</a></p>

<p><a href="https://i.sstatic.net/ttT6l.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/ttT6l.png" alt="enter image description here"></a></p>

CC BY-SA 3.0
2017-08-03T22:32:43.370Z
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.
Konrad Rudolph
CC BY-SA 3.0
2017-08-03T18:55:25.983Z
2017-08-03T18:55:25.983Z
Bioathlete
<p>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.</p>

<p>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?</p>

<p>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 <a href="https://www.ncbi.nlm.nih.gov/refseq/about/nonredundantproteins/" rel="noreferrer">NCBI-NR</a>, that are restricted to single phyla or classes, or to bacteria.</p>

<p>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.</p>

CC BY-SA 3.0
2017-08-01T16:46:41.617Z
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?
terdon
CC BY-SA 3.0
2017-08-01T16:50:45.040Z
@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.
Konrad Rudolph
CC BY-SA 3.0
2017-08-01T19:16:24.470Z
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
Chris_Rands
<rna-seq><metagenome>
How do I efficiently perform a metagenome screen of “all” species?
# How do I efficiently perform a metagenome screen of “all” species?

(Asked by: username_0 on 2017-08-01)

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](https://www.ncbi.nlm.nih.gov/refseq/about/nonredundantproteins/), 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.

**Tags:** rna-seq, metagenome
> **Comments:**
> **username_1**: 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?  
> **username_0**: @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.  
> **username_2**: 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  

---

## 5 answers

### Answer by username_5 on 2017-08-01. Score: 7

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](https://www.arb-silva.de/projects/ssu-ref-nr/), 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](http://ccb.jhu.edu/software/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](https://github.com/infphilo/centrifuge/) would be a better option now; the authors even claim that [they can store the entirety of nr in a 70 GB index](http://genome.cshlp.org/content/26/12/1721.full) with Centrifuge (which would be better / more comprehensive than just a ribosomal screen).


---

### Answer by username_3 on 2017-09-15. Score: 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](https://academic.oup.com/sysbio/article/64/1/e26/2847641). Since your samples sound vaguely soil-ish, I will also link [this recent Nature paper](https://www.nature.com/articles/s41559-017-0091) using placement to analyze rainforest soil samples.


---

### Answer by username_1 on 2017-08-01. Score: 4

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.

> **Comments:**
> **username_0**: According to the [BLAST API Usage guidelines](https://blast.ncbi.nlm.nih.gov/Blast.cgi?CMD=Web&PAGE_TYPE=BlastDocs&DOC_TYPE=DeveloperInfo), 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.  
> **username_1**: @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.  

---

### Answer by username_4 on 2017-08-03. Score: 4

You might want to look into [minihashes](https://en.wikipedia.org/wiki/MinHash) as a more efficient way of clustering and comparing your sequences to known species and samples. [In the mash paper](https://genomebiology.biomedcentral.com/articles/10.1186/s13059-016-0997-x) 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](https://sourmash.readthedocs.io/en/latest/) (a tool using the same approach) provide [precomputed sketches for the whole of NCBI's RefSeq database](https://sourmash.readthedocs.io/en/latest/databases.html) (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 by username_6 on 2017-08-03. Score: 1

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>

[![enter image description here](https://i.sstatic.net/ttT6l.png)](https://i.sstatic.net/ttT6l.png)

> **Comments:**
> **username_0**: 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.  

---

