Answers.Body: <p>SILVA is intended for metagenomic samples and thus does not contain the human ribosomal subunits.  While the 5.8S subunit is in the ensembl, ncbi, and ucsc annotations, and there are numerous rRNA "biotype" (or "biomol") entries, the 5S, 18S, and 28S subunit sequences are absent.  Even downloading rRNA entries in the rmsk table from ucsc will not return all 4 major human subunits.</p>
  
  <p>You can however (as I discovered) download all rDNA gene sequences (including those 4 subunits) from the nucleotide db on ncbi using the search term "biomol_rrna[PROP]".  This returns 240,240 coding sequences (or 30,775 gene sequences).  You can then use this as a reference for mapping, blasting, or as input to rRNA QC tools.</p>
  
  <p>The following e-Utilities commands from ncbi will download the gene sequences:</p>
  
  <pre><code>esearch -db nuccore -query "biomol_rRNA[prop]" | efetch -format fasta &gt; all_rrnas.fasta
  </code></pre>
  
  <p>This is not a common analysis.  Usually, rRNA reads can just be ignored and analyses are unaffected, which is likely why there is a paucity of records for these sequences in annotation files, but mapping to just these sequences could make for a simple QC measure as a gauge for rRNA contamination.  While the sequences may be present in the genome, they're not in the annotations and may even be masked as repetitive, since there are so many copies.  So you have to be careful to count unique reads and not mappings.</p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-01-15T19:57:30.487Z
Answers.LastActivityDate: 2020-01-15T19:57:30.487Z
Answers.UserDisplayName: hepcat72
Answers.Body: <p>I'd recommend using <a href="https://ccb.jhu.edu/software/kraken2/" rel="nofollow noreferrer">kraken2</a> for this. There's a pre-built kraken2 database for SILVA (16S only) which can be downloaded from here:</p>
  
  <p><a href="https://ccb.jhu.edu/software/kraken2/downloads.shtml" rel="nofollow noreferrer">https://ccb.jhu.edu/software/kraken2/downloads.shtml</a></p>
  
  <p>For other ribosomal subunits, it's necessary to create a kraken2 database from scratch. 16S/18S subunits are in the SSU database, while 23S/28S are in the LSU database:</p>
  
  <p><a href="https://www.arb-silva.de/download/arb-files/" rel="nofollow noreferrer">https://www.arb-silva.de/download/arb-files/</a></p>
  
  <p>Note that SILVA is a database of collected rRNA sequences. As far as I'm aware there's no canonical "human" 18S or 28S sequence in SILVA, but it includes a collection of multiple 18S &amp; 28S sequences that have been annotated with the human taxon.</p>
  
  <p>Kraken2 has a command-line interface for creating reports:</p>
  
  <pre><code>kraken2 --db &lt;location_of_database&gt; --gzip-compressed input.fastq.gz --report output.kraken.report.txt &gt; output.kraken.txt
  </code></pre>
  
  <p><a href="https://ccb.jhu.edu/software/pavian/" rel="nofollow noreferrer">Pavian</a> can be used for visualising the results and creating proportion tables.</p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-01-15T22:05:03.817Z
