### R code from vignette source 'tutorial.Rnw'

###################################################
### code chunk number 1: setup
###################################################
options(width=50)
library(isoformExprTutorial)


###################################################
### code chunk number 2: aldoa
###################################################
aldoa_eg <- org.Hs.egSYMBOL2EG$ALDOA
txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene
aldoa_gr <- exons(txdb, vals = list(gene_id = aldoa_eg),
                  columns = c("tx_id", "gene_id"))
aldoa_range <- range(aldoa_gr)


###################################################
### code chunk number 3: forming-models
###################################################
aldoa_vals <- values(aldoa_gr)
values(aldoa_gr) <- NULL
tx <- multisplit(aldoa_gr, aldoa_vals$tx_id)
tx_to_val <- match(names(tx), unlist(aldoa_vals$tx_id))
values(tx)$gene_id <- 
  rep(unlist(aldoa_vals$gene_id), 
      elementLengths(aldoa_vals$tx_id))[tx_to_val]
values(tx)$tx_id <- names(tx)


###################################################
### code chunk number 4: solution-exonsBy
###################################################
exons_grl <- exonsBy(txdb)
ans <- subsetByOverlaps(exons_grl, aldoa_range)
values(ans)$tx_id <- names(ans)
tx_gr <- transcripts(txdb, columns = c("tx_id", "gene_id"))
values(ans)$gene_id <- 
  drop(values(tx_gr)$gene_id)[match(names(ans), 
                                    values(tx_gr)$tx_id)]


###################################################
### code chunk number 5: merge-gene-symbols
###################################################
gene_id_keys <- 
  values(tx)$gene_id[!is.na(values(tx)$gene_id)]
tx_gene_sym <- rep.int(NA, length(tx))
tx_gene_sym[!is.na(values(tx)$gene_id)] <- 
  unlist(mget(gene_id_keys, org.Hs.egSYMBOL, 
              ifnotfound = NA),
         use.names = FALSE)
values(tx)$gene_sym <- tx_gene_sym


###################################################
### code chunk number 6: ccds-only
###################################################
tx <- tx[c(2, 4, 5, 7)]


###################################################
### code chunk number 7: normal-bam
###################################################
extdatadir <- system.file("extdata", 
                          package = "isoformExprTutorial")
files <- tools::list_files_with_exts(extdatadir, "bam")
names(files) <- tools::file_path_sans_ext(basename(files))
bamFiles <- Rsamtools::BamFileList(files)
bam <- bamFiles$normal


###################################################
### code chunk number 8: read-ga
###################################################
param <- ScanBamParam(tag = "XS", which = aldoa_range)
ga <- readGappedAlignments(path(bam), 
                           use.names = TRUE, 
                           param = param)


###################################################
### code chunk number 9: read-ga2
###################################################
reads <- grglist(ga)
metadata(reads)$bamfile <- bam


###################################################
### code chunk number 10: solution-metadata
###################################################
metadata(reads)$param <- param


###################################################
### code chunk number 11: elementGaps-real
###################################################
elementGaps <- function(x) {
  x_flat <- unlist(x, use.names = FALSE)
  egaps <- gaps(ranges(x))
  first_segment <- start(PartitioningByWidth(x))
  sn <- seqnames(x_flat)[first_segment][togroup(egaps)]
  strand <- strand(x_flat)[first_segment][togroup(egaps)]
  relist(GRanges(sn, unlist(egaps, use.names = FALSE), 
                 strand, seqlengths = seqlengths(x)), 
         egaps)
}


###################################################
### code chunk number 12: splices
###################################################
splices <- elementGaps(reads)
values(splices)$XS <- values(reads)$XS


###################################################
### code chunk number 13: elementGaps-simple
###################################################
elementGaps <- function(reads) {
  psetdiff(unlist(range(reads), use.names=FALSE), reads)
}


