## ----prelim-------------------------------------------------------------- suppressMessages(library(CNPBayes)) suppressMessages(library(SummarizedExperiment)) se <- readRDS(system.file("extdata", "simulated_se.rds", package="CNPBayes")) grl <- readRDS(system.file("extdata", "grl_deletions.rds", package="CNPBayes")) ## ----plot-cnvs----------------------------------------------------------- cnv.region <- consensusCNP(grl[1], max.width=5e6) i <- subjectHits(findOverlaps(cnv.region, rowRanges(se))) xlim <- c(min(start(se)), max(end(se))) par(las=1) plot(0, xlim=xlim, ylim=c(0, 26), xlab="Mb", ylab="sample index", type="n", xaxt="n") at <- pretty(xlim, n=10) axis(1, at=at, labels=round(at/1e6, 1), cex.axis=0.8) rect(start(grl), seq_along(grl)-0.2, end(grl), seq_along(grl)+0.2, col="gray", border="gray") ## ----median-------------------------------------------------------------- cnv.region <- consensusCNP(grl, max.width=5e6) i <- subjectHits(findOverlaps(cnv.region, rowRanges(se))) med.summary <- matrixStats::colMedians(assays(se)[["cn"]][i, ], na.rm=TRUE) par(las=1) plot(0, xlim=xlim, ylim=c(0, 26), xlab="Mb", ylab="sample index", type="n", xaxt="n") at <- pretty(xlim, n=10) axis(1, at=at, labels=round(at/1e6, 1), cex.axis=0.8) rect(start(grl), seq_along(grl)-0.2, end(grl), seq_along(grl)+0.2, col="gray", border="gray") abline(v=c(start(cnv.region), end(cnv.region)), lty=2) ## ----pca----------------------------------------------------------------- cnv.region2 <- reduce(unlist(grl)) i.pc <- subjectHits(findOverlaps(cnv.region2, rowRanges(se))) x <- oligoClasses::lrr(se)[i.pc, ] nas <- rowSums(is.na(x)) na.index <- which(nas > 0) x <- x[-na.index, , drop=FALSE] pc.summary <- prcomp(t(x))$x[, 1] meds.for.pc <- matrixStats::colMedians(x, na.rm=TRUE) if(cor(pc.summary, meds.for.pc) < 0) pc.summary <- -1*pc.summary par(las=1) plot(0, xlim=xlim, ylim=c(0, 26), xlab="Mb", ylab="sample index", type="n", xaxt="n") at <- pretty(xlim, n=10) axis(1, at=at, labels=round(at/1e6, 1), cex.axis=0.8) rect(start(grl), seq_along(grl)-0.2, end(grl), seq_along(grl)+0.2, col="gray", border="gray") abline(v=c(start(cnv.region2), end(cnv.region2)), lty=2) ## ----summary-plots------------------------------------------------------- par(mfrow=c(1,2), las=1) plot(med.summary, main="median summary of\nconsensus CNP", cex.main=0.7, pch=20) plot(pc.summary, main="PC summary of\nreduced CNV ranges", cex.main=0.7, pch=20) ## ----shapiro------------------------------------------------------------- shapiro.test(med.summary)