Skip to contents

Count Splice Junction Coverage from BAM Files

coverageSJ.Rd

Extracts splice junction information from RNA-seq alignments for all genes in the experiment. For each junction, reports the donor position, acceptor position, and the number of reads supporting the junction (splicedIn count).

Usage

coverageSJ(data, nbThreads = 1)

Arguments

data

List. Experiment object from newExperiment containing:

orthogroups

GRanges with metadata columns geneID, orthogroupID, speciesID.

sampleInfo

data.frame with columns sampleID, bam, isPairedEnd, strandSpecific, speciesID.

nbThreads

Integer. Number of genes processed in parallel. Default: 1.

Value

A data.frame with one row per junction per sample:

geneID

Gene identifier.

sampleID

Sample identifier.

donor

Genomic position of 3' end of the donor exon.

acceptor

Genomic position of the 5' end of the acceptor exon.

splicedIn

Number of reads spanning this junction.

Details

The function performs the following steps for each gene:

  1. Loads aligned reads from BAM files for samples of the matching species.

  2. Adjusts strand based on library type (strandSpecific parameter).

  3. Extracts splice junctions from CIGAR strings (N operations).

  4. Filters junctions to those strictly within the gene boundaries.

  5. Returns donor/acceptor positions and read counts.

Strand-Specific Libraries

The strandSpecific column in sampleInfo controls strand handling:

  • 0: Unstranded library. Junctions from both strands are kept.

  • 1: Stranded. Read 1 matches the transcript strand (e.g., ligation-based protocols; TopHat fr-secondstrand).

  • 2: Reversely stranded. Read 1 is antisense to the transcript (e.g., dUTP, Illumina TruSeq Stranded; TopHat fr-firststrand). Read strand is inverted.

Parallelization

Genes are processed in parallel using nbThreads workers.

See also

newExperiment to create the experiment object, findOSJ to aggregate junctions across species (next step).

Examples

if (FALSE) { # \dontrun{
data <- newExperiment(sampleInfo, speciesInfo, orthogroups,
                      dnaDir, msaDir, speciesTree)

# Quantify junctions (8 threads)
coverage <- coverageSJ(data, nbThreads = 8)

head(coverage)
#   geneID    sampleID  donor  acceptor  splicedIn
# 1 ENSG001   sample_1  999    1501      42
# 2 ENSG001   sample_1  999    2001      15
# Note: donor/acceptor are exon boundary coordinates (not intron coordinates)
} # }