RBP workflow
library(SummarizedExperiment)
library(GenomicFeatures)
library(dplyr)
library(ggplot2)
library(satuRn)
library(limma)
library(tibble)
library(here)
devtools::load_all(here("../splicelogic"))# dir <- system.file("extdata", package = "splicelogic")
# se <- readRDS(file.path(dir, "GTEx_v9_ONT_SummarizedExperiment.rds"))
se <- readRDS(here("workflows/data/GTEx_v9_ONT_SummarizedExperiment.rds"))
dim(se)## [1] 385669 11keep <- rowSums(assay(se) >= 10) >= 6
se <- se[keep, ]
dim(se)## [1] 18985 11# filters: 385669 -> 18985 transcripts
library(DESeq2)
library(tximeta)
colsums <- colSums(assay(se, "counts"))
cpm <- t(t(assay(se, "counts")) / colsums * 1e6)
assay(se, "abundance") <- cpm
# dummy lengths... no correction for length
assay(se, "length", withDimnames=FALSE) <- matrix(1000, nrow=nrow(se), ncol=ncol(se))DTU testing with satuRn
We compare the experiment condition (exp, RBP knockdown)
against control (ctrl) using satuRn.
colData(se)$condition <- dplyr::case_when(
colData(se)$condition == "ctrl" ~ "WT",
colData(se)$condition == "exp" ~ "KD"
)
# KD as reference level
colData(se)$condition <- factor(colData(se)$condition, levels = c("KD", "WT"))
# build design from the factor directly (matches satuRn vignette convention)
condition <- colData(se)$condition
donor <- factor(colData(se)$donor)
design <- model.matrix(~ 0 + condition + donor)
L <- limma::makeContrasts(
WT_vs_KD = conditionWT - conditionKD,
levels = design
)
L## Contrasts
## Levels WT_vs_KD
## conditionKD -1
## conditionWT 1
## donorQV44 0
## donorS4Z8 0
## donorS95S 0rowData(se)$isoform_id <- rowData(se)$tx_name
# flatten CharacterList gene_id to plain character
rowData(se)$gene_id <- vapply(rowData(se)$gene_id, `[`, character(1), 1)
se <- se[
duplicated(rowData(se)$gene_id) |
duplicated(rowData(se)$gene_id, fromLast = TRUE),
]
# 18985 transcripts -> 11874 transcripts after removing single-isoform genes
se <- satuRn::fitDTU(se, formula = ~0 + condition + donor)se <- satuRn::testDTU(
object = se,
contrasts = L,
diagplot1 = TRUE,
diagplot2 = TRUE
)
gene_ids <- rowData(se)$gene_id
# every tested transcript; `sig_res` below is the significant subset, and the
# non-significant rows are the background pool for slSampleNonEventExons()
dtu_res <- rowData(se)[["fitDTUResult_WT_vs_KD"]] |>
as.data.frame() |>
tibble::rownames_to_column("tx_name") |>
tibble::as_tibble() |>
dplyr::mutate(gene_id = gene_ids) |>
dplyr::select(tx_name, gene_id, estimates, pval, regular_FDR) |>
dplyr::arrange(pval)
sig_res <- dtu_res |> dplyr::filter(regular_FDR < 0.1)
sig_res <- sig_res |> mutate( tx_id = tx_name)
sig_res## # A tibble: 1,987 × 6
## tx_name gene_id estimates pval regular_FDR tx_id
## <chr> <chr> <dbl> <dbl> <dbl> <chr>
## 1 ENST00000378292.9 ENSG00000198467.16 2.98 1.00e-14 1.19e-10 ENST000…
## 2 ENST00000355572.6 ENSG00000196923.14 -3.01 5.09e-12 3.02e- 8 ENST000…
## 3 ENST00000355841.7 ENSG00000196923.14 3.63 7.77e-12 3.08e- 8 ENST000…
## 4 ENST00000403994.9 ENSG00000140416.26 -3.92 1.66e-11 4.94e- 8 ENST000…
## 5 ENST00000311785.11 ENSG00000138674.17 1.84 2.18e-11 5.17e- 8 ENST000…
## 6 ENST00000358278.7 ENSG00000140416.26 2.82 3.85e-11 7.63e- 8 ENST000…
## 7 ENST00000292807.9 ENSG00000161203.14 -2.68 6.14e-11 8.47e- 8 ENST000…
## 8 ENST00000382456.7 ENSG00000161203.14 2.68 6.14e-11 8.47e- 8 ENST000…
## 9 ENST00000613865.5 ENSG00000138326.21 1.53 6.42e-11 8.47e- 8 ENST000…
## 10 ENST00000645482.3 ENSG00000198467.16 -3.50 7.22e-11 8.58e- 8 ENST000…
## # ℹ 1,977 more rows# save sig_res for later use
saveRDS(sig_res, here("workflows/data/sig_res.rds"))# load the saved DTU results
sig_res <- readRDS(here("workflows/data/sig_res.rds"))
sig_res## # A tibble: 1,987 × 6
## tx_name gene_id estimates pval regular_FDR tx_id
## <chr> <chr> <dbl> <dbl> <dbl> <chr>
## 1 ENST00000378292.9 ENSG00000198467.16 2.98 1.00e-14 1.19e-10 ENST000…
## 2 ENST00000355572.6 ENSG00000196923.14 -3.01 5.09e-12 3.02e- 8 ENST000…
## 3 ENST00000355841.7 ENSG00000196923.14 3.63 7.77e-12 3.08e- 8 ENST000…
## 4 ENST00000403994.9 ENSG00000140416.26 -3.92 1.66e-11 4.94e- 8 ENST000…
## 5 ENST00000311785.11 ENSG00000138674.17 1.84 2.18e-11 5.17e- 8 ENST000…
## 6 ENST00000358278.7 ENSG00000140416.26 2.82 3.85e-11 7.63e- 8 ENST000…
## 7 ENST00000292807.9 ENSG00000161203.14 -2.68 6.14e-11 8.47e- 8 ENST000…
## 8 ENST00000382456.7 ENSG00000161203.14 2.68 6.14e-11 8.47e- 8 ENST000…
## 9 ENST00000613865.5 ENSG00000138326.21 1.53 6.42e-11 8.47e- 8 ENST000…
## 10 ENST00000645482.3 ENSG00000198467.16 -3.50 7.22e-11 8.58e- 8 ENST000…
## # ℹ 1,977 more rowslibrary(txdbmaker)
# load TxDb for GRanges conversion
# download once (warning: ~50MB compressed)
wget https://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_48/gencode.v48.annotation.gtf.gz
txdb <- txdbmaker::makeTxDbFromGFF(
"/Users/beacm/Documents/splicelogic-paper/workflows/gencode.v48.annotation.gtf.gz",
format = "gtf"
)
# save to not rebuild every time
AnnotationDbi::saveDb(txdb, here("workflows/data/gencode.v48.sqlite"))
# check how many transcripts in se are in txdb
rowData(se)$tx_name %in% AnnotationDbi::keys(txdb, "TXNAME") |> sum() / length(rowData(se)$tx_name)Splice event detection with splicelogic
# build exon GRangesList and rename from internal TXID to TXNAME
txdb <- AnnotationDbi::loadDb(here("workflows/data/gencode.v48.sqlite"))
exons <- exonsBy(txdb, by = "tx")
txps <- AnnotationDbi::select(txdb, keys(txdb, "TXID"), c("TXNAME", "GENEID"), "TXID") |>
tibble::as_tibble() |>
dplyr::mutate(TXID = as.character(TXID))
# rename exons to match tx_name in se, matching by TXID
all.equal(names(exons), txps$TXID)
names(exons) <- txps$TXNAME[match(names(exons), txps$TXID)]
saveRDS(exons, here("workflows/data/glinos_exons.rds"))# load from saved sqlite
txdb <- AnnotationDbi::loadDb(here("workflows/data/gencode.v48.sqlite"))
# load exons from saved rds
exons <- readRDS(here("workflows/data/glinos_exons.rds"))# prepare exons for splicelogic
sig_exons <- prepare_exons(
txdb, sig_res, coef_col = "estimates", verbose = TRUE
)## Extracting exons from TxDb...## Mapping transcript IDs...## Merging DTU results onto exons...## Done. Returned 22481 exon ranges from 1987 unique transcripts.#find se
se_events <- sig_exons |> preprocess(coef_col = "estimates")|> find_se()
#find all events
events <- sig_exons |> preprocess(coef_col = "estimates")|> find_all_events()## Calculating skipped exon events...## Calculating included exon events...## Calculating mutually exclusive exon events...## Calculating retained intron events...## Calculating spliced intron events...## Calculating alternative 5' splice site events...## Calculating alternative 3' splice site events...## Calculating alternative transcription start site events...## Calculating alternative transcription end site events...## Done! 2863 events detected.# define a helper function
my_paste <- \(x) paste(unique(x), collapse = ",")
# rows are collapsed when they share the same exon and the same event type. The
# exon is identified by `exon_key` (gene, chromosome, start, end and strand),
# built as in the vignette's recipe for an exon id, rather than by `exon_id`:
# the annotation gives the same exon more than one id when it is coding in one
# transcript and non-coding in another
collapsed_events <- events |>
dplyr::mutate(exon_key = paste(gene_id, seqnames, start, end, strand, sep = ":")) |>
dplyr::group_by(exon_key, event_type) |>
plyranges::reduce_ranges_directed(
n = plyranges::n(),
exon_id = my_paste(exon_id),
tx_id = my_paste(tx_id),
exon_rank = my_paste(exon_rank),
event_tx_id = my_paste(event_tx_id),
)event_counts <- events |>
as_tibble() |>
dplyr::group_by(event_type) |>
dplyr::summarise(count = dplyr::n())
collapsed_event_counts <- collapsed_events |>
as_tibble() |>
dplyr::group_by(event_type) |>
dplyr::summarise(count = dplyr::n())
# combine both count tables into one, labelled by source
counts <- dplyr::bind_rows(
"All events" = event_counts,
"Collapsed events" = collapsed_event_counts,
.id = "set"
)
# event types from most to fewest events
counts$event_type <- factor(
counts$event_type,
levels = event_counts$event_type[order(-event_counts$count)]
)
ggplot(counts, aes(x = event_type, y = count, fill = set)) +
geom_col(position = position_dodge(width = 0.8)) +
scale_fill_manual(values = c(
"All events" = "#B0CFFF", # light blue
"Collapsed events" = "#3A86FF" # darker blue
)) +
labs(title = "Splice event counts",
x = "Event Type", y = "Count", fill = NULL) +
theme_bw()
# use splicelogic's built-in plotting function to visualize events for one gene,
# picked at random among the genes with events (the seed keeps the same pick on
# every run, as long as the event set does not change)
set.seed(14)
gene_of_interest <- sample(unique(events$gene_id), 1)
# filter events for the gene of interest
events_for_gene <- events[events$gene_id == gene_of_interest]
# visualize the events for the gene of interest
splicelogic::plot_transcripts(events_for_gene, sig_exons, rescale_introns = TRUE) +
ggtitle(paste("Splice events for gene:", gene_of_interest))
# make a plot for each gene with events and save to files
# plot_transcripts_by_gene returns a list of ggplot objects, one for each gene
plots <- plot_transcripts_by_gene(events, sig_exons)
# save each plot to the same pdf file - takes a while: -- 1725 genes
pdf(here("workflows/figures/splice_events_by_gene.pdf"), width = 8, height = 6)
for (gene in names(plots)) {
print(plots[[gene]])
}
dev.off()Collapse redundant events
find_se() returns an exon once per transcript pair it is
part of: an exon in two down transcripts and skipped by two up ones
comes back four times. The event counts above already collapse these
(collapsed_events). For the sequence analyses below we also
keep one row per exon, so that an exon shared by many transcripts does
not weigh more in the profiles and logos, and the background set is
sized by exons rather than pairs. n records how many pairs
each exon came from. The collapse is the same as in
collapse_redundant_events, applied to the SE events, and
like it collapses rows only for the same exon (gene, chromosome, start,
end, strand) and the same event type.
# TRUE: the sequence analyses below use one row per exon. FALSE: they use every
# SE event, one row per transcript pair, as before the collapse, to compare
collapse_se <- TRUE
se_collapsed <- se_events |>
dplyr::mutate(exon_key = paste(gene_id, seqnames, start, end, strand, sep = ":")) |>
dplyr::group_by(exon_key, event_type) |>
plyranges::reduce_ranges_directed(
n = plyranges::n(),
gene_id = my_paste(gene_id),
exon_id = my_paste(exon_id),
tx_id = my_paste(tx_id),
exon_rank = my_paste(exon_rank),
event_tx_id = my_paste(event_tx_id),
)
data.frame(
set = c("SE"),
events = c(length(se_events)),
exons = c(length(se_collapsed)),
max_n = c(max(se_collapsed$n))
)## set events exons max_n
## 1 SE 303 227 12# the set the sequence analyses below use
se_exons <- if (collapse_se) se_collapsed else se_eventsSE exons per gene
How the SE events spread across genes: the number of distinct SE
exons in each gene. se_collapsed already has one row per
exon and gene, with n the events it came from, so the
counts come straight from it.
se_per_gene <- se_collapsed |>
dplyr::group_by(gene_id) |>
dplyr::summarise(n_events = sum(n),
n_exons = plyranges::n())
data.frame(
se_events = length(se_events),
se_exons = length(se_collapsed),
se_genes = nrow(se_per_gene)
)## se_events se_exons se_genes
## 1 303 227 190# ggplot needs a data.frame, so the small per-gene summary is converted here
as.data.frame(se_per_gene) |>
dplyr::count(n_exons) |>
ggplot(aes(x = factor(n_exons), y = n)) +
geom_col(fill = "#3A86FF", width = 0.6) +
geom_text(aes(label = n), vjust = -0.4, size = 3.5, colour = "grey20") +
scale_y_continuous(expand = expansion(mult = c(0, 0.1))) +
labs(title = "Distinct SE exons per gene",
subtitle = sprintf("%d SE events, %d distinct exons, %d genes",
length(se_events), length(se_collapsed),
nrow(se_per_gene)),
x = "Distinct SE exons in the gene", y = "Genes") +
theme_bw()
Sequence context around the event exons, focusing on skipped exons (SE)
The analysis below uses a set of sequence-context helpers
consolidated in helper.R (nucleotide-composition profiles,
background exon sampling, splice-site logos and RBP motif scanning). The
splicelogic package itself only exports get_seq() for
sequence work, so these stay as paper-repo helpers (for now).
source(here("workflows/helper.R"))slEventFlankProfiles() takes a finder’s output straight
away: it flanks each exon (strand-aware, so “upstream” is 5’ in
transcription order), slides overlapping windows across the flank, and
returns the per-window nucleotide fractions for each side.
se_profiles <- slEventFlankProfiles(se_exons, genome = "hg38", width = 100)
# rows = exons, cols = window_nucleotide (e.g. w1_A)
dim(se_profiles$upstream)## [1] 227 76dim(se_profiles$downstream)## [1] 227 76Background set
We want to compare the sequence context of the events with that of a
control set. For the control set (non-event exons), we sample as many
exons as there are in se_exons (one per SE exon, or one per
SE event with collapse_se = FALSE), drawn at random from
the internal exons of transcripts with no evidence of differential usage
(pval > 0.5).
non_event <- slSampleNonEventExons(
dtu_res, exons, n = se_exons, pval_min = 0.5, seed = 1
)non_event_profiles <- slEventFlankProfiles(
non_event, genome = "hg38", width = 100
)
dim(non_event_profiles$upstream)## [1] 227 76The SE profiles stacked over the non-SE ones, on the same axes:
patchwork::wrap_plots(
slPlotFlankProfiles(se_profiles, exon_label = "SE"),
slPlotFlankProfiles(non_event_profiles, exon_label = "non-SE"),
ncol = 1
)
The two sets, one panel per nucleotide (mean ± standard error), and next to them the difference between them, SE minus non-SE:
se_vs_non_event <- list(SE = se_profiles, "non-SE" = non_event_profiles)
patchwork::wrap_plots(
slPlotFlankProfilesByNt(se_vs_non_event),
slPlotFlankProfilesByNt(se_vs_non_event, difference = TRUE),
nrow = 1
)
An alternative, within-transcript control: for each event, the internal exon of the same transcript that is shared with the event partner and closest in rank to the event exon.
matched <- slMatchedNonEventExons(se_events, exons)
slEventFlankProfiles(matched, genome = "hg38", width = 100) |>
slPlotFlankProfiles(exon_label = "matched exon")Splice-site sequence logos
The window profiles above average over 10 bp, which smooths the
splice site away. slFlankLogoMatrices() keeps every base
position separate and returns a position frequency matrix per side,
running width bases into the intron and
exon_bases into the exon so the boundary is inside the
matrix.
Note which site each side carries: the region
upstream of an exon is the end of the preceding intron
and holds the 3’ splice site (acceptor,
...ag | EXON); the region downstream is
the start of the following intron and holds the 5’ splice
site (donor, EXON | gu...).
se_logo <- slFlankLogoMatrices(se_exons, genome = "hg38",
width = 25, exon_bases = 3)
# letter height = information content, the classic logo scaling: unconstrained
# positions shrink to nothing and the consensus stands out
slPlotFlankLogo(se_logo, method = "bits")
The same for the background set, to see whether the event splice sites differ from those of exons in transcripts with no differential usage.
non_event_logo <- slFlankLogoMatrices(non_event, genome = "hg38",
width = 25, exon_bases = 3)
slPlotFlankLogo(non_event_logo, method = "bits")
# where the two sets differ most, per position and nucleotide
delta <- se_logo$downstream - non_event_logo$downstream
round(delta[, order(-apply(abs(delta), 2, max))[1:6]], 3)## 5 19 6 -3 25 4
## A 0.079 -0.075 0.013 -0.093 -0.044 -0.070
## C 0.026 0.132 0.013 0.022 0.088 0.044
## G -0.159 -0.004 0.084 0.048 -0.053 0.048
## U 0.053 -0.053 -0.110 0.022 0.009 -0.022