Skip to content

Tutorial

Jump to bottom
Tianyuan Liu edited this page Apr 28, 2024 · 3 revisions

SQANTI-SIM Tutorial

Tutorial 1: Running SQANTI-SIM on an example dataset

In this tutorial, we will walk you through the process of using SQANTI-SIM with the example toy dataset provided. SQANTI-SIM allows you to simulate long-read RNA-Seq data and evaluate the performance of your transcript reconstruction methods. You can execute SQANTI-SIM with a single command, full-sim, which wraps all three main simulation steps (classif, design and sim).

Prerequisites

Before you begin, make sure you have SQANTI-SIM installed. If not, please refer to the Requirements and Installation guide.

First, create a working directory and activate the SQANTI-SIM conda environment:

mkdir sqanti-sim_tutorial_1
cd sqanti-sim_tutorial_1
conda activate SQANTI-SIM.env

After activating the Anaconda environment, your prompt should look like this:

(SQANTI-SIM.env)$

Keep the environment active while using the SQANTI-SIM scripts. To complete the setup, add SQANTI-SIM to the $PYTHONPATH and $PATH:

export PATH="$PATH:<path_to>/SQANTI-SIM/"

Simulate ONT data

We will use the example datasets located in <path_to>/SQANTI-SIM/example/mouse (mm39_refGene_subset.gtf and mm39_chr1_subsequence.fa to simulate long reads from 1000 transcripts. Of these, 150 are novel transcripts (50 ISM, 50 NIC, and 50 NNC). This simulation will generate ONT (--ont) cDNA reads (you can change this using --read_type dRNA) from 1000 different isoforms (--trans_number 1000), with 50 ISM, 50 NIC, and 50 NNC transcripts. Since we are running the simulation in sample mode and not specifying the number of long reads to simulate (not using --long_count), expression values will be assigned based on empirical pre-analyzed expression distributions.

sqanti-sim.py full-sim sample \
    --gtf <path_to>/SQANTI-SIM/example/mouse/mm39_refGene_subset.gtf \
    --genome <path_to>/SQANTI-SIM/example/mouse/mm39_chr1_subsequence.fa \
    --trans_number 1000 \
    --ISM 50 --NIC 50 --NNC 50 \
    --ont \
    --seed 1234

This will generate an ONT cDNA reads file (ONT_simulated.fastq) and a reduced GTF annotation file (sqanti-sim_modified.gtf) to simulate novelty. These files are essential to incorporate controlled novelty into your dataset.

In this example, we provide a reconstructed transcriptome using this simulated data (example_detected_isoforms.gtf). To run the SQANTI-SIM evaluation, you need to provide your reconstructed transcriptome, the reference genome, the reduced reference annotation, and the SQANTI-SIM index file, which contains structural annotation and metadata of this simulation process.

sqanti-sim.py eval \
    --transcriptome <path_to>/SQANTI-SIM/example/mouse/example_detected_isoforms_mm39.gtf \
    --gtf sqanti-sim_modified.gtf \
    --genome <path_to>/SQANTI-SIM/example/mouse/mm39_chr1_subsequence.fa \
    -i sqanti-sim_index.tsv \
    -d tutorial_eval

Once finished it will generate a report with the accuracy of the reconstructed transcriptome in recovering the simulated dataset.

Congratulations! 👏 You've completed your first SQANTI-SIM analysis. To explore a more advanced and controlled simulation process, check out the next tutorial. For a complete understanding of all features, consult the SQANTI-SIM wiki.

Tutorial 2: Running SQANTI-SIM step by step to simulate multi-strategy dataset and evaluate reconstructed transcriptome

This tutorial guides you through the complete SQANTI-SIM pipeline, covering the following topics:

In this tutorial, we will evaluate the performance of a transcript identification pipeline using SQANTI-SIM. We'll start by classifying reference transcripts from an example dataset into SQANTI3 structural categories. Then, we'll decide how many novel and known transcripts to simulate and proceed with the simulation. After that, we'll run a transcript identification pipeline using the modified reference annotation and simulated reads. Finally, we'll evaluate its performance and generate the SQANTI-SIM report.

Before you begin, make sure you have SQANTI-SIM installed, the Anaconda environment is activated and cDNA_Cupcake is available.

mkdir sqanti-sim_tutorial_2
cd sqanti-sim_tutorial_2
conda activate SQANTI-SIM.env

export PYTHONPATH=$PYTHONPATH:<path_to>/cDNA_Cupcake/sequence/
export PYTHONPATH=$PYTHONPATH:<path_to>/cDNA_Cupcake/
export PATH="$PATH:<path_to>/SQANTI-SIM/"

Step 1. Classify reference transcripts

To simulate transcripts based on SQANTI3 structural category definitions, we need to classify all the reference transcripts into their potential structural categories, assuming that the query transcript cannot match with itself. We'll use the classif mode of SQANTI-SIM, which implements an adapted version of the SQANTI3 structural category classification algorithm.

Please note that the reference annotation must be supplied in GTF format. An example of a reference transcriptome and its required format is available from GENCODE. In this tutorial, we'll use a subset of the Homo sapiens Ensembl genome assembly and gene annotation. We'll start by downloading an example dataset for human chromosome 12:

wget http://ftp.ensembl.org/pub/release-107/fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna.chromosome.12.fa.gz
gunzip Homo_sapiens.GRCh38.dna.chromosome.12.fa.gz 

wget http://ftp.ensembl.org/pub/release-107/gtf/homo_sapiens/Homo_sapiens.GRCh38.107.chr.gtf.gz
gunzip Homo_sapiens.GRCh38.107.chr.gtf.gz
grep "^12" Homo_sapiens.GRCh38.107.chr.gtf > Homo_sapiens.chr12.gtf
rm Homo_sapiens.GRCh38.107.chr.gtf

To classify the reference transcripts into potential SQANTI3 structural categories, run SQANTI-SIM in classif mode:

sqanti-sim.py classif \
    --gtf Homo_sapiens.chr12.gtf

This will generate an index file called prefix_index.tsv with all the transcripts of the reference annotation classified into their potential structural categories. In the terminal a summary table with the number of transcripts classified in each structural category will be printed.

_______________________________________________________________________________
S Q A N T I S I M 📊

Classification summary Table 🔎
_______________________________________________________________________________
| full-splice_match: 255
| incomplete-splice_match: 2566
| novel_in_catalog: 4385
| novel_not_in_catalog: 4552
| fusion: 121
| antisense: 407
| genic_intron: 3
| genic: 92
| intergenic: 1290

Step 2. Design experiment for simulation

Before simulating novel reads, we need to decide how many different transcripts from different structural categories we want to simulate and at what expression levels. This will generate a modified reference annotation without the simulated novel transcripts and will add three new columns to the index file: (sim_type) if the read is in the reference annotation known or if we don't have knowledge about it novel, (requested_counts) the number of reads to be simulated for that transcript, and (requested_tpm) the TPM requested for that transcript.

SQANTI-SIM implements 3 different modes to assign this expression values to the transcripts:

  • equal: Assigns the same expression value to every transcript. Expression levels won't have any variance but you will have a very controlled scenario in which every transcript will be simulated with the same coverage.
  • custom: Allows you to specify different expression values for known and novel transcripts using negative binomial distributions.
  • sample: Generates expression distributions based on real sample data. In this mode you are able to simulate the same or different expression between novel and known transcripts, and also to reproduce the isoform complexity profile (number of different expressed isoforms per gene).

In this example, we'll simulate in sample mode using default pre-analyzed expression distributions from PacBio data. A total of 1800 known transcripts and 1200 novel transcripts from the three main structural categories of SQANTI3 will be simulated: 400 ISM, 400 NIC and 400 NNC.

sqanti-sim.py design sample \
    -i sqanti-sim_index.tsv \
    --gtf Homo_sapiens.chr12.gtf \
    --genome Homo_sapiens.GRCh38.dna.chromosome.12.fa \
    --pb \
    --trans_number 3000 \
    --ISM 400 --NIC 400 --NNC 400 \
    --seed 1234

This generates the modified reference annotation (prefix_modified.gtf) and adds new columns to the index file for expression levels.

Step 3. Simulate reads

Now, choose the sequencing platform for simulating your long reads. We'll simulate PacBio, Illumina and CAGE data in this tutorial. To keep the simulation manageable, we'll create a small sample:

sqanti-sim.py sim \
    -i sqanti-sim_index.tsv \
    --gtf Homo_sapiens.chr12.gtf \
    --genome Homo_sapiens.GRCh38.dna.chromosome.12.fa \
    --pb --pass_num 1 \
    --illumina --CAGE \
    --short_count 10000 \
    --seed 1234

Step 4. Run a transcriptome reconstruction pipeline

This step is independent of the SQANTI-SIM pipeline. You can run your preferred transcriptome reconstruction pipeline and provide the modified reference annotation and simulated reads to it. In this example, we'll use FLAIR to demonstrate the process. IMPORTANT: Remember to provide the modified reference annotation and the simulated reads to your transcriptome reconstruction pipeline.

In order to run our transcriptome reconstruction pipeline, we have to ensure that we are using the simulated reads as input and the modified reference annotation as reference. In this example, we'll use FLAIR to demonstrate the process.

Install FLAIR by referring to its GitHub. Then, run FLAIR using only long-read data and generate the long-read-defined transcriptome:

flair 123 \
    -r PBSIM3_simulated.fasta \
    -g Homo_sapiens.GRCh38.dna.chromosome.12.fa \
    -f sqanti-sim_modified.gtf \
    -o example_detected

Step 5. Evaluate identification and discovery performance

Finally, generate a report with the evaluation metrics of the performance of the reconstruction algorithm by running SQANTI-SIM in eval mode:

sqanti-sim.py eval \
    --transcriptome example_detected.isoforms.gtf \
    --gtf sqanti-sim_modified.gtf \
    --genome Homo_sapiens.GRCh38.dna.chromosome.12.fa \
    -i sqanti-sim_index.tsv \
    -d tutorial_eval \
    -o example \
    --CAGE_peak CAGE_simulated.bed \
    --coverage STAR_mapping/Illumina_simulated_1SJ.out.tab

This will produce a report detailing the performance of your transcriptome reconstruction pipeline in terms of transcript identification and discovery.

🏁By following these steps, you've learned how to use SQANTI-SIM to classify reference transcripts, design an experiment for simulation, simulate reads, and evaluate the performance of a transcriptome reconstruction pipeline. Feel free to explore more advanced features of SQANTI-SIM and adapt this tutorial to your specific analysis needs.

Clone this wiki locally