### R code from vignette source 'SplicingGraphs.Rnw' ################################################### ### code chunk number 1: settings ################################################### options(width=90) .precomputed_results_path <- "precomputed_results" .loadPrecomputed <- function(objname) { filename <- paste0(objname, ".rda") path <- file.path(.precomputed_results_path, filename) tempenv <- new.env(parent=emptyenv()) load(path, envir=tempenv) updateObject(get(objname, envir=tempenv)) } ################################################### ### code chunk number 2: txdb ################################################### library(TxDb.Hsapiens.UCSC.hg19.knownGene) txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene ################################################### ### code chunk number 3: isActiveSeq ################################################### isActiveSeq(txdb)[1:25] ################################################### ### code chunk number 4: keep_chr14_only ################################################### isActiveSeq(txdb)[-match("chr14", names(isActiveSeq(txdb)))] <- FALSE names(which(isActiveSeq(txdb))) ################################################### ### code chunk number 5: SplicingGraphs_constructor ################################################### library(SplicingGraphs) sg <- SplicingGraphs(txdb) # should take between 5 and 10 sec on # a modern laptop sg ################################################### ### code chunk number 6: load_sg_with_bubbles ################################################### ## We load a precomputed 'sg' that contains all the bubbles in ## 'sg@.bubbles_cache'. sg_with_bubbles <- .loadPrecomputed("sg_with_bubbles") sg@.bubbles_cache <- sg_with_bubbles@.bubbles_cache ## Replace NAs with FALSE in circularity flag (because having the flag set ## to NA instead of FALSE (or vice-versa) is not considered a significant ## difference between the 2 objects): isCircular(sg) <- isCircular(sg) %in% TRUE isCircular(sg_with_bubbles) <- isCircular(sg_with_bubbles) %in% TRUE if (!identical(sg, sg_with_bubbles)) stop("'sg' is not identical to precomputed version") ################################################### ### code chunk number 7: names_sg ################################################### names(sg)[1:20] ################################################### ### code chunk number 8: seqnames_sg ################################################### seqnames(sg)[1:20] strand(sg)[1:20] table(strand(sg)) ################################################### ### code chunk number 9: elementNROWS_sg ################################################### elementNROWS(sg)[1:20] ################################################### ### code chunk number 10: sg_double_bracket_3183 ################################################### sg[["3183"]] ################################################### ### code chunk number 11: mcols_sg_double_bracket_3183 ################################################### mcols(sg[["3183"]]) ################################################### ### code chunk number 12: mcols_sg_double_bracket_3183 ################################################### mcols(sg[["3183"]])$txpath ################################################### ### code chunk number 13: txpath_sg_3183 ################################################### txpath(sg[["3183"]]) ################################################### ### code chunk number 14: plotTranscripts_sg_3183 (eval = FALSE) ################################################### ## plotTranscripts(sg[["3183"]]) ################################################### ### code chunk number 15: plotTranscripts_sg_3183_as_pdf ################################################### pdf("3183-transcripts.pdf", width=6, height=3) plotTranscripts(sg[["3183"]]) dev.off() ################################################### ### code chunk number 16: unlist_sg ################################################### ex_by_tx <- unlist(sg) head(names(ex_by_tx)) ################################################### ### code chunk number 17: unlist_sg ################################################### ex_by_tx[names(ex_by_tx) %in% c("10001", "100129075")] ################################################### ### code chunk number 18: single_bracket_subsetting ################################################### sg[strand(sg) == "-"] sg[1:20] tail(sg) # equivalent to 'sg[tail(seq_along(sg))]' sg["3183"] ################################################### ### code chunk number 19: sgedges_sg_3183 ################################################### sgedges(sg["3183"]) ################################################### ### code chunk number 20: sgnodes_sg_3183 ################################################### sgnodes(sg["3183"]) ################################################### ### code chunk number 21: sgedgesByGene_sg ################################################### edges_by_gene <- sgedgesByGene(sg) ################################################### ### code chunk number 22: edges_by_gene_3183 ################################################### edges_by_gene[["3183"]] ################################################### ### code chunk number 23: plot_sg_3183 (eval = FALSE) ################################################### ## plot(sg["3183"]) ## plot(sgraph(sg["3183"], tx_id.as.edge.label=TRUE)) ################################################### ### code chunk number 24: plot_sg_3183_as_pdf ################################################### pdf("3183-sgraph.pdf", width=3, height=5) plot(sgraph(sg["3183"])) dev.off() pdf("3183-sgraph-labelled.pdf", width=3, height=5) plot(sgraph(sg["3183"], tx_id.as.edge.label=TRUE)) dev.off() ################################################### ### code chunk number 25: plot_sg_7597_as_pdf ################################################### pdf("7597-sgraph-labelled.pdf", width=5, height=5) plot(sgraph(sg["7597"], tx_id.as.edge.label=TRUE)) dev.off() ################################################### ### code chunk number 26: bubbles_sg_7597 ################################################### bubbles(sg["7597"]) ################################################### ### code chunk number 27: ASCODE2DESC ################################################### codes <- bubbles(sg["7597"])$AScode data.frame(AScode=codes, description=ASCODE2DESC[codes], row.names=NULL) codes <- bubbles(sg["10202"])$AScode data.frame(AScode=codes, description=ASCODE2DESC[codes], row.names=NULL) ################################################### ### code chunk number 28: AScode_summary ################################################### AScode_list <- lapply(seq_along(sg), function(i) bubbles(sg[i])$AScode) names(AScode_list) <- names(sg) AScode_table <- table(unlist(AScode_list)) AScode_table <- sort(AScode_table, decreasing=TRUE) AScode_summary <- data.frame(AScode=names(AScode_table), NbOfEvents=as.vector(AScode_table), Description=ASCODE2DESC[names(AScode_table)]) nrow(AScode_summary) ################################################### ### code chunk number 29: head_AScode_summary ################################################### head(AScode_summary, n=10) ################################################### ### code chunk number 30: nb_bubbles_per_gene ################################################### nb_bubbles_per_gene <- elementNROWS(AScode_list) ################################################### ### code chunk number 31: head_sort_nb_bubbles_per_gene ################################################### head(sort(nb_bubbles_per_gene, decreasing=TRUE)) ################################################### ### code chunk number 32: nb_unique_bubbles_per_gene ################################################### nb_unique_bubbles_per_gene <- elementNROWS(unique(CharacterList(AScode_list))) ################################################### ### code chunk number 33: head_sort_nb_unique_bubbles_per_gene ################################################### head(sort(nb_unique_bubbles_per_gene, decreasing=TRUE)) ################################################### ### code chunk number 34: plot_sg_4287_as_pdf ################################################### pdf("4287-sgraph.pdf", width=5, height=5) plot(sgraph(sg["4287"])) dev.off() ################################################### ### code chunk number 35: bam_files ################################################### library(RNAseqData.HNRNPC.bam.chr14) bam_files <- RNAseqData.HNRNPC.bam.chr14_BAMFILES names(bam_files) # the names of the runs ################################################### ### code chunk number 36: quickBamFlagSummary ################################################### quickBamFlagSummary(bam_files[1], main.groups.only=TRUE) ################################################### ### code chunk number 37: readGAlignmentPairs ################################################### param <- ScanBamParam(which=GRanges("chr14", IRanges(1, 20000000))) galp <- readGAlignmentPairs(bam_files[1], param=param) length(galp) # nb of alignment pairs galp ################################################### ### code chunk number 38: ScanBamParam ################################################### flag0 <- scanBamFlag(isSecondaryAlignment=FALSE, isNotPassingQualityControls=FALSE, isDuplicate=FALSE) param0 <- ScanBamParam(flag=flag0) ################################################### ### code chunk number 39: assignReads (eval = FALSE) ################################################### ## ## The following loop takes about 7 minutes on a modern laptop/desktop... ## for (i in seq_along(bam_files)) { ## bam_file <- bam_files[i] ## cat("Processing run ", names(bam_file), " ... ", sep="") ## galp <- readGAlignmentPairs(bam_file, use.names=TRUE, param=param0) ## sg <- assignReads(sg, galp, sample.name=names(bam_file)) ## cat("OK\n") ## } ################################################### ### code chunk number 40: load_sg_with_reads ################################################### sg_with_reads <- .loadPrecomputed("sg_with_reads") ## Remove the reads from 'sg_with_reads' and compare with 'sg'. if (!isTRUE(all.equal(removeReads(sg_with_reads), sg))) stop("after removal of the hits, precomputed version of 'sg_with_reads' ", "is not identical to 'sg'") sg <- sg_with_reads ################################################### ### code chunk number 41: sg_3183_assignments ################################################### edges_by_tx <- sgedgesByTranscript(sg["3183"], with.hits.mcols=TRUE) edge_data <- mcols(unlist(edges_by_tx)) colnames(edge_data) ################################################### ### code chunk number 42: head_edge_data_sgedge_id_ERR127306.hits ################################################### head(edge_data[ , c("sgedge_id", "ERR127306.hits")]) ################################################### ### code chunk number 43: load_reads_in_3183_region ################################################### param <- ScanBamParam(flag=flag0, which=range(unlist(sg[["3183"]]))) reads <- readGAlignmentPairs(bam_files[1], use.names=TRUE, param=param) junction_reads <- reads[njunc(first(reads)) + njunc(last(reads)) != 0L] ################################################### ### code chunk number 44: plotTranscripts_sg_3183_and_reads (eval = FALSE) ################################################### ## plotTranscripts(sg[["3183"]], reads=junction_reads, from=21675000, to=21702000) ################################################### ### code chunk number 45: plotTranscripts_sg_3183_and_reads_as_pdf ################################################### pdf("3183-transcripts-and-reads.pdf", width=12, height=12) plotTranscripts(sg[["3183"]], reads=junction_reads, from=21675000, to=21702000) dev.off() ################################################### ### code chunk number 46: plotTranscripts_sg_3183_and_reads_zoom (eval = FALSE) ################################################### ## plotTranscripts(sg[["3183"]], reads=junction_reads, from=21698400, to=21698600) ################################################### ### code chunk number 47: plotTranscripts_sg_3183_and_reads_zoom_as_pdf ################################################### pdf("3183-transcripts-and-reads-zoom.pdf", width=12, height=12) plotTranscripts(sg[["3183"]], reads=junction_reads, from=21698400, to=21698600) dev.off() ################################################### ### code chunk number 48: countReads_sg ################################################### sg_counts <- countReads(sg) ################################################### ### code chunk number 49: head_sg_counts ################################################### dim(sg_counts) head(sg_counts[1:5]) ################################################### ### code chunk number 50: sapply_sg_counts ################################################### sapply(sg_counts[-(1:2)], sum) ################################################### ### code chunk number 51: SessionInfo ################################################### sessionInfo()