rnaQUAST is a tool for evaluating RNA-Seq assemblies using reference genome and gene database. In addition, rnaQUAST is also capable of estimating gene database coverage by raw reads and de novo quality assessment using third-party software.
rnaQUAST version 2.3 was released under GPLv2 on June 21st, 2024 and can be downloaded from https://github.com/ablab/rnaquast/releases.
For impatient people:
You will need Python, gffutils, matplotlib and joblib. Also you will need GMAP (or BLAT) and BLASTN installed on your machine and added to the $PATH
variable.
You may also install rnaQUAST via conda
conda install -c bioconda rnaquast
To verify your installation run
python rnaQUAST.py --test
To run rnaQUAST on your data use the following command
python rnaQUAST.py \
--transcripts /PATH/TO/transcripts1.fasta /PATH/TO/ANOTHER/transcripts2.fasta /PATH/TO/MULTIPLE/*.fasta [...] \
--reference /PATH/TO/reference_genome.fasta --gtf /PATH/TO/gene_coordinates.gtf
rnaQUAST can be installed via conda:
conda install -c bioconda rnaquast
If you wish to run rnaQUAST from the release archive you need:
rnaQUAST still works under Python2 (2.5+), but since Python2 is outdated, its support is not maintained since version 2.0.
Note, that due to the limitations of BLAT
, in order to work with reference genomes of size more than 4 Gb a pslSort
is also required.
Paths to blastn
and GMAP
(or BLA
T) should be added to the $PATH
environmental variable. To check that everything is installed correctly we recommend to run:
python rnaQUAST.py --test
Note that gffutils
is used to complete gene coordinates in case of missing transcripts / genes records. For more information, see advanced options.
When reference genome and gene database are unavailable, we recommend to run BUSCO and GeneMarkS-T in rnaQUAST pipeline.
BUSCO requirements
BUSCO allows to detect core genes in the assembled transcripts. To use it you should install BUSCO v4+, tblastn, HMMER and transeq and add these tools to the $PATH
variable.
To run BUSCO provide lineage-specific database name via --busco
option. You may also download the appropriate database from http://busco.ezlab.org manually and provide it using the same option (see options for details).
GeneMarkS-T requirements
GeneMarkS-T allows to predict genes in the assembled transcripts without reference genome. If you wish to use it in rnaQUAST pipeline, GeneMarkS-T should be properly installed and added to the $PATH
variable.
rnaQUAST is also capable of calculating various statistics using raw reads (e.g. database coverage by reads). To obtain them you need to install STAR aligner and add it to the $PATH
variable. To learn more see input options.
To run rnaQUAST you need to provide either FASTA files with transcripts (recommended), or align transcripts to the reference genome manually and provide the resulting PSL files.
-r <REFERENCE>, --reference <REFERENCE>
Single file with reference genome containing all chromosomes/scaffolds in FASTA format (preferably with *.fasta, *.fa, *.fna, *.ffn or *.frn
extension) OR
*`.txt`** file containing the one-per-line list of FASTA files with reference sequences.
--gtf <GENE_COORDINATES>
File with gene coordinates in GTF/GFF format (needs information about parent relations). We recommend to use files downloaded from GENCODE or Ensembl.
--gene_db <GENE_DB>
Path to the gene database generated by gffutils. The database is created during the first run. This option is not compatible with --gtf
option. We recommend to use this option once the database is created in order to speed up the run.
--gmap_index <INDEX FOLDER>,
Folder containing pre-built GMAP index for the reference genome. Using previously constructed index decreases running time. Note, that you still need to provide the reference genome that was used for index construction when this option is used.
-c <TRANSCRIPTS ...>, --transcripts <TRANSCRIPTS, ...>
File(s) with transcripts in FASTA format separated by space. Wildcards can be used, e.g. --transcripts */*.fasta
.
-psl <TRANSCRIPTS_ALIGNMENT ...>, --alignment <TRANSCRIPTS_ALIGNMENT, ...>
File(s) with transcript alignments to the reference genome in PSL format separated by space.
-sam <READS_ALIGNMENT>, --reads_alignment <READS_ALIGNMENT>
File with read alignments to the reference genome in SAM format.
-1 <LEFT_READS>, --left_reads <LEFT_READS>
File with forward paired-end reads in FASTQ or gzip-compressed fastq format.
-2 <RIGHT_READS>, --right_reads <RIGHT_READS>
File with reverse paired-end reads in FASTQ or gzip-compressed fastq format.
-s <SINGLE_READS>, --single_reads <SINGLE_READS>
File with single reads in FASTQ or gzip-compressed fastq format.
-o <OUTPUT_DIR>, --output_dir <OUTPUT_DIR>
Directory to store all results. Default is rnaQUAST_results/results_<datetime>
.
--test
Run rnaQUAST on the test data from the test_data
folder, output directory is rnaOUAST_test_output
.
-d, --debug
Report detailed information, typically used only for detecting problems.
-h, --help
Show help message and exit.
-t <INT>, --threads <INT>
Maximum number of threads. Default is min(number of CPUs / 2, 16).
-l <LABELS ...>, --labels <LABELS ...>
Name(s) of assemblies that will be used in the reports separated by space and given in the same order as files with transcripts / alignments.
--prokaryote
Use this option if the genome is prokaryotic.
-ss, --strand_specific
Set if transcripts were assembled using strand-specific RNA-Seq data in order to benefit from knowing whether the transcript originated from the + or - strand.
--min_alignment <MIN_ALIGNMENT>
Minimal alignment length to be used, default value is 50.
--no_plots
Do not draw plots (makes rnaQUAST run a bit faster).
--blat
Run with BLAT alignment tool instead of GMAP.
--busco
Run BUSCO tool, which detects core genes in the assembly (see Installation & requirements). Use this option to provide BUSCO database name to use or path to the local database. Also, you can set auto-lineage
for automated lineage selection.
--gene_mark
Run with GeneMarkS-T gene prediction tool. Use --prokaryote
option if the genome is prokaryotic.
--disable_infer_genes
Use this option if your GTF file already contains genes records, otherwise gffutils will fix it. Note that gffutils may work for quite a long time.
--disable_infer_transcripts
Is option if your GTF file already contains transcripts records, otherwise gffutils will fix it. Note that gffutils may work for quite a long time.
--lower_threshold
Lower threshold for x-assembled/covered/matched metrics, default: 50%.
--upper_threshold
Upper threshold for x-assembled/covered/matched metrics, default: 95%.
In this section we describe metrics, statistics and plots generated by rnaQUAST. Metrics highlighted with bold italic are considered as the most important and are included in the short summary report (short_report.txt
).
For the simplicity of explanation, transcript is further referred to as a sequence generated by the assembler and isoform denotes a sequence from the gene database. Figure below demonstrates how rnaQUAST classifies transcript and isoform sequences using alignment information.
The following text files with reports are contained in comparison_output
directory and include results for all input assemblies. In addition, these reports are contained in <assembly_label>_output
directories for each assembly separately.
database_metrics.txt
Gene database metrics.
Coverage by reads. The following metrics are calculated only when --left_reads
, --right_reads
, --single_reads
or --sam
options are used (see options for details).
--lower_threshold / --upper_threshold
options (50% / 95% by default).basic_mertics.txt
Basic transcripts metrics are calculated without reference genome and gene database.
alignment_metrics.txt
Alignment metrics are calculated with reference genome but without using gene database. To calculate the following metrics rnaQUAST filters all short partial alignments (see --min_alignment
option) and attempts to select the best hits for each transcript.
<assembly_label>.paralogs.fasta
file.<assembly_label>.unaligned.fasta
file.Number of assembled transcripts = Unaligned + Aligned = Unaligned + (Uniquely aligned + Multiply aligned + Misassembly candidates reported by GMAP (or BLAT)).
Alignment metrics for non-misassembled transcripts
<assembly_label>.misassembled.fasta
file.sensitivity.txt
Assembly completeness (sensitivity). For the following metrics (calculated with reference genome and gene database) rnaQUAST attempts to select best-matching database isoforms for every transcript. Note that a single transcript can contribute to multiple isoforms in the case of, for example, paralogous genes or genomic repeats. At the same time, an isoform can be covered by multiple transcripts in the case of fragmented assembly or duplicated transcripts in the assembly.
--lower_threshold / --upper_threshold
options (50% / 95% by default). 95%-assembled isoforms are stored in <assembly_label>.95%assembled.fasta
file.--lower_threshold / --upper_threshold
options (50% / 95% by default).--lower_threshold / --upper_threshold
options (50% / 95% by default). For each isoform rnaQUAST calculates the number of x%-covered exons divided by the total number of exons. Afterwards it computes average value for all covered isoforms.BUSCO metrics. The following metrics are calculated only when --busco
option is used (see options for details).
GeneMarkS-T metrics. The following metrics are calculated when reference and gene database are not provided or --gene_mark
option is used (see options for details).
specificity.txt
Assembly specificity. To compute the following metrics we use only transcripts that have at least one significant alignment and are not misassembled.
<assembly_label>.unannotated.fasta
file.--lower_threshold / --upper_threshold
options (50% / 95% by default).--lower_threshold / --upper_threshold
options (50% / 95% by default).relative_database_coverage.txt
Relative database coverage metrics are calculated only when raw reads (or read alignments) are provided. rnaQUAST uses read alignments to estimate the upper bound of the database coverage and the number of x-covered genes / isoforms / exons (see read coverage) and computes the following metrics:
These files are contained in <assembly_label>_output
directories for each assembly separately.
<assembly_label>.unaligned.fasta
– transcripts without any significant alignments.<assembly_label>.paralogs.fasta
– transcripts having 2 or more significant alignments.<assembly_label>.misassembled.fasta
– misassembly candidates detected by methods described above. See misassemblies.txt
description for details.<assembly_label>.correct.fasta
– transcripts with exactly 1 significant alignment that do not contain misassemblies.<assembly_label>.x%-assembled.list
– IDs of the isoforms from the database that have at least x% captured by a single assembled transcript, where x is specified by the user with an option --upper_threshold
(95% by default).<assembly_label>.unannotated.fasta
– transcripts that do not cover any isoform from the database.The following text file is contained in comparison_output
directory and <assembly_label>_output
directories for each assembly separately.
reads.x%-covered.list
– IDs of the isoforms from the database that have at least x% bases covered by all reads, where x is specified with --lower_threshold / --upper_threshold
options (50% / 95% by default).The following plots are similarly contained in both comparison_output
directory and <assembly_label>_output
directories. Please note, that most of the plots represent cumulative distributions and some plots are given in logarithmic scale.
Basic
transcript_length.png
_ – assembled transcripts length distribution (+ database isoforms length distribution).block_length.png
– alignment blocks length distribution (+ database exons length distribution).x-aligned.png
– transcript aligned fraction distribution.blocks_per_alignment.png
– distribution of number of blocks per alignment (+ distribution of number of database exons per isoform).alignment_multiplicity.png
– distribution for the number of significant alignment for each multiply-aligned transcript.mismatch_rate.png
_ – substitution errors per alignment distribution.Nx.png
– Nx plot for transcripts. Nx is a maximal number N, such that the total length of all transcripts longer than N bp is at least x% of the total length of all transcripts.NAx.png
– Nx plot for alignments.Sensitivity
x-assembled.png
– a histogram in which each bar represents the number of isoforms from the database that have at least x% captured by a single assembled transcript.x-covered.png
– a histogram in which each bar represents the number of isoforms from the database that have at least x% of bases covered by all alignments.x-assembled_exons.png
– a histogram in which each bar represents the number of exons from the database that have at least x% captured by a single assembled transcript.x-covered_exons.png
– a histogram in which each bar represents the number of exons from the database that have at least x% of bases covered by all alignments.alignments_per_isoform.png
– plot showing number of transcript alignments per isoformSpecificity
x-matched.png
– a histogram in which each bar represents the number of transcripts that have at least x% matched to an isoform from the database.x-matched_blocks.png
– a histogram in which each bar represents the number of all blocks from all transcript alignments that have at least x% matched to an isoform from the database.Your comments, bug reports, and suggestions are very welcomed. They will help us to further improve rnaQUAST. If you have any troubles running rnaQUAST, please send us logs/rnaQUAST.log
from the output directory.
Submit your issues and comments to our GitHub repository.