3
$\begingroup$

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

Improve this question
$\endgroup$
  • $\begingroup$ Have you thought about using an existing tool to do this QC? RNAseQC (github.com/broadinstitute/rnaseqc) for example has has a line that tells you how many reads mapped to rRNA genes. $\endgroup$ – PPK Jan 15 '20 at 13:37
  • $\begingroup$ 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. $\endgroup$ – hepcat72 Jan 15 '20 at 14:40
  • $\begingroup$ 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. $\endgroup$ – hepcat72 Jan 15 '20 at 16:36
  • $\begingroup$ 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. $\endgroup$ – hepcat72 Jan 15 '20 at 16:38
  • $\begingroup$ As far as I understand it uses the genes with the "rRNA" biotype in the Ensembl gene annotation. $\endgroup$ – PPK Jan 15 '20 at 17:13
1
$\begingroup$

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.

Improve this answer
$\endgroup$
0
$\begingroup$

I'd recommend using 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 can be used for visualising the results and creating proportion tables.

Improve this answer
$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – hepcat72 Jan 15 '20 at 22:05
  • $\begingroup$ 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. $\endgroup$ – hepcat72 Jan 16 '20 at 2:14
  • $\begingroup$ Autocorrect. “Silva fasta files”. Not SULVA fastq dukes”. $\endgroup$ – hepcat72 Jan 16 '20 at 4:32
  • $\begingroup$ 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. $\endgroup$ – hepcat72 Jan 16 '20 at 4:38
0
$\begingroup$

SortMeRNA has fasta files from RFAM and Silva for 5S, 5.8S, 18S, and 28S. You can download them from GitHub

Improve this answer
$\endgroup$
  • $\begingroup$ 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... $\endgroup$ – hepcat72 Jan 16 '20 at 22:57
  • $\begingroup$ 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. $\endgroup$ – hepcat72 Jan 17 '20 at 18:54

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.