Answers.Text: But SILVA doesn’t have the 5S, 18S, or 28S subunits in it. Any tool suggestion that uses the SILVA database is going to have the same problem. Besides, the question is about what databases of rRNAs to use, not about tools. I can use any tool to count rRNA reads and do other QC, but regardless of the tool: garbage in, garbage out. A tool can’t fill in missing data and SILVA simply doesn’t have all the human rDNA sequences.
Answers.UserDisplayName: hepcat72
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-01-16T02:14:40.540Z
Answers.Text: Does what? I’ve downloaded all the SULVA Fastq dukes and blasted all 4 subunits against it. They’re not there. I don’t think they’re intended to be there. If they’re there, give me the accessions.
Answers.UserDisplayName: hepcat72
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-01-16T04:32:31.697Z
Answers.Text: Autocorrect. “Silva fasta files”. Not SULVA fastq dukes”.
Answers.UserDisplayName: hepcat72
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-01-16T04:38:08.500Z
Answers.Text: No, I mean, show me the SILVA accession number for the specific human 5S, 18S, or 28S ribosomal subunits. They’re not there. I have the sequences for them. I blasted the SILVA fasta files with them. They’re simply not there. And when I use the SILVA fasta file as a reference, I do not get as many reads mapping to that reference as I do when I use the 5S, 5.8S, 18S, and 28S human ribosomal subunits that I manually retrieved from GenBank as the reference.
Answers.UserDisplayName: hepcat72
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-01-15T21:57:58.117Z
Answers.LastActivityDate: 2020-01-16T04:37:25.810Z
Answers.LastEditDate: 2020-01-16T04:37:25.810Z
Answers.UserDisplayName: gringer
Answers.Body: <p>SortMeRNA has fasta files from RFAM and Silva for 5S, 5.8S, 18S, and 28S.  You can download them from <a href="https://github.com/biocore/sortmerna/tree/master/data/rRNA_databases" rel="nofollow noreferrer">GitHub</a></p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-01-16T22:57:04.560Z
Answers.Text: I did see that there is a galaxy tool wrapper for SortMeRNA that can filter rRNAs. I haven’t tried it yet. I also had quickly looked through rfam, but I didn’t spend time looking to see what it had for download. I assume the tool doesn’t include the data? I’d be interested to see how the SILVA/rfam data combo compares to my 30k rDNA download from NCBI nuccore...
Answers.UserDisplayName: hepcat72
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-01-17T18:54:18.133Z
Answers.Text: There does not appear to be a convenient way of downloading only rRNA entries from RFAM.  They appear to have gene sequences in fasta format for roughly 240k rRNAs, but neither their ftp site nor browsing interface provide a way of compiling and downloading them specifically.
Answers.UserDisplayName: hepcat72
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-01-16T21:02:05.700Z
Answers.LastActivityDate: 2020-01-16T21:02:05.700Z
Answers.UserDisplayName: Chris S.
Body: <p>I would like to create a simple QC check of RNA-Seq data that simply maps the data to rRNA and then counts the number of reads that map.</p>
  
  <p>I've done this manually for human using sequences for 5S, 5.8S, 18S, and 28S rRNA (rDNA, actually).  I gathered those sequences manually to create a small fasta file to use as a reference and it worked well.</p>
  
  <p>I would like to have a more definitive source for this reference though, and I don't necessarily want to use a species-specific database, so I thought that using a fasta file downloaded from <a href="https://www.arb-silva.de/no_cache/download/archive/release_138/Exports/" rel="nofollow noreferrer">SILVA</a> would be great, but I don't know if I'm using the wrong fasta file or what.  I'm not very familiar with SILVA, but I noticed that when I tested it using the same sequencing data, I got unexpectedly many fewer hits than I did with my manually created 4-sequence reference than with the SILVA NR99 fasta file.</p>
  
  <p>So I blasted my 4 rDNA sequences against the SILVA NR99 data and only the 18S was found.  The others (5S, 5.8S, and 28S) don't appear to be there.  Am I misunderstanding the nature of the SILVA project or could I be doing something wrong?</p>
  
  <p><em>UPDATE</em>: Please note, this question is about <strong>where to find rRNA and/or rDNA sequences to use as a reference</strong>, not about tools that can be employed in quantifying rRNA read content.  While I appreciate hearing about tools that can generate nice reports about rRNA reads in an RNA-Seq run, the question is not answered unless it addresses the missing rRNA/rDNA sequences that are necessary for comprehensive quantification.  Any tool can count rRNA/rDNA mappings, but the results will only be as good as the rRNA reference data supplied (unless it is a reference-less tool).  And it's the lack of 5S, 18S, and 28S ribosomal subunits in SILVA that is my main concern.  It's the difference between 1k reads and 500k reads mapping to rRNA/rDNA sequences.</p>
  
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2020-01-15T13:37:58.870Z
Comments.Text: Have you thought about using an existing tool to do this QC? RNAseQC (https://github.com/broadinstitute/rnaseqc) for example has has a line that tells you how many reads mapped to rRNA genes.
Comments.UserDisplayName: PPK
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2020-01-15T14:40:25.423Z
Comments.Text: Using what reference & annotations? The 5S, 18S, and 28S are nowhere to be found. I couldn’t find them in any GTF files. I’m not an expert in rRNA, but I assumed they would be there. The researchers complained they weren’t, & when I looked, I was shocked I could not find them. I had to gather them up manually. I didn’t know about that output of rnaseqc, but I assume it’s dependent on the annotation you provide it.
Comments.UserDisplayName: hepcat72
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2020-01-15T16:36:56.103Z
Comments.Text: I was using bwa in my galaxy workflow, but any relevant tool is fine. rnaseqc has a `-BWArRNA <arg>` which takes a fasta reference for rRNA sequences. It mentions that if that flag is absent, it uses the GTF file, which is says is fast, but not robust, which I assume alludes to the fact that the GTF files don't have the rRNA gene entries necessary, which is the whole point of this question. When I simply looked at the rRNA entries that are in the gtf for hg38p12, there were very few reads reported, but against my 4 rRNA seq ref, there were 5 thousand times more.
Comments.UserDisplayName: hepcat72
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2020-01-15T16:38:43.490Z
Comments.Text: And BTW, since I'm using galaxy, do you know of an rnaseqc wrapper that includes that -BWArRNA flag? I can't find one that includes it.
Comments.UserDisplayName: hepcat72
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2020-01-15T17:13:20.503Z
Comments.Text: As far as I understand it uses the genes with the "rRNA" biotype in the Ensembl gene annotation.
Comments.UserDisplayName: PPK
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2020-01-15T17:55:33.060Z
Comments.Text: The doc says it uses the supplied GTF. Is the documentation wrong or am I misinterpreting it?
Comments.UserDisplayName: hepcat72
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2020-01-15T19:25:37.837Z
Comments.Text: I just double-checked, and while there are a number of 5.8S entries in the ensemble rRNA biotype annotations for human, none are for the 5S, 18S, or 28S subunits, which is the problem I cited in my question. I tried out SILVA to try and find those rRNAs, not only for human, but for all species.
Comments.UserDisplayName: hepcat72
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2020-01-16T04:46:11.763Z
Comments.Text: You appear to misunderstand the purpose of SILVA. It is "a comprehensive on-line resource for quality checked and aligned ribosomal RNA sequence data." It contains multiple rRNA sequences that have been assigned to the human taxon. As far as I am aware, it doesn't contain consensus sequences for any taxa.
Comments.UserDisplayName: gringer
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2020-01-16T04:51:27.940Z
Comments.Text: And that’s exactly what my question was - whether I was misunderstanding SILVA and if I could use it to count rRNA sequences in an RNA-Seq run. The answer is yes, I was misunderstanding it and no, you can’t use it to precisely quantify specific species’ rRNA contamination.
Comments.UserDisplayName: hepcat72
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2020-01-16T04:53:25.397Z
Comments.Text: The SILVA NR99 data only contains 16S/18S sequences. If you want 23S/28S sequences, then try the LSU data from release 132:
  
  https://www.arb-silva.de/no_cache/download/archive/release_132/Exports/
Comments.UserDisplayName: gringer
Tags: <ribosomal>
Title: Ribosomal RNA QC quantification using SILVA
ThreadText: # Ribosomal RNA QC quantification using SILVA
  
  (Asked by: username_0 on 2020-01-13)
  
  I would like to create a simple QC check of RNA-Seq data that simply maps the data to rRNA and then counts the number of reads that map.
  
  I've done this manually for human using sequences for 5S, 5.8S, 18S, and 28S rRNA (rDNA, actually). I gathered those sequences manually to create a small fasta file to use as a reference and it worked well.
  
  I would like to have a more definitive source for this reference though, and I don't necessarily want to use a species-specific database, so I thought that using a fasta file downloaded from [SILVA](https://www.arb-silva.de/no_cache/download/archive/release_138/Exports/) would be great, but I don't know if I'm using the wrong fasta file or what. I'm not very familiar with SILVA, but I noticed that when I tested it using the same sequencing data, I got unexpectedly many fewer hits than I did with my manually created 4-sequence reference than with the SILVA NR99 fasta file.
  
  So I blasted my 4 rDNA sequences against the SILVA NR99 data and only the 18S was found. The others (5S, 5.8S, and 28S) don't appear to be there. Am I misunderstanding the nature of the SILVA project or could I be doing something wrong?
  
  *UPDATE*: Please note, this question is about **where to find rRNA and/or rDNA sequences to use as a reference**, not about tools that can be employed in quantifying rRNA read content. While I appreciate hearing about tools that can generate nice reports about rRNA reads in an RNA-Seq run, the question is not answered unless it addresses the missing rRNA/rDNA sequences that are necessary for comprehensive quantification. Any tool can count rRNA/rDNA mappings, but the results will only be as good as the rRNA reference data supplied (unless it is a reference-less tool). And it's the lack of 5S, 18S, and 28S ribosomal subunits in SILVA that is my main concern. It's the difference between 1k reads and 500k reads mapping to rRNA/rDNA sequences.
  
  **Tags:** ribosomal
  
  > **Comments:**
  > **username_1**: Have you thought about using an existing tool to do this QC? RNAseQC (https://github.com/broadinstitute/rnaseqc) for example has has a line that tells you how many reads mapped to rRNA genes.  
  > **username_0**: Using what reference & annotations? The 5S, 18S, and 28S are nowhere to be found. I couldn’t find them in any GTF files. I’m not an expert in rRNA, but I assumed they would be there. The researchers complained they weren’t, & when I looked, I was shocked I could not find them. I had to gather them up manually. I didn’t know about that output of rnaseqc, but I assume it’s dependent on the annotation you provide it.  
  > **username_0**: I was using bwa in my galaxy workflow, but any relevant tool is fine. rnaseqc has a `-BWArRNA <arg>` which takes a fasta reference for rRNA sequences. It mentions that if that flag is absent, it uses the GTF file, which is says is fast, but not robust, which I assume alludes to the fact that the GTF files don't have the rRNA gene entries necessary, which is the whole point of this question. When I simply looked at the rRNA entries that are in the gtf for hg38p12, there were very few reads reported, but against my 4 rRNA seq ref, there were 5 thousand times more.  
  > **username_0**: And BTW, since I'm using galaxy, do you know of an rnaseqc wrapper that includes that -BWArRNA flag? I can't find one that includes it.  
  > **username_1**: As far as I understand it uses the genes with the "rRNA" biotype in the Ensembl gene annotation.  
  > **username_0**: The doc says it uses the supplied GTF. Is the documentation wrong or am I misinterpreting it?  
  > **username_0**: I just double-checked, and while there are a number of 5.8S entries in the ensemble rRNA biotype annotations for human, none are for the 5S, 18S, or 28S subunits, which is the problem I cited in my question. I tried out SILVA to try and find those rRNAs, not only for human, but for all species.  
  > **username_2**: You appear to misunderstand the purpose of SILVA. It is "a comprehensive on-line resource for quality checked and aligned ribosomal RNA sequence data." It contains multiple rRNA sequences that have been assigned to the human taxon. As far as I am aware, it doesn't contain consensus sequences for any taxa.  
  > **username_0**: And that’s exactly what my question was - whether I was misunderstanding SILVA and if I could use it to count rRNA sequences in an RNA-Seq run. The answer is yes, I was misunderstanding it and no, you can’t use it to precisely quantify specific species’ rRNA contamination.  
  > **username_2**: The SILVA NR99 data only contains 16S/18S sequences. If you want 23S/28S sequences, then try the LSU data from release 132:  https://www.arb-silva.de/no_cache/download/archive/release_132/Exports/  
  
  ---
  
  ## 3 answers
  
  ### Answer by username_0 on 2020-01-15. Score: 1 | Accepted answer
  
  SILVA is intended for metagenomic samples and thus does not contain the human ribosomal subunits. While the 5.8S subunit is in the ensembl, ncbi, and ucsc annotations, and there are numerous rRNA "biotype" (or "biomol") entries, the 5S, 18S, and 28S subunit sequences are absent. Even downloading rRNA entries in the rmsk table from ucsc will not return all 4 major human subunits.
  
  You can however (as I discovered) download all rDNA gene sequences (including those 4 subunits) from the nucleotide db on ncbi using the search term "biomol\_rrna[PROP]". This returns 240,240 coding sequences (or 30,775 gene sequences). You can then use this as a reference for mapping, blasting, or as input to rRNA QC tools.
  
  The following e-Utilities commands from ncbi will download the gene sequences:
  
  ```
  esearch -db nuccore -query "biomol_rRNA[prop]" | efetch -format fasta > all_rrnas.fasta
  
  ```
  
  This is not a common analysis. Usually, rRNA reads can just be ignored and analyses are unaffected, which is likely why there is a paucity of records for these sequences in annotation files, but mapping to just these sequences could make for a simple QC measure as a gauge for rRNA contamination. While the sequences may be present in the genome, they're not in the annotations and may even be masked as repetitive, since there are so many copies. So you have to be careful to count unique reads and not mappings.
  
  
  ---
  
  ### Answer by username_2 on 2020-01-15. Score: 0
  
  I'd recommend using [kraken2](https://ccb.jhu.edu/software/kraken2/) for this. There's a pre-built kraken2 database for SILVA (16S only) which can be downloaded from here:
  
  <https://ccb.jhu.edu/software/kraken2/downloads.shtml>
  
  For other ribosomal subunits, it's necessary to create a kraken2 database from scratch. 16S/18S subunits are in the SSU database, while 23S/28S are in the LSU database:
  
  <https://www.arb-silva.de/download/arb-files/>
  
  Note that SILVA is a database of collected rRNA sequences. As far as I'm aware there's no canonical "human" 18S or 28S sequence in SILVA, but it includes a collection of multiple 18S & 28S sequences that have been annotated with the human taxon.
  
  Kraken2 has a command-line interface for creating reports:
  
  ```
  kraken2 --db <location_of_database> --gzip-compressed input.fastq.gz --report output.kraken.report.txt > output.kraken.txt
  
  ```
  
  [Pavian](https://ccb.jhu.edu/software/pavian/) can be used for visualising the results and creating proportion tables.
  
  > **Comments:**
  > **username_0**: But SILVA doesn’t have the 5S, 18S, or 28S subunits in it. Any tool suggestion that uses the SILVA database is going to have the same problem. Besides, the question is about what databases of rRNAs to use, not about tools. I can use any tool to count rRNA reads and do other QC, but regardless of the tool: garbage in, garbage out. A tool can’t fill in missing data and SILVA simply doesn’t have all the human rDNA sequences.  
  > **username_0**: Does what? I’ve downloaded all the SULVA Fastq dukes and blasted all 4 subunits against it. They’re not there. I don’t think they’re intended to be there. If they’re there, give me the accessions.  
  > **username_0**: Autocorrect. “Silva fasta files”. Not SULVA fastq dukes”.  
  > **username_0**: No, I mean, show me the SILVA accession number for the specific human 5S, 18S, or 28S ribosomal subunits. They’re not there. I have the sequences for them. I blasted the SILVA fasta files with them. They’re simply not there. And when I use the SILVA fasta file as a reference, I do not get as many reads mapping to that reference as I do when I use the 5S, 5.8S, 18S, and 28S human ribosomal subunits that I manually retrieved from GenBank as the reference.  
  
  ---
  
  ### Answer by username_3 on 2020-01-16. Score: 0
  
  SortMeRNA has fasta files from RFAM and Silva for 5S, 5.8S, 18S, and 28S. You can download them from [GitHub](https://github.com/biocore/sortmerna/tree/master/data/rRNA_databases)
  
  > **Comments:**
  > **username_0**: I did see that there is a galaxy tool wrapper for SortMeRNA that can filter rRNAs. I haven’t tried it yet. I also had quickly looked through rfam, but I didn’t spend time looking to see what it had for download. I assume the tool doesn’t include the data? I’d be interested to see how the SILVA/rfam data combo compares to my 30k rDNA download from NCBI nuccore...  
  > **username_0**: There does not appear to be a convenient way of downloading only rRNA entries from RFAM.  They appear to have gene sequences in fasta format for roughly 240k rRNAs, but neither their ftp site nor browsing interface provide a way of compiling and downloading them specifically.  
  
  ---
  
  