### 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)