###################################################
### code chunk number 14: elementGaps-optimized
###################################################
elementGaps <- function(x) {
  x_flat <- unlist(x, use.names = FALSE)
  egaps <- gaps(ranges(x))
  first_segment <- start(PartitioningByWidth(x))
  sn <- seqnames(x_flat)[first_segment][togroup(egaps)]
  strand <- strand(x_flat)[first_segment][togroup(egaps)]
  relist(GRanges(sn, unlist(egaps, use.names = FALSE), 
                 strand, seqlengths = seqlengths(x)), 
         egaps)
}


###################################################
### code chunk number 15: solution-partitioning
###################################################
tx_part <- PartitioningByWidth(tx)
tx_flat <- unlist(tx, use.names = FALSE)
start(tx_flat)[start(tx_part)]


###################################################
### code chunk number 16: pair-reads
###################################################
pairs <- split(unlist(reads, use.names=FALSE),
               factor(names(reads)[togroup(reads)], 
                      unique(names(reads))))
metadata(pairs) <- metadata(reads)


###################################################
### code chunk number 17: pair-reads-xs
###################################################
xs <- values(reads)$XS
has_xs <- !is.na(xs)
pair_xs <- setNames(rep.int(NA, length(pairs)), 
                    names(pairs))
pair_xs[names(reads)[has_xs]] <- xs[has_xs]
values(pairs)$XS <- unname(pair_xs)


###################################################
### code chunk number 18: resolve-xs
###################################################
xs <- values(pairs)$XS
strand <- ifelse(!is.na(xs) & xs != "?", xs, "*")
strand(pairs) <- relist(Rle(strand, elementLengths(pairs)), 
                        pairs)


###################################################
### code chunk number 19: pairReadRanges
###################################################
pairReadRanges <- function(reads) {
pairs <- split(unlist(reads, use.names=FALSE),
               factor(names(reads)[togroup(reads)], 
                      unique(names(reads))))
metadata(pairs) <- metadata(reads)
xs <- values(reads)$XS
has_xs <- !is.na(xs)
pair_xs <- setNames(rep.int(NA, length(pairs)), 
                    names(pairs))
pair_xs[names(reads)[has_xs]] <- xs[has_xs]
values(pairs)$XS <- unname(pair_xs)
  pairs
}  


###################################################
### code chunk number 20: strandFromXS
###################################################
strandFromXS <- function(pairs) {
xs <- values(pairs)$XS
strand <- ifelse(!is.na(xs) & xs != "?", xs, "*")
strand(pairs) <- relist(Rle(strand, elementLengths(pairs)), 
                        pairs)
  pairs
}


###################################################
### code chunk number 21: pair-resolve-splices
###################################################
splices <- pairReadRanges(splices)
splices <- strandFromXS(splices)


###################################################
### code chunk number 22: readReadRanges
###################################################
readReadRanges <- function(bam) {
param <- ScanBamParam(tag = "XS", which = aldoa_range)
ga <- readGappedAlignments(path(bam), 
                           use.names = TRUE, 
                           param = param)
reads <- grglist(ga)
metadata(reads)$bamfile <- bam
splices <- elementGaps(reads)
values(splices)$XS <- values(reads)$XS
pairs <- pairReadRanges(reads)
pairs <- strandFromXS(pairs)
splices <- pairReadRanges(splices)
splices <- strandFromXS(splices)
values(pairs)$splices <- splices
pairs
}


###################################################
### code chunk number 23: read-tumor-normal
###################################################
normal <- readReadRanges(bamFiles$normal)
tumor <- readReadRanges(bamFiles$tumor)


###################################################
### code chunk number 24: findOverlaps
###################################################
hits <- findOverlaps(pairs, tx)


###################################################
### code chunk number 25: extract-hits
###################################################
hit_pairs <- ranges(pairs)[queryHits(hits)]
hit_splices <- ranges(splices)[queryHits(hits)]
hit_tx <- ranges(tx)[subjectHits(hits)]


###################################################
### code chunk number 26: compatibility-checking
###################################################
read_within <- 
  elementLengths(setdiff(hit_pairs, hit_tx)) == 0L
