## ----setup, echo=FALSE, results="hide"---------------------------------------- knitr::opts_chunk$set(tidy=FALSE, cache=FALSE, dev="png", message=FALSE, error=FALSE, warning=FALSE) ## ----echo=FALSE--------------------------------------------------------------- knitr::include_graphics("images/SEESAW.png") ## ----------------------------------------------------------------------------- suppressPackageStartupMessages(library(ensembldb)) library(EnsDb.Hsapiens.v86) library(fishpond) edb <- EnsDb.Hsapiens.v86 t2g <- makeTx2Tss(edb) # GRanges object mcols(t2g)[,c("tx_id","group_id")] ## ----eval=FALSE--------------------------------------------------------------- # n <- ncol(y)/2 # rep1a1 <- assay(y, "infRep1")[,y$allele == "a1"] # rep1a2 <- assay(y, "infRep1")[,y$allele == "a2"] # mcols(y)$someInfo <- rowSums(abs(rep1a1 - rep1a2) < 1) < n # y <- y[ mcols(y)$someInfo, ] ## ----------------------------------------------------------------------------- suppressPackageStartupMessages(library(SummarizedExperiment)) ## ----------------------------------------------------------------------------- set.seed(1) y <- makeSimSwishData(allelic=TRUE) colData(y) levels(y$allele) # a1/a2 allelic fold changes ## ----echo=FALSE--------------------------------------------------------------- # hidden chunk to add ranges to the `se` tss <- t2g[!duplicated(t2g$group_id)] tss$symbol <- mapIds(edb, tss$gene_id, "SYMBOL", "GENEID") names(tss) <- paste0(tss$symbol, "-", tss$tss) mcols(tss) <- mcols(tss)[,c("tx_id","gene_id","tss","group_id","symbol")] # slow... #tx_id <- CharacterList(split(t2g$tx_id,t2g$group_id)) #tss$tx_id <- tx_id[names(tss)] tab <- table(tss$gene_id) tss$ntss <- as.numeric(tab[tss$gene_id]) tss <- tss[tss$ntss > 1 & tss$ntss < 5 & seqnames(tss) == "1"] tss <- tss[order(tss$gene_id),] tss <- tss[43:1042] # swap 2nd and 3rd isoform of first gene tss[2:3] <- tss[3:2] rowRanges(y) <- tss ## ----fig.dim=c(7,3.5)--------------------------------------------------------- y <- computeInfRV(y) # for posterior mean, variance gene <- rowRanges(y)$gene_id[1] idx <- mcols(y)$gene_id == gene plotAllelicHeatmap(y, idx=idx) ## ----------------------------------------------------------------------------- y <- labelKeep(y) y <- swish(y, x="allele", pair="sample") ## ----fig.dim=c(8,4)----------------------------------------------------------- dat <- data.frame(minusLogQ=-log10(mcols(y)$qvalue[idx]), row.names=rownames(y)[idx]) plotAllelicHeatmap(y, idx=idx, annotation_row=dat) ## ----fig.dim=c(5,5)----------------------------------------------------------- par(mfrow=c(2,1), mar=c(1,4.1,2,2)) plotInfReps(y, idx=1, x="allele", cov="sample", xaxis=FALSE, xlab="") plotInfReps(y, idx=2, x="allele", cov="sample", xaxis=FALSE, xlab="") ## ----fig.dim=c(8,7)----------------------------------------------------------- gene <- rowRanges(y)$gene_id[1] plotAllelicGene(y, gene, edb) ## ----fig.dim=c(8,7)----------------------------------------------------------- plotAllelicGene(y, symbol="ADSS", db=edb) ## ----fig.dim=c(8,7)----------------------------------------------------------- plotAllelicGene(y, gene, edb, transcriptAnnotation="transcript", labels=list(a2="maternal",a1="paternal")) ## ----fig.dim=c(5,4)----------------------------------------------------------- y$allele_new <- y$allele # note a2 is non-effect, a1 is effect: levels(y$allele) # replace a2 then a1: levels(y$allele_new) <- c("maternal","paternal") plotInfReps(y, idx=1, x="allele_new", legend=TRUE, legendPos="bottom") ## ----------------------------------------------------------------------------- set.seed(1) y <- makeSimSwishData(dynamic=TRUE) colData(y) ## ----echo=FALSE--------------------------------------------------------------- rowRanges(y) <- tss ## ----------------------------------------------------------------------------- y <- labelKeep(y) y <- swish(y, x="allele", pair="sample", cov="time", cor="pearson") ## ----------------------------------------------------------------------------- mcols(y)[1:2,c("stat","qvalue")] ## ----------------------------------------------------------------------------- y <- computeInfRV(y) ## ----fig.dim=c(7,7)----------------------------------------------------------- par(mfrow=c(2,1), mar=c(2.5,4,2,2)) plotInfReps(y, idx=1, x="time", cov="allele", shiftX=.01, xaxis=FALSE, xlab="", main="") par(mar=c(4.5,4,0,2)) plotInfReps(y, idx=2, x="time", cov="allele", shiftX=.01, main="") ## ----fig.dim=c(7,5)----------------------------------------------------------- plotInfReps(y, idx=1, x="time", cov="allele", shiftX=.01) dat <- data.frame( time = y$time[1:10], a2 = assay(y, "mean")[1,y$allele=="a2"], a1 = assay(y, "mean")[1,y$allele=="a1"]) lines(lowess(dat[,c(1,2)]), col="dodgerblue") lines(lowess(dat[,c(1,3)]), col="goldenrod4") ## ----fig.dim=c(8,4)----------------------------------------------------------- idx <- c(1:4) row_dat <- data.frame(minusLogQ=-log10(mcols(y)$qvalue[idx]), row.names=rownames(y)[idx]) col_dat <- data.frame(time=y$time[1:10], row.names=paste0("s",1:10)) plotAllelicHeatmap(y, idx=idx, annotation_row=row_dat, annotation_col=col_dat) ## ----------------------------------------------------------------------------- y$time_bins <- cut(y$time,breaks=c(0,.25,.75,1), include.lowest=TRUE, labels=FALSE) y$time_bins <- paste0("time-",y$time_bins) table(y$time_bins[ y$allele == "a2" ]) ## ----fig.dim=c(8,7)----------------------------------------------------------- gene <- rowRanges(y)$gene_id[1] plotAllelicGene(y, gene, edb, cov="time_bins", qvalue=FALSE, log2FC=FALSE) ## ----fig.dim=c(8,7)----------------------------------------------------------- plotAllelicGene(y, gene, edb, cov="time_bins", covFacetIsoform=TRUE, qvalue=FALSE, log2FC=FALSE) ## ----------------------------------------------------------------------------- sessionInfo()