### R code from vignette source 'vignettes/GenomicFeatures/inst/doc/spliceGraph.Rnw' ################################################### ### code chunk number 1: spliceGraph.Rnw:32-33 ################################################### reCnt <- FALSE ################################################### ### code chunk number 2: loadTxdb ################################################### library("TxDb.Hsapiens.UCSC.hg18.knownGene") txdb <- TxDb.Hsapiens.UCSC.hg18.knownGene ################################################### ### code chunk number 3: deactChr ################################################### activeChr <- ! isActiveSeq(txdb) activeChr[names(activeChr) == "chr19"] <- TRUE isActiveSeq(txdb) <- activeChr ################################################### ### code chunk number 4: loadGenomicFeatures ################################################### library("GenomicFeatures") sg <- spliceGraph(txdb) ################################################### ### code chunk number 5: exbyedges ################################################### sg["1008"] ################################################### ### code chunk number 6: spliceevents ################################################### head(metadata(sg)$spliceEvents, n=2) ################################################### ### code chunk number 7: bubbles ################################################### spliceEvents <- metadata(sg)$spliceEvents sp <- split(spliceEvents$bubbleCode, factor(spliceEvents$eventType)) eventSum <- data.frame(event=names(sp), Nr=elementLengths(sp), row.names=NULL) eventSum <- eventSum[order(eventSum$Nr, decreasing=TRUE), ] ################################################### ### code chunk number 8: splicEventSum ################################################### head(eventSum) ################################################### ### code chunk number 9: varietyOfEvents ################################################### allEvents <- lapply(spliceEvents$code, paste, collapse=" ") allEvents <- unlist(allEvents, use.names=FALSE) nrOfdiffEvents <- length(unique(allEvents)) nrOfdiffEvents ################################################### ### code chunk number 10: freqEvents ################################################### events <- factor(unlist(allEvents, use.names=FALSE)) occurance <- tabulate(events) idx <- !duplicated(unlist(allEvents, use.names=FALSE)) uniCodes <- spliceEvents$code[idx] size <- elementLengths(uniCodes) fn <- "power-law.png" png(fn, width=7, height=7, res=300, units="in") plot(log(size), log(occurance), frame.plot=FALSE, pch=16, cex=0.7, xlab="log( size of alternative splicing events )", ylab="log( frequency of alternative splicing events )", main="Complexity of splice events (Chr 19)") dev.off() ################################################### ### code chunk number 11: spliceGraph.Rnw:435-439 ################################################### idx <- ! duplicated(spliceEvents$bubbleCode) nrOfBubs <- elementLengths(split(spliceEvents$bubbleCode[idx], spliceEvents$gn[idx])) nrOfBubs <- nrOfBubs[order(nrOfBubs,decreasing=TRUE)] ################################################### ### code chunk number 12: spliceGraph.Rnw:444-445 ################################################### head(nrOfBubs) ################################################### ### code chunk number 13: spliceGraph.Rnw:451-456 ################################################### uniEventsPerGn <- lapply(split(allEvents, spliceEvents$gn), unique) nrOfUniEvents <- elementLengths(uniEventsPerGn) nrOfUniEvents <- nrOfUniEvents[order(nrOfUniEvents, decreasing=TRUE)] ################################################### ### code chunk number 14: spliceGraph.Rnw:462-463 ################################################### head(nrOfUniEvents) ################################################### ### code chunk number 15: origExToDJEx ################################################### sg.flat <- unlist(sg, use.names=FALSE) orig.Ex <- values(sg.flat)[["exon_ids"]] disJ.Ex <- values(sg.flat)[["disJ_exon_id"]] newExNames <- rep(disJ.Ex, elementLengths(orig.Ex)) origExToNewEx <- data.frame(orig.Ex=unlist(orig.Ex), disJ.Ex=newExNames, row.names=NULL) head(origExToNewEx) ################################################### ### code chunk number 16: getGeneToEdge ################################################### gnIDs <- values(sg)[["gene_id"]] edgeIDs <- names(sg) gnToEdge <- setNames(gnIDs, edgeIDs) ################################################### ### code chunk number 17: countingReads ################################################### if(reCnt) { require("Rsamtools") bamPath <- "/shared/labs/EDI/users/mfitzgib/Solexa" fls <- sub(".bai$", "", list.files(bamPath, recursive=TRUE, pattern="accepted_hits.*bai$", full=TRUE)) fls <- fls[c(1:3, 62:64)] bfs <- BamFileList(fls) names(bfs) <- gsub("/shared/labs/EDI/users/mfitzgib/Solexa/tophat_", "",fls) } ################################################### ### code chunk number 18: CountReads ################################################### if(reCnt) { library(parallel) resEx.byEdge <- summarizeOverlaps(features = sg, reads = bfs, mc.cores = getOption("mc.cores", 3L)) cD.exByEdge <- assays(resEx.byEdge)$counts colnames(cD.exByEdge) <- sub("/accepted_hits.bam", "", colnames(cD.exByEdge)) save(cD.exByEdge, file="cD.exByEdge-SG-Vig.Rda") } else { fn <- system.file("extdata", "cD.exByEdge-SG-Vig.Rda", package="GenomicFeatures") load(fn) } ################################################### ### code chunk number 19: readsPerEx ################################################### exsByGenes <- exonsBy(txdb, "gene") gnToEx <- rep(names(exsByGenes), elementLengths(exsByGenes)) exsByGenes.flat <- unlist(exsByGenes, use.names=FALSE) # remove duplicates notDupl <- !duplicated(values(exsByGenes.flat)[["exon_id"]]) exsByGenes.flat <- exsByGenes.flat[notDupl] names(exsByGenes.flat) <- values(exsByGenes.flat)[["exon_id"]] gnToEx <- gnToEx[notDupl] names(gnToEx) <- names(exsByGenes.flat) ################################################### ### code chunk number 20: edgeVsEx ################################################### if(reCnt) { res.exsByGenes <- summarizeOverlaps(features = exsByGenes.flat, reads = bfs, mc.cores = getOption("mc.cores", 3L)) cD.exsByGenes <- assays(res.exsByGenes)$counts colnames(cD.exsByGenes) <- sub("/accepted_hits.bam", "", colnames(cD.exsByGenes)) save(cD.exsByGenes, file="cD.exsByGenes-SG-Vig.Rda") } else { fn <- system.file("extdata", "cD.exsByGenes-SG-Vig.Rda", package="GenomicFeatures") load(fn) } ################################################### ### code chunk number 21: prepDiffExp ################################################### edIds <- names(sg) gnIds.edges <- gnToEdge[rownames(cD.exByEdge)] nrOfEdgesPerGns <- elementLengths(split(edIds, factor(gnIds.edges))) sglEdgeGns <- names(nrOfEdgesPerGns)[nrOfEdgesPerGns < 2] exIds <- names(exsByGenes.flat) gnIds.exs <- gnToEx[names(exsByGenes.flat)] nrOfExsPerGns <- elementLengths(split(exIds, factor(gnIds.exs))) sglExGns <- names(nrOfExsPerGns)[nrOfExsPerGns < 2] ################################################### ### code chunk number 22: subsetCounts ################################################### cD.exByEdge <- cD.exByEdge[, colnames(cD.exsByGenes)] ################################################### ### code chunk number 23: DEXSeq ################################################### benign <- grepl("Benign",colnames(cD.exByEdge ) ) design <- factor(ifelse(benign, "Benign", "SOC"), levels=c("Benign", "SOC")) names(design) <- ifelse(benign, "Benign", "SOC") design ################################################### ### code chunk number 24: resDexSeq ################################################### sg.flat <- unlist(sg, use.names=FALSE) edgeIDs <- rep(names(sg), elementLengths(sg)) rle <- strand(sg.flat) start <- start(sg.flat) end <- end(sg.flat) temp <- data.frame(start, end, edgeIDs, strand=rep(rle@values, rle@lengths), chr=rep(seqnames(sg.flat)@values, seqnames(sg.flat)@lengths)) temp <- temp[order(temp$edgeIDs, -temp$end), ] end <- temp$end[!duplicated(temp$edgeIDs)] names(end) <- temp$edgeIDs[!duplicated(temp$edgeIDs)] temp <- temp[order(temp$edgeIDs, temp$start), ] edgeAnnot <- temp[!duplicated(temp$edgeIDs),] rownames(edgeAnnot) <- edgeAnnot$edgeIDs edgeAnnot$end <- end[rownames(edgeAnnot)] edgeAnnot <- edgeAnnot[order(gnToEdge[edgeAnnot$edgeIDs], edgeAnnot$start), ] ################################################### ### code chunk number 25: newData ################################################### library(DEXSeq) eData.Edge <- newExonCountSet(countData = cD.exByEdge[edgeAnnot$edgeIDs,] , design = design, geneIDs = gnToEdge[edgeAnnot$edgeIDs], exonIDs = edgeAnnot$edgeIDs, exonIntervals = edgeAnnot) idx <- ! gnToEdge[edgeAnnot$edgeIDs] %in% sglEdgeGns eData.Edge <- eData.Edge[idx, ] ################################################### ### code chunk number 26: table ################################################### exAnnot <- data.frame(chr=rep(seqnames(exsByGenes.flat)@values, seqnames(exsByGenes.flat)@lengths), start=start(exsByGenes.flat), end=end(exsByGenes.flat), strand=rep(strand(exsByGenes.flat)@values, strand(exsByGenes.flat)@lengths), row.names=values(exsByGenes.flat)[["exon_id"]]) ################################################### ### code chunk number 27: qqq ################################################### eData.Ex <- newExonCountSet(countData = cD.exsByGenes, design = design, geneIDs = gnToEx[rownames(cD.exsByGenes)], exonIDs = rownames(cD.exsByGenes), exonIntervals=exAnnot[rownames(cD.exsByGenes),]) eData.Ex <- eData.Ex[! gnToEx[rownames(cD.exsByGenes)] %in% sglExGns, ] ################################################### ### code chunk number 28: ttEx ################################################### eData.Ex <- estimateSizeFactors(eData.Ex) eData.Ex <- estimateDispersions(eData.Ex) eData.Ex <- fitDispersionFunction(eData.Ex) eData.Ex <- testForDEU(eData.Ex) eData.Ex <- estimatelog2FoldChanges(eData.Ex) tt.Ex <- DEUresultTable(eData.Ex) ################################################### ### code chunk number 29: ttEdge ################################################### eData.Edge <- estimateSizeFactors(eData.Edge) eData.Edge <- estimateDispersions(eData.Edge) eData.Edge <- fitDispersionFunction(eData.Edge) eData.Edge <- testForDEU(eData.Edge) eData.Edge <- estimatelog2FoldChanges(eData.Edge) tt.Edge <- DEUresultTable(eData.Edge) ################################################### ### code chunk number 30: plots ################################################### plotDisp <- function(eData, case="") { meanvalues <- rowMeans(counts(eData)) plot(meanvalues, fData(eData)$dispBeforeSharing, main = paste(case, " mean vs CR dispersion", sep=":"), frame.plot=FALSE, pch=16, cex=0.8, log = "xy", ylab="Dispersion", xlab="Mean counts") x <- 0.01:max(meanvalues) y <- eData@dispFitCoefs[1] + eData@dispFitCoefs[2]/x lines(x, y, col = "purple", lwd=2) } fn <- "Disp-plots.png" png(fn, width=11, height=5.5, res=300, units="in") par(mfrow=c(1,2)) plotDisp(eData.Ex, "Exons") plotDisp(eData.Edge, "Edges") dev.off() ################################################### ### code chunk number 31: se ################################################### library(DESeq) cds.Ex <- newCountDataSet(countData=cD.exsByGenes, conditions=design) cds.Ex <- cds.Ex[gnToEx[rownames(cD.exsByGenes)] %in% sglExGns,] sel <- gnToEdge[rownames(cD.exByEdge)] %in% sglEdgeGns cds.Edge <- newCountDataSet(countData=cD.exByEdge[sel,], conditions=design) ################################################### ### code chunk number 32: difExEx ################################################### cds.Ex <- estimateSizeFactors( cds.Ex ) cds.Ex <- estimateDispersions( cds.Ex ) ttSingle.Ex <- nbinomTest( cds.Ex,"SOC", "Benign" ) ################################################### ### code chunk number 33: diffExEdge ################################################### cds.Edge <- estimateSizeFactors( cds.Edge ) cds.Edge <- estimateDispersions( cds.Edge ) ttSingle.Edge <- nbinomTest( cds.Edge,"SOC", "Benign" ) ################################################### ### code chunk number 34: DexSeqSum ################################################### int.DEXSeq <- c("exonID", "pvalue", "padjust", "meanBase", "log2fold(SOC/Benign)") int.DESeq <- c("id", "pval", "padj", "baseMean", "log2FoldChange") temp <- tt.Edge[, int.DEXSeq] colnames(temp) <- int.DESeq finTT.Edge <- rbind(ttSingle.Edge[, int.DESeq], temp) finTT.Edge$gnID <- gnToEdge[finTT.Edge$id] temp <- tt.Ex[, int.DEXSeq] colnames(temp) <- int.DESeq finTT.Ex <- rbind(ttSingle.Ex[, int.DESeq], temp) finTT.Ex$gnID <- gnToEx[finTT.Ex$id] ################################################### ### code chunk number 35: compEdgeEx ################################################### sum(cD.exByEdge) sum(cD.exsByGenes) sum(cD.exsByGenes*100)/sum(cD.exByEdge) ################################################### ### code chunk number 36: MAplots ################################################### fn <- "MA-plots.png" png(fn, width=11, height=5.5, res=300, units="in") par(mfrow=c(1,2)) plot(finTT.Ex$baseMean, finTT.Ex$log2FoldChange, log = "x", col = ifelse(finTT.Ex$padj < 0.1, "purple", "black"), ylim = c(-5, 5), main = "Exons MvsA", pch=16, cex=0.7, frame.plot=FALSE) plot(finTT.Edge$baseMean, finTT.Edge$log2FoldChange, col = ifelse(finTT.Edge$padj < 0.1, "purple", "black"), ylim = c(-5, 5), main = "Edges MvsA", pch=16, cex=0.7, log = "x", frame.plot=FALSE) dev.off() ################################################### ### code chunk number 37: sigDiff ################################################### ex.test <- nrow(finTT.Ex) edge.test <- nrow(finTT.Edge) sig.diff.ex <- sum(finTT.Ex$padj < 0.05, na.rm=TRUE) sig.diff.edge <- sum(finTT.Edge$padj < 0.05, na.rm=TRUE) sig.diff.ex*100/ex.test sig.diff.edge*100/edge.test ################################################### ### code chunk number 38: pcutSum ################################################### p.cuts <- 10^seq( 0, -5, length.out=100 ) ps.Ex <- sapply(p.cuts, function(p.cut) { sum(finTT.Ex$pval < p.cut, na.rm=TRUE)*100/ length(finTT.Ex$pval) }) ps.Edge <- sapply(p.cuts, function(p.cut) { sum(finTT.Edge$pval < p.cut, na.rm=TRUE)*100/ length(finTT.Edge$pval) }) fn <- "Comparison-edges-exons.png" png(fn, width=7, height=7, res=300, units="in") plot(ps.Edge, p.cuts, type="l", col="purple", xlab="% of differentially expressed elements", ylab="P-value cut off", lwd=2, frame.plot=FALSE, log="xy", main="Exons versus edges") lines(ps.Ex, p.cuts, lty=2, lwd=2) grid(lwd=2) legend("bottomright", legend=c("Edges", "Exons"), col=c("purple", "black"), lty=c(1,2), lwd=2, border=NA, box.col=NA, bg=NA) dev.off() ################################################### ### code chunk number 39: pp ################################################### library(xtable) T <- xtable(head(finTT.Ex[order(finTT.Ex$padj),], n=10), caption="Top table of the exon model", label="table:ex", display=c("s","d","E","E","f","f","d")) print(T, include.rownames=FALSE, table.placement="H") ################################################### ### code chunk number 40: rr ################################################### library(xtable) T <- xtable(head(finTT.Edge[order(finTT.Edge$padj),], n=10), caption="Top table of the edge model", label="table:edge", display=c("s","d","E","E","f","f","d")) print(T, include.rownames=FALSE, table.placement="H") ################################################### ### code chunk number 41: ww ################################################### plotExp <- function(fn, eDat) { png(fn, width=11, height=5.5, res=300, units="in") par(mfrow=c(1,2)) COL <- c("#3399FF", "#FF3333") plotDEXSeq(eDat, geneID="90522", cex.axis = 1.2, cex = 1.3, lwd = 2, legend = TRUE, color=COL, color.samples=COL[design], splicing=TRUE) plotDEXSeq(eDat, geneID="90522", expression = FALSE, norCounts = TRUE, cex.axis = 1.2, cex = 1.3, lwd = 2, legend = TRUE, color=COL, color.samples=COL[design], splicing=TRUE) dev.off() } fn <- "Edge-plot.png" plotExp(fn, eData.Edge) fn <- "Ex-plot.png" plotExp(fn, eData.Ex) ################################################### ### code chunk number 42: spliceGraph.Rnw:1143-1152 ################################################### bubID <- values(sg)[["bubble_ids"]] eid <- rep(names(sg), elementLengths(bubID)) rownames(spliceEvents) <- spliceEvents$bubbleCode bubID <- unlist(bubID, use.names=FALSE) eidToBubID <- data.frame(bubID, eid, eventType=spliceEvents[bubID,]$eventType) ################################################### ### code chunk number 43: spliceGraph.Rnw:1158-1171 ################################################### eidToBubID <- eidToBubID[!is.na(eidToBubID$eventType) & !is.na(eidToBubID$eventType ), ] finTT.Edge <- finTT.Edge[finTT.Edge$id %in% eidToBubID$eid, ] rownames(finTT.Edge) <- as.character(finTT.Edge$id) eidToBubID <- eidToBubID[eidToBubID$eid %in% rownames(finTT.Edge),] Res <- data.frame(eidToBubID, finTT.Edge[as.character(eidToBubID$eid),]) Res <- Res[order(Res$pval), ] ################################################### ### code chunk number 44: rr ################################################### T <- head(Res[, c("gnID","padj", "pval", "log2FoldChange", "eventType")], n=15) colnames(T) <- c("gnID", "padj", "pval", "log2FC", "event") T <- xtable(T, caption="The 15 most significant annotated splice events affacted by differential expression.", label="table:events", display=c("s", "s","E", "E","f","s")) align(T) <- "lrrrrp{5cm}" print(T, include.rownames=FALSE, table.placement="H") ################################################### ### code chunk number 45: SessionInfo ################################################### sessionInfo()