Annotation file / assembled transcripts evaluation (GFF)#

You have RNA-Seq reads and are running an RNA-Seq analysis pipeline. You have performed quality control, aligned the reads to the genome, and assembled the alignments into transcripts. Now, you have obtained a GFF or GTF file generated by a transcriptome assembler. The next question you want to address is the quality of these assembled transcripts.

A common approach in the pipeline involves using gffcompare to compare the assembled transcripts with a reference annotation file, such as RefSeq, Gencode, or CHESS. However, this method assumes the correctness of the annotation file and evaluates how closely our assembled transcripts align with the annotations.

It is important to acknowledge that annotation files can contain misannotations and noise. Even the most meticulously curated human annotations are not exempt from errors. Moreover, there are even more errors in other species' annotation files.

Important

We provide a new approach to evaluate transcripts by counting how many bad splice junctions are in a transcript!

You can run Splam on (1) annotation files or (2) assembled transcripts. Splam outputs the score of every donor and acceptor site. By using these scores,

  1. You can get an overview of the quality of all introns (splice junctions) inside your annotation file.

  2. You can assess each transcript by checking how many bad splice junctions a transcript has, which can be a good filtering criteria for assembled transcripts.

  3. There are more potential usages to be explored!


Step 1: Preparing your input files#

The first step is to prepare three files for Splam analysis. In this tutorial, we will be using a toy dataset.

Input files

  1. An annotation file in GFF or GTF format [example file: refseq_110_GRCh38_chr_fixed.gff].

  2. A reference genome in FASTA format [example file: chr9_subset.fa].

  3. The Splam model, which you can find it here: splam.pt


Step 2: Extracting introns in your annotation file#

In this step, you take an annotation file (1) and run

splam extract refseq_40_GRCh38.p14_chr_fixed.gff -o tmp_out_annotation

Splam uses gffutils to extract introns from all transcripts. It starts by creating an sqlite3 database and then iterates through all exons in the GFF file and writes intron coordinates into a BED file.

By default, the BED is written into tmp_out/junction.bed. The BED file consists of six columns: CHROM, START, END, JUNC_NAME, INTRON_NUM, and STRAND. Here are a few entries from the BED file:

Output: junction.bed

1chr9    4849549 4860125 JUNC00000007    3       +
2chr9    5923308 5924658 JUNC00000008    6       -
3chr9    5924844 5929044 JUNC00000009    8       -

Here are some optional arguments:

-o / --outdir DIR

The directory where the output file is written to. The default output directory is tmp_out. You can set your own output directory using this argument.

-f / --file-format FILE_FORMAT

Splam automatically detects whether your input file is a BAM or GFF file based on its extension. In this section, we are using Splam to evaluate a given annotation file, so please ensure that your input file has a .gff, .gtf, .GFF, or .GTF extension.

-d / --database DATABASE

The path to the annotation database built using gffutils. If this argument is provided, Splam loads the database instead of creating a new one.


Step 3: Scoring extracted introns#

In this step, the goal is to score all the extracted splice junctions. To accomplish this, you will need 3 essential files. (1) The BED file that was generated in Step 2, (2) the reference genome (2) which shares coordinates with the junction BED file, and (3) the Splam model (3). Once you have these files in place, you can run the following command:

splam score -G chr9_subset.fa -m ../model/splam_script.pt -o tmp_out_annotation tmp_out_annotation/junction.bed

In this step, a new BED file is produced, featuring eight columns. Two extra columns, namely DONOR_SCORE and ACCEPTOR_SCORE, are appended to the file. It is worth noting that any unstranded introns are excluded from the output. (P.S. They might be from unstranded transcripts assembled by StringTie).

Output: junction_score.bed

1chr9    4849549 4860125 JUNC00000007    3       +       0.7723698       0.5370769
2chr9    5923308 5924658 JUNC00000008    6       -       0.9999831       0.9999958
3chr9    5924844 5929044 JUNC00000009    8       -       0.9999883       0.9999949

Here are the required arguments:

-G / --reference-genome REF.fasta

The path to the reference genome in FASTA format. Please ensure that this file shares the same coordinates as your input alignment file, which is where you align your RNA-Seq reads. Splam will handle the indexing process for you if the reference genome has not been indexed yet.

-m / --model MODEL.pt

This argument is the path to the trained Splam model. If you haven't downloaded the Splam model yet, here is the link.

Here are some optional arguments:

-A / --assembly-report REPORT

The path to an assembly report file in tsv format which contains the chromosome identifiers and lengths. This information is built into Splam if running on a human genome (defaults to human GRCh38, patch 14). However, this argument is required if running on non-human species. See our mouse example for reference.

-d / --device pytorch_DEV

By default, Splam automatically detects your environment and runs in cuda mode if CUDA is available. However, if your computer is running macOS, Splam will check if mps mode is available. If neither cuda nor mps are available, Splam will run in cpu mode. You can explicitly specify the mode using the -d / --device argument.

-b / --batch-size BATCH

Additionally, you can adjust the batch size using the -b / --batch-size argument. This argument defines the number of samples that will be propagated through the Splam network. By default, the batch size is set to 10. We recommend setting a small batch size (for instance 2) when running Splam in cpu mode.

-o / --outdir DIR

The directory where the output file is written to. The default output directory is tmp_out. This argument is same as the one in Step 2. Note that if you set your own output directory, you have to set the same output directory for this step as well. Otherwise, Splam will not be able to find some essential temporary files. We recommend users not to set this argument and use the default value.


Step 4: Evaluating isoforms by Splam scores#

To summarize the quality of each isoform, users can run Splam to count how many spurious splice junctions are in each transcript.

splam clean -o tmp_out_annotation -t 0.8

Output: cleaned.gff

The output file of this step is a sorted Splam-cleaned GFF file. You can replace the original GFF file with the cleaned GFF file. Also, reference the report.tsv for a list of transcripts, and their ratio and number of spurious/removed junctions.

Splam score threshold suggestion

For evaluating the accuracy of GFF annotation files, we advise using a stricter score threshold of 0.8.

Here are some optional arguments:

-t / --threshold threshold

This is the score cutoff threshold for Splam to determine whether a given splice junction is spurious (discarded) or not. It is a floating-point value between 0 and 1. If the score of either the donor or acceptor site falls below this value, then any spliced alignments containing this junction will be removed. The default threshold is set to 0.1.

-o / --outdir DIR

The directory where the output file is written to. The default output directory is tmp_out. This argument is same as the one in Step 2 and Step 3. Note that if you set your own output directory, you have to set the same output directory for this step as well, or otherwise, Splam will not be able to find some essential temporary files. We recommend users not to set this argument and use the default value.


What's next?#

Congratulations! You have finished this tutorial.

See also