tx_within <- 
  elementLengths(intersect(hit_tx, hit_splices)) == 0L
compatible <- read_within & tx_within


###################################################
### code chunk number 27: uniquely-mapped
###################################################
compat_hits <- hits[compatible]
reads_unique <- tabulate(queryHits(compat_hits), 
                         queryLength(compat_hits)) == 1L
unique <- logical(length(hits))
unique[compatible] <- reads_unique[queryHits(compat_hits)]


###################################################
### code chunk number 28: overlap-annotation
###################################################
strand_specific <- 
  all(strand(pairs) != "*")[queryHits(hits)]
values(hits) <- DataFrame(strand_specific,
                          compatible,
                          unique)


###################################################
### code chunk number 29: solution-spliced
###################################################
values(hits)$spliced <- 
  (elementLengths(pairs) > 2)[queryHits(hits)]


###################################################
### code chunk number 30: grkey
###################################################
gr2key <- function(x) {
  paste(seqnames(x), start(x), end(x), strand(x), 
        sep = ":")
}


###################################################
### code chunk number 31: key2gr
###################################################
key2gr <- function(x, ...) {
  key_mat <- matrix(unlist(strsplit(x, ":", fixed=TRUE)), 
                    nrow = 4)
  GRanges(key_mat[1,],
          IRanges(as.integer(key_mat[2,]), 
                  as.integer(key_mat[3,])),
          key_mat[4,], ...)
}


###################################################
### code chunk number 32: txkey
###################################################
introns <- elementGaps(tx)
introns_flat <- unlist(introns, use.names = FALSE)
tx_keys <- gr2key(introns_flat)


###################################################
### code chunk number 33: splice-summary
###################################################
splices_flat <- unlist(splices, use.names = FALSE)
splice_table <- table(gr2key(splices_flat))
splice_summary <- 
  key2gr(names(splice_table), 
         score = as.integer(splice_table),
         novel = !names(splice_table) %in% tx_keys,
         seqlengths = seqlengths(splices))


###################################################
### code chunk number 34: aggregate-isoforms
###################################################
countByTx <- function(x) {
  tabulate(subjectHits(hits)[x], subjectLength(hits))
}
compatible_strand <- 
  countByTx(with(values(hits), 
                 compatible & strand_specific))
counts <- DataFrame(compatible_strand,
                    lapply(values(hits)[-1], countByTx))


###################################################
### code chunk number 35: findIsoformOverlaps
###################################################
findIsoformOverlaps <- function(pairs) {
  splices <- values(pairs)$splices
hits <- findOverlaps(pairs, tx)
hit_pairs <- ranges(pairs)[queryHits(hits)]
hit_splices <- ranges(splices)[queryHits(hits)]
hit_tx <- ranges(tx)[subjectHits(hits)]
read_within <- 
  elementLengths(setdiff(hit_pairs, hit_tx)) == 0L
tx_within <- 
  elementLengths(intersect(hit_tx, hit_splices)) == 0L
compatible <- read_within & tx_within
compat_hits <- hits[compatible]
reads_unique <- tabulate(queryHits(compat_hits), 
                         queryLength(compat_hits)) == 1L
unique <- logical(length(hits))
unique[compatible] <- reads_unique[queryHits(compat_hits)]
strand_specific <- 
  all(strand(pairs) != "*")[queryHits(hits)]
values(hits) <- DataFrame(strand_specific,
                          compatible,
                          unique)
  hits
}


###################################################
### code chunk number 36: countIsoformHits
###################################################
countIsoformHits <- function(hits) {
countByTx <- function(x) {
  tabulate(subjectHits(hits)[x], subjectLength(hits))
}
compatible_strand <- 
  countByTx(with(values(hits), 
                 compatible & strand_specific))
counts <- DataFrame(compatible_strand,
                    lapply(values(hits)[-1], countByTx))
  counts
}


###################################################
### code chunk number 37: summarizeSplices
###################################################
summarizeSplices <- function(reads) {
  splices <- values(reads)$splices
splices_flat <- unlist(splices, use.names = FALSE)
splice_table <- table(gr2key(splices_flat))
splice_summary <- 
  key2gr(names(splice_table), 
         score = as.integer(splice_table),
         novel = !names(splice_table) %in% tx_keys,
         seqlengths = seqlengths(splices))
  splice_summary
}


###################################################
### code chunk number 38: process-all
###################################################
normal_hits <- findIsoformOverlaps(normal)
normal_counts <- countIsoformHits(normal_hits)
normal_splices <- summarizeSplices(normal)
tumor_hits <- findIsoformOverlaps(tumor)
tumor_counts <- countIsoformHits(tumor_hits)
tumor_splices <- summarizeSplices(tumor)


###################################################
### code chunk number 39: combine-samples
###################################################
assays <- mapply(cbind, normal_counts, tumor_counts, 
                 SIMPLIFY = FALSE)
colData <- DataFrame(tumorStatus = c("tumor", "normal"))
rownames(colData) <- colData$tumorStatus
se <- SummarizedExperiment(assays, tx, colData)


###################################################
### code chunk number 40: order-se
###################################################
uc <- assay(se, "unique")
uc_ord <- order(rowSums(uc), decreasing = TRUE)
uc_top <- uc[head(uc_ord, 2),]
fisher.test(uc_top)$estimate


###################################################
### code chunk number 41: get-unique-reads
###################################################
getUniqueReads <- function(reads, hits) {
  sel <- values(hits)$unique & 
         subjectHits(hits) %in% c(1, 4)
  reads[unique(queryHits(hits)[sel])]
}
normal_uniq <- getUniqueReads(normal, normal_hits)
tumor_uniq <- getUniqueReads(tumor, tumor_hits)
both_uniq <- mstack(normal = unlist(normal_uniq),
                    tumor = unlist(tumor_uniq))


###################################################
### code chunk number 42: combine-splices
###################################################
normal_uniq_splices <- summarizeSplices(normal_uniq)
tumor_uniq_splices <- summarizeSplices(tumor_uniq)
uniq_splices <- mstack(normal = normal_uniq_splices,
                       tumor = tumor_uniq_splices)


###################################################
### code chunk number 43: ggbio-isoform-coverage
###################################################
read_track <- autoplot(uniq_splices, geom = "arch",
                       aes(size = score, 
                           height = width / 5), 
                       color = "deepskyblue3",
                       ylab = "coverage",
                       facets = name ~ .) +
  stat_coverage(both_uniq, facets = name ~ .)
tx_16 <- keepSeqlevels(tx, "chr16")
tx_track <- autoplot(tx_16, geom = "alignment", ylab = "")
tracks(read_track, tx_track, heights = c(3, 1))


###################################################
### code chunk number 44: ggbio-zoom
###################################################
tracks(read_track,
       tx_track, heights = c(3, 1),
       xlim = c(30075000, 30080000))


###################################################
### code chunk number 45: ggbio-zoom-coverage
###################################################
cov_chr16 <- coverage(both_uniq)$chr16
roi <- range(ranges(slice(cov_chr16, 1000)))
roi <- roi + 500
tracks(read_track,
       tx_track, heights = c(3, 1),
       xlim = roi)


###################################################
### code chunk number 46: novel-splices
###################################################
all_splices <- mstack(normal = normal_splices,
                      tumor = tumor_splices)
novel_splices <- 
  all_splices[values(all_splices)$novel &
              values(all_splices)$score == 9]
uniq_novel_splices <- c(uniq_splices, novel_splices)


###################################################
### code chunk number 47: ggbio-novel-junctions
###################################################
novel_track <- autoplot(uniq_novel_splices, geom = "arch",
                        aes(size = score, 
                            height = width / 5,
                            color = novel), 
                        ylab = "coverage",
                        facets = name ~ .) +
  stat_coverage(both_uniq, facets = name ~ .)
tracks(novel_track, tx_track, heights = c(3, 1),
       xlim = roi)