SLRanger

An integrated approach for spliced leader detection and operon prediction in eukaryotes using long RNA reads

Workflow

Workflow of SLRanger. (A) The direct RNA sequencing workflow for species exhibiting trans-splicing. In such species, spliced leader (SL) sequences are located at the 5′ ends of RNA reads due to the trans-splicing mechanism. Mature mRNAs containing SL sequences are sequenced using Oxford Nanopore direct RNA sequencing, and the resulting long reads are subsequently aligned to the reference genome. (B) The SL detection module of SLRanger. After aligning long RNA reads to the reference genome, the unaligned 5′-end fragments are extracted and aligned against a reference set of SL (including SL1 sequence and all SL2 variant sequences) and random sequences (as controls). Using SLRanger’s scoring scheme, each candidate read is assigned an “SLRanger score.” Reads will be assigned to the corresponding SL-type according to the highest “SLRanger score”. Reads that can’t be classified as SL1 will be considered as the SL2 type. If SL1 and SL2 can’t be distinguished, the read is labeled SL_unknown. If variants of SL2 can’t be distinguished, the read is labeled SL2_unknown. A dynamic cutoff is then applied to identify “high-confidence SL reads”. (C) The principle underlying operon prediction by SLRanger. Based on the presence of high-confidence SL sequences in each read and their genomic mapping positions relative to gene annotations, operon structures are inferred. Genes with a high proportion of SL1-type reads are predicted to be upstream operon genes, whereas genes with a high proportion of SL2-type reads or supported by multiple SL2-type reads are predicted to be downstream operon genes.

Installation

The pipeline is invoked using a CLI written in Python(3.9-3.11) and requires a Unix-based operating system. For conda method, we provided installation from pypi and git

i. Conda method

1. Prepare a new conda env

conda create -n SLRanger_env python=3.9
conda activate SLRanger_env
conda install -c bioconda bedtools minimap2 samtools

2a. Install from PyPI

# install from pypi
pip install SLRanger

2b. Install from Github

#  install from git
git clone https://github.com/lrslab/SLRanger.git
cd SLRanger/
python -m pip install .

ii. Docker method

docker pull zhihaguo/slranger_env

Manual

SLRanger encompasses two primary functions, spliced leader (SL) detection and operon prediction, used to determine whether long RNA reads carry SL sequences and predict the operon structure based on the SL information.

1. Preprocessing

Reference and annotation selection

The long RNA reads will be mapped to the genome reference. The genome reference (fasta/fa/fna file) and annotation file (GFF file) should be determined before running SLRanger. These can be downloaded from NCBI or assembled independently.

In our sample folder,we provided C. elegans annotation file.

Long reads alignment

Additionally, we require users to provide their own alignment file (BAM file). For long reads, minimap2 is the recommended software. In the sample folder, we have provided test.bam, which was generated using the following command.

minimap2 -ax splice -uf -t 80 -k14 --MD --secondary=no $reference $basecall_file > tmp.sam
samtools view -hbS tmp.sam | samtools sort -@ 32 -F 260 -o test.bam
samtools index test.bam

2. Spliced Leader detection

SL_detect.py is designed to detect spliced leaders.

Command options

Available options can be viewed by running SL_detect.py -h in the command line.

SL_detect.py -h
usage: SL_detect.py [-h] -r REF -b BAM [-o OUTPUT] [--visualization] [-t CPU]

help to know spliced leader and distinguish SL1 and SL2

options:
  -h, --help            show this help message and exit
  -r REF, --ref REF     SL reference (fasta file recording SL sequence, required)
  -i BAM, --input BAM     input the bam file (required)
  -m , --mode           RNA or cDNA
  -o OUTPUT, --output OUTPUT
                        output file (default: SLRanger.txt)
  -t CPU, --cpu CPU     CPU number (default: 4)
  -c CUTOFF, --cutoff CUTOFF
                        The value used to filter high confident SL reads. 
                        The higher the value, the stricter it is. 
                        The range is between 0-10. (default: 4)                      

Output description

i. result table

ii. visualization result

The summary table and figures, including the Data Summary Table and the pictures including Cumulative Counts (SW), Cumulative Counts (SL), Query Length Distribution, Aligned Length Distribution, SL Type Distribution. will be output in a webpage format. An example is provided here.

Example

We provided test data to run as below.

git clone https://github.com/lrslab/SLRanger.git
cd sample/
unzip data.zip
SL_detect.py --ref SL_list_cel.fa --input RNA_test.bam -o SLRanger.txt -t 4 
SL_detect.py --ref SL_list_cel.fa --input cDNA_test.bam -o SLRanger_cDNA.txt -t 4

3. Operon prediction

operon_predict.py predicts operons and always writes a per-gene SL1/SL2 count table. Genes without high-confidence SL reads are omitted from that table. A read mapped to multiple semicolon-delimited genes contributes once to each gene.

Operon prediction can use either the original BAM input (-b/--bam) or an existing two-column read-to-gene mapping file (-m/--mapping). If only SL1 reads are detected, the count table is still written and operon prediction is skipped. SL2-only input is allowed to continue through operon prediction.

Command options

Available options can be viewed by running operon_predict.py -h in the command line.

operon_predict.py  -h
usage: operon_predict.py [-h] -g GFF (-b BAM | -m MAPPING) -i INPUT
                         [-o OUTPUT] [--gene-sl-table GENE_SL_TABLE]
                         [--sl1-map SL1_MAP] [--sl2-map SL2_MAP]
                         [-d DISTANCE] [-c CUTOFF]
help to know spliced leader and distinguish SL1 and SL2

options:
  -h, --help            show this help message and exit
  -g GFF, --gff GFF, -r GFF, --refer GFF
                        GFF annotation file
  -b BAM, --bam BAM     BAM file
  -m MAPPING, --mapping MAPPING
                        existing read-to-gene mapping file
  -i INPUT, --input INPUT
                        input the SL detection file
  -o OUTPUT, --output OUTPUT
                        output operon detection file (default: SLRanger.gff)
  --gene-sl-table GENE_SL_TABLE
                        per-gene SL1/SL2 count table
  --sl1-map SL1_MAP     comma-separated SL_type values treated as SL1
  --sl2-map SL2_MAP     comma-separated SL_type values treated as SL2
  -d DISTANCE, --distance DISTANCE
                        promoter scope (default: 5000)
  -c CUTOFF, --cutoff CUTOFF
                        cutoff of high-confidence SL reads (default: 4)

Output description

When operon prediction runs, a GFF file is returned. The per-gene count table is always written. By default it is placed next to the GFF and named <output_stem>_gene_sl1_sl2.tsv. It contains SL1, SL2, sum_count, sl2_ratio, and the gene classification columns used for prediction.

Without --sl1-map or --sl2-map, the legacy SLRanger classification is preserved: exact SL1 is treated as SL1, random, SL_unknown, and SL1_unknown are excluded, and the remaining recognized types are treated as SL2. Supply both map options to classify custom reference names explicitly.

Example

We provided test data to run as below (should be run after SL_detect.py).

cd sample/
operon_predict.py -g cel_wormbase.gff -b RNA_test.bam -i SLRanger.txt  -o test.gff

To reuse an existing read-to-gene mapping and custom SL names:

operon_predict.py -g annotation.gff -m reads_to_genes.tsv -i SLRanger.txt \
  -o operons.gff --sl1-map SL1 --sl2-map SL2,SL3,SL4

Cite our work

Our paper is online now. Please cite our work -- SLRanger: an integrated approach for spliced leader detection and operon prediction using long RNA reads on Briefings in Bioinformatics.