Download FASTQ files (optional)

SRA is a repository of biological sequences that stores data from many published articles with the potential to answer pressing biological questions.

The latest versions of psichomics support automatic downloading of SRA data from recount, a resource of pre-processed data for thousands of SRA projects (including gene read counts, splice junction quantification and sample metadata). First, check if the project of your interest is available in recount, thus making it quicker to analyse gene expression and alternative splicing for your samples of interest.

Data from SRA can be downloaded using the fasterq-dump command from sra-tools. For instance, to retrieve samples from the SRP126561 project:

# List SRA samples
samples=(SRR6368612 SRR6368613 SRR6368614 SRR6368615 SRR6368616 SRR6368617)

# Download samples
fasterq-dump --split-3 ${samples}

--split-3 allows to output one or two FASTQ files for single-end or paired-end sequencing, respectively (a third FASTQ file may also be returned containing orphaned single-end reads obtained from paired-end sequencing data)

Sample-associated data is also available from the Run Selector page. Click RunInfo Table to download the whole metadata table for all samples (usually downloaded in a file named SraRunTable.txt).

Align RNA-seq data to quantify splice junctions

The quantification of each alternative splicing event is based on the proportion of junction reads that support the inclusion isoform, known as percent spliced-in or PSI (Wang et al., 2008).

To estimate this value for each splicing event, both alternative splicing annotation and quantification of RNA-Seq reads aligning to splice junctions (junction quantification) are required. While alternative splicing annotation is provided by the package, junction quantification will need to be prepared from user-provided data by aligning the RNA-seq reads from FASTQ files to a genome of reference. As junction reads are required to quantify alternative splicing, a splice-aware aligner will be used.

psichomics currently supports STAR output.

mkdir hg19_STAR
STAR --runMode genomeGenerate \
     --genomeDir hg19_STAR \
     --genomeFastaFiles /path/to/hg19.fa \
     --sjdbGTFfile /path/to/hg19.gtf \
     --runThreadN 4

# Arguments:
#     --runMode             Generate the genome index
#     --genomeDir           Path to genome index (output)
#     --genomeFastaFiles    Path to genome FASTA file(s)
#     --sjdbGTFfile         Path to junction GTF annotation
#     --runThreadN 4        Run in parallel using 4 threads

align () {
    echo "Aligning ${1} using STAR..."
    STAR --readFilesIn ${1}_1.fastq ${1}_2.fastq \
         --runThreadN 16 \
         --genomeDir hg19_STAR \
         --readFilesCommand zcat \
         --quantMode GeneCounts \
         --outFileNamePrefix ${1}
}
# Arguments:
#     --readFilesIn              FASTQ files to align
#     --runThreadN 16            Run in parallel using 16 threads
#     --genomeDir                Path to genome index (input)
#     --readFilesCommand zcat    Use zcat to extract compressed files
#     --quantMode                Return gene read counts
#     --outFileNamePrefix        Prefix for output files
         
for each in ${samples}; do
    align "${each}"
done

Prepare output for psichomics

To process the resulting data files, type in R:

# Change working directory to where the STAR output is
setwd("/path/to/aligned/output/")

library(psichomics)
prepareGeneQuant(
    "SRR6368612ReadsPerGene.out.tab", "SRR6368613ReadsPerGene.out.tab",
    "SRR6368614ReadsPerGene.out.tab", "SRR6368615ReadsPerGene.out.tab",
    "SRR6368616ReadsPerGene.out.tab", "SRR6368617ReadsPerGene.out.tab")
prepareJunctionQuant("SRR6368612SJ.out.tab", "SRR6368613SJ.out.tab", 
                     "SRR6368614SJ.out.tab", "SRR6368615SJ.out.tab",
                     "SRR6368616SJ.out.tab", "SRR6368617SJ.out.tab")

To load the data, move the files (including the SRA metadata) to a new folder and follow the instructions in Load user-provided data into psichomics.

Sample information

If you are working with sample metadata from SRA Run Selector, see how to prepare SRA Run Selector data.

• Tab-separated values (TSV)
• Sample identifiers (rows) and their attributes (columns)
• The first column must contain sample identifiers and be named Sample ID
• Optionally, indicate the subject associated to each sample in a column named Subject ID (subject identifiers must be the same as the ones used in subject information)
• Example of a valid sample information dataset:

Subject information

• Tab-separated values (TSV)
• Subject identifiers (rows) and their attributes (columns)
• The first column must contain subject identifiers and be named Subject ID
• Example of a valid subject information dataset:

Gene expression

• Tab-separated values (TSV)
• Read counts of genes (rows) across sample (columns) (sample identifiers must be the same as the ones used in sample information)
• The first column must contain unique gene names (symbols, Ensembl ID, etc.) and be named Gene ID
• Example of a valid gene expression dataset:

Exon-exon junction quantification

• Tab-separated values (TSV)
• Read counts of exon-exon junctions (rows) across samples (columns) (sample identifiers must be the same as the ones used in sample information)
• The first column must contain junction identifiers and be named Junction ID
• Only chromosome number and capital letters X, Y, Z, W, and M are supported, followed by the genomic regions; acceptable junction identifiers include:
  • 10_18748_21822
  • chromosome 10 (18748 to 21822)
  • chr10:18748-21822
• Optionally, indicate the strand with + or - at the end of the junction identifier:
  • 10:3213:9402:+
  • chr10:3213-9402 -
• Junction identifiers whose chromosomes are alt, random or Un (i.e. alternative sequences) are discarded
• Example of a valid exon-exon junction quantification dataset:

Alternative splicing quantification (also known as inclusion levels)

Note that psichomics cannot currently parse alternative splicing events (e.g. identify the cognate gene and coordinates) from generic, user-provided tables.

• Tab-separated values (TSV)
• Quantification values of alternative splicing events (rows) across samples (columns) (sample identifiers must be the same as the ones used in sample information)
• The first column must be named AS Event ID
• Values between 0 and 1 or between 0 and 100: if the latter, values are automatically scaled between 0 and 1
• Example of a valid alternative splicing quantification dataset:

To load the data, move the files to a new folder and follow the instructions in Load user-provided data into psichomics.