## ---- include = FALSE----------------------------------------------------------------------------- knitr::opts_chunk$set( collapse = TRUE, comment = "#>", width = 100 ) options(width=100) # devtools::load_all(".") # delete later ## ---- fig.width=10, fig.height=4, out.width="100%", out.height="100%"----------------------------- library(ramr) data(ramr) head(ramr.samples) ramr.data[1:10,ramr.samples[1:3]] plotAMR(ramr.data, ramr.samples, ramr.tp.unique[1]) plotAMR(ramr.data, ramr.samples, ramr.tp.nonunique[c(1,6,11)]) ## ---- fig.width=10, fig.height=6, out.width="100%", out.height="100%"----------------------------- # set the seed if reproducible results required set.seed(999) # unique random AMRs amrs.unique <- simulateAMR(ramr.data, nsamples=25, regions.per.sample=2, min.cpgs=5, merge.window=1000, dbeta=0.2) # non-unique AMRs outside of regions with unique AMRs amrs.nonunique <- simulateAMR(ramr.data, nsamples=4, exclude.ranges=amrs.unique, samples.per.region=2, min.cpgs=5, merge.window=1000) # random noise outside of AMR regions noise <- simulateAMR(ramr.data, nsamples=25, regions.per.sample=20, exclude.ranges=c(amrs.unique, amrs.nonunique), min.cpgs=1, max.cpgs=1, merge.window=1, dbeta=0.5) # "smooth" methylation data without AMRs (negative control) smooth.data <- simulateData(ramr.data, nsamples=25, cores=2) # methylation data with AMRs and noise noisy.data <- simulateData(ramr.data, nsamples=25, amr.ranges=c(amrs.unique, amrs.nonunique, noise), cores=2) # that's how regions look like library(gridExtra) do.call("grid.arrange", c(plotAMR(noisy.data, amr.ranges=amrs.unique[1:4]), ncol=2)) do.call("grid.arrange", c(plotAMR(noisy.data, amr.ranges=sort(amrs.nonunique)[1:8]), ncol=2)) do.call("grid.arrange", c(plotAMR(noisy.data, amr.ranges=noise[1:4]), ncol=2)) # can we find them? system.time(found <- getAMR(noisy.data, ramr.method="beta", min.cpgs=5, merge.window=1000, qval.cutoff=1e-2, cores=2)) # all possible regions all.ranges <- getUniverse(noisy.data, min.cpgs=5, merge.window=1000) # true positives tp <- sum(found %over% c(amrs.unique, amrs.nonunique)) # false positives fp <- sum(found %outside% c(amrs.unique, amrs.nonunique)) # true negatives tn <- length(all.ranges %outside% c(amrs.unique, amrs.nonunique)) # false negatives fn <- sum(c(amrs.unique, amrs.nonunique) %outside% found) # accuracy, MCC acc <- (tp+tn) / (tp+tn+fp+fn) mcc <- (tp*tn - fp*fn) / (sqrt(tp+fp)*sqrt(tp+fn)*sqrt(tn+fp)*sqrt(tn+fn)) setNames(c(tp, fp, tn, fn), c("TP", "FP", "TN", "FN")) setNames(c(acc, mcc), c("accuracy", "MCC")) ## ---- fig.width=10, fig.height=6, out.width="100%", out.height="100%"----------------------------- # identify AMRs amrs <- getAMR(ramr.data, ramr.samples, ramr.method="beta", min.cpgs=5, merge.window=1000, qval.cutoff=1e-3, cores=2) # inspect sort(amrs) do.call("grid.arrange", c(plotAMR(ramr.data, ramr.samples, amrs[1:10]), ncol=2)) ## ------------------------------------------------------------------------------------------------- # annotating AMRs using R library annotatr library(annotatr) annotation.types <- c("hg19_cpg_inter", "hg19_cpg_islands", "hg19_cpg_shores", "hg19_cpg_shelves", "hg19_genes_intergenic", "hg19_genes_promoters", "hg19_genes_5UTRs", "hg19_genes_firstexons", "hg19_genes_3UTRs") annotations <- build_annotations(genome='hg19', annotations=annotation.types) amrs.annots <- annotate_regions(regions=amrs, annotations=annotations, ignore.strand=TRUE, quiet=FALSE) summarize_annotations(annotated_regions=amrs.annots, quiet=FALSE) ## ----session-------------------------------------------------------------------------------------- sessionInfo()