################################################### ### chunk number 1: options ################################################### options(width=60) options(SweaveHooks = list(fig = function() par(mar=c(5,5,1.5,1)))) ################################################### ### chunk number 2: packages ################################################### require(yeastRNASeq) require(ShortRead) require(Genominator) data(yeastAligned) ################################################### ### chunk number 3: echo ################################################### sapply(yeastAligned, length) yeastAligned[[1]] ################################################### ### chunk number 4: biomart1 eval=FALSE ################################################### ## require(biomaRt) ## mart <- useMart("ensembl", "scerevisiae_gene_ensembl") ## attributes.gene <- c("ensembl_gene_id", "chromosome_name", "start_position", ## "end_position", "strand", "gene_biotype") ## attributes.tr <- c("ensembl_gene_id", "ensembl_transcript_id", "ensembl_exon_id", "chromosome_name", "start_position", ## "end_position", "strand", "gene_biotype", "exon_chrom_start", "exon_chrom_end", "rank") ## yAnno.gene <- getBM(attributes = attributes.gene, mart = mart) ## yAnno.tr <- getBM(attributes = attributes.tr, mart = mart) ################################################### ### chunk number 5: biomart2 eval=FALSE ################################################### ## save(yAnno.tr, yAnno.gene, file = "data/yAnno.rda") ################################################### ### chunk number 6: biomart3 ################################################### load(file = "data/yAnno.rda") ################################################### ### chunk number 7: biomart4 ################################################### head(yAnno.gene, 2) head(yAnno.tr, 2) dim(yAnno.gene) dim(yAnno.tr) subset(yAnno.gene, ensembl_gene_id == "YPR098C") subset(yAnno.tr, ensembl_gene_id == "YPR098C") length(unique(yAnno.tr$ensembl_transcript_id)) ################################################### ### chunk number 8: biomart5 ################################################### yAnno <- yAnno.tr ################################################### ### chunk number 9: eval=FALSE ################################################### ## require(rtracklayer) ## session <- browserSession() ## genome(session) <- "sacCer2" ## ucsc.sgd <- getTable(ucscTableQuery(session, "sgdGene")) ## ucsc.ens <- getTable(ucscTableQuery(session, "ensGene")) ################################################### ### chunk number 10: eval=FALSE ################################################### ## save(ucsc.sgd, ucsc.ens, file = "data/ucsc.rda") ################################################### ### chunk number 11: ################################################### load("data/ucsc.rda") ################################################### ### chunk number 12: ################################################### head(ucsc.sgd, 1) head(ucsc.ens, 1) subset(ucsc.sgd, name == "YPR098C") subset(ucsc.ens, name == "YPR098C") ################################################### ### chunk number 13: ################################################### data("geneLevelData") ################################################### ### chunk number 14: ################################################### head(geneLevelData) ################################################### ### chunk number 15: annoIR1 ################################################### chrMap <- levels(chromosome(yeastAligned[[1]])) names(chrMap) <- c(as.character(as.roman(1:16)), NA) head(chrMap) yAnno$chrom <- chrMap[yAnno$chromosome] yAnno <- yAnno[!is.na(yAnno$chrom), ] ################################################### ### chunk number 16: annoIR2 ################################################### annoByChr <- split(yAnno, yAnno$chrom) annoIR <- lapply(annoByChr, function(d) { IRanges(start = d$exon_chrom_start, end = d$exon_chrom_end) }) annoIR <- do.call(RangesList, annoIR) annoIR ################################################### ### chunk number 17: readsIR1 ################################################### toRangesList <- function(aln) { alignedByChr <- split(aln, chromosome(aln)) rngs <- lapply(alignedByChr, function(alnChr) { IRanges(start = position(alnChr), width = 1) }) do.call(RangesList, rngs) } alnAsRanges <- lapply(yeastAligned, toRangesList) alnAsRanges[[1]] ################################################### ### chunk number 18: readsOverlap ################################################### exonNames <- split(yAnno$ensembl_exon_id, yAnno$chrom) oCounts <- sapply(alnAsRanges, function(aln) { do.call(c, lapply(findOverlaps(annoIR, aln), as.table)) }) rownames(oCounts) <- do.call(c, exonNames) head(oCounts) ################################################### ### chunk number 19: GenImport ################################################### library(Genominator) chrMap <- levels(chromosome(yeastAligned[[1]])) chrMap eData <- importFromAlignedReads(yeastAligned, chrMap = chrMap, filename = "my.db", tablename = "raw", overwrite = TRUE, deleteIntermediates = FALSE) head(eData) ################################################### ### chunk number 20: ExpData1 ################################################### eData2 <- ExpData(db = "my.db", tablename = "mut_1_f", mode = "w") head(eData2, 3) ################################################### ### chunk number 21: ExpData2 ################################################### getRegion(eData, chr = 1, strand = 0, start = 10000, end = 12000) laneCounts <- summarizeExpData(eData) laneCounts summarizeExpData(eData, fxs = "MAX") ################################################### ### chunk number 22: annoExp ################################################### chrMap <- c(as.character(as.roman(1:16)), "MT", "2-micron") yAnno$chr <- match(yAnno$chromosome, chrMap) yAnno$start <- yAnno$exon_chrom_start yAnno$end <- yAnno$exon_chrom_end yAnno <- yAnno[, c("ensembl_gene_id", "ensembl_exon_id", "chr", "strand", "start", "end", "gene_biotype")] rownames(yAnno) <- yAnno$ensembl_exon_id head(yAnno,2) ################################################### ### chunk number 23: sumByAnno1 ################################################### exonCounts <- summarizeByAnnotation(eData, yAnno, ignoreStrand = TRUE) head(exonCounts, 3) ################################################### ### chunk number 24: sumByAnno2 ################################################### head(summarizeByAnnotation(eData, yAnno, bindAnno = TRUE, fxs = "COUNT"), 3) ################################################### ### chunk number 25: sumByAnno2 ################################################### geneCounts <- summarizeByAnnotation(eData, yAnno, ignoreStrand = TRUE, meta.id = "ensembl_gene_id") head(geneCounts,3) ################################################### ### chunk number 26: splitBy ################################################### exonSplit <- splitByAnnotation(eData, yAnno[1:100,], ignoreStrand = TRUE) exonSplit2 <- splitByAnnotation(eData, yAnno[1:100,], expand = TRUE, ignoreStrand = TRUE) exonSplit3 <- splitByAnnotation(eData, yAnno[1:100,], expand = TRUE, addOverStrands = TRUE, ignoreStrand = TRUE) ################################################### ### chunk number 27: applyMapped ################################################### countsPerBase <- applyMapped(exonSplit, yAnno, FUN = function(map, anno) { colSums(map, na.rm = TRUE) / (anno$end - anno$start + 1) }) ################################################### ### chunk number 28: makeUI eval=FALSE ################################################### ## yAnnoUI <- Genominator:::makeUIgenes(yAnno, gene.id = "ensembl_gene_id", transcript.id = "ensembl_transcript_id", verbose = TRUE) ## subset(yAnnoUI, ensembl_gene_id == "YAL005C") ## subset(yAnno, ensembl_gene_id == "YAL005C") ################################################### ### chunk number 29: makeUIsave eval=FALSE ################################################### ## save(yAnnoUI, file = "data/yAnnoUI.rda") ################################################### ### chunk number 30: makeUIload ################################################### load("data/yAnnoUI.rda") ################################################### ### chunk number 31: genecountsUI ################################################### geneCountsUI <- summarizeByAnnotation(eData, yAnnoUI, ignoreStrand = TRUE, meta.id = "ensembl_gene_id") ################################################### ### chunk number 32: gofplot1 ################################################### plot(regionGoodnessOfFit(geneCountsUI, groups = rep(c("mut", "wt"), times = c(2,2))), chisq = TRUE) ################################################### ### chunk number 33: gofplot2 ################################################### plot(regionGoodnessOfFit(geneCountsUI, groups = rep("all", 4)), chisq = TRUE) ################################################### ### chunk number 34: collapse ################################################### eData2 <- collapseExpData(eData, tablename = "collapsed", collapse = "sum", groups = rep(c("mut", "wt"), c(2,2)), verbose = TRUE) ################################################### ### chunk number 35: upper.quartile ################################################### notZero <- which(rowSums(geneCountsUI) != 0) upper.quartiles <- apply(geneCountsUI[notZero,], 2, function(x) quantile(x, 0.75)) uq.scaled <- upper.quartiles / sum(upper.quartiles) * sum(laneCounts) laneCounts uq.scaled sum(laneCounts) sum(uq.scaled) ################################################### ### chunk number 36: LR eval=FALSE ################################################### ## groups <- factor(rep(c("mut", "wt"), times = c(2,2))) ## pvalues <- apply(geneCountsUI[notZero,], 1, function(y) { ## fit <- glm(y ~ groups, family = poisson(), offset = log(uq.scaled)) ## fit0 <- glm(y ~ 1, family = poisson(), offset = log(uq.scaled)) ## anova(fit0, fit, test = "Chisq")[2,5] ## }) ################################################### ### chunk number 37: eval=FALSE ################################################### ## save(pvalues, file = "data/pvalues.rda") ################################################### ### chunk number 38: ################################################### load("data/pvalues.rda") ################################################### ### chunk number 39: multtest ################################################### library(multtest) adj <- mt.rawp2adjp(pvalues, proc = "BH") adj <- adj$adjp[order(adj$index),] rownames(adj) <- names(pvalues) ################################################### ### chunk number 40: plotting ################################################### annoFactory <- Genominator:::makeAnnoFactory.AnnoData(cbind(yAnnoUI, feature = "gene"), featureColumnName = "feature", groupColumnName = NULL, idColumnName = "ensembl_gene_id", dp = DisplayPars(plotId = TRUE, idColor = "black")) rp <- Genominator:::makeRegionPlotter( list("mut" = list(expData = eData2, what = "mut", fxs = Genominator:::makeConvolver(26), dp = DisplayPars(lwd = 0.2, color = "red")), "wt" = list(expData = eData2, what = "wt", fxs = Genominator:::makeConvolver(26), dp = DisplayPars(lwd = 0.3, color = "blue", type = "p"))), annoFactory = annoFactory) rp(1, 100000, 102000) ################################################### ### chunk number 41: sessionInfo ################################################### toLatex(sessionInfo())