## ----setup, include=FALSE----------------------------------------------------- knitr::opts_chunk$set(echo = TRUE) knitr::opts_chunk$set(fig.width = 7) knitr::opts_chunk$set(fig.height = 5) knitr::opts_chunk$set(dpi = 72) # preload to avoid loading messages library(NanoMethViz) library(dplyr) exon_tibble <- get_exons_mm10() ## ----eval = FALSE------------------------------------------------------------- # if (!requireNamespace("BiocManager", quietly = TRUE)) # install.packages("BiocManager") # # BiocManager::install("NanoMethViz") ## ----------------------------------------------------------------------------- library(NanoMethViz) library(dplyr) ## ----------------------------------------------------------------------------- # methylation data stored in tabix file methy <- system.file(package = "NanoMethViz", "methy_subset.tsv.bgz") # tabix is just a special gzipped tab-separated-values file read.table(gzfile(methy), col.names = methy_col_names(), nrows = 6) ## ----------------------------------------------------------------------------- # helper function extracts exons from TxDb package exon_tibble <- get_exons_mm10() head(exon_tibble) ## ----------------------------------------------------------------------------- sample <- c( "B6Cast_Prom_1_bl6", "B6Cast_Prom_1_cast", "B6Cast_Prom_2_bl6", "B6Cast_Prom_2_cast", "B6Cast_Prom_3_bl6", "B6Cast_Prom_3_cast" ) group <- c("bl6", "cast", "bl6", "cast", "bl6", "cast") sample_anno <- data.frame(sample, group, stringsAsFactors = FALSE) sample_anno ## ----------------------------------------------------------------------------- nmr <- NanoMethResult(methy, sample_anno, exon_tibble) ## ----------------------------------------------------------------------------- plot_gene(nmr, "Peg3") ## ----------------------------------------------------------------------------- # loading saved results from previous bsseq analysis bsseq_dmr <- read.table( system.file(package = "NanoMethViz", "dmr_subset.tsv.gz"), sep = "\t", header = TRUE, stringsAsFactors = FALSE ) ## ----------------------------------------------------------------------------- plot_gene(nmr, "Peg3", anno_regions = bsseq_dmr) ## ----warning = FALSE---------------------------------------------------------- # warnings have been turned off in this vignette, but this will generally # generate many warnings as the smoothing for many reads will fail plot_gene(nmr, "Peg3", anno_regions = bsseq_dmr, spaghetti = TRUE) ## ----warning = FALSE---------------------------------------------------------- plot_gene(nmr, "Peg3", anno_regions = bsseq_dmr, heatmap = TRUE) ## ----------------------------------------------------------------------------- methy_calls <- system.file(package = "NanoMethViz", c("sample1_nanopolish.tsv.gz", "sample2_nanopolish.tsv.gz")) # have a look at the first 10 rows of methy_data methy_calls_example <- read.table( methy_calls[1], sep = "\t", header = TRUE, nrows = 6) methy_calls_example ## ----message=F---------------------------------------------------------------- # create a temporary file to store the converted data methy_tabix <- file.path(tempdir(), "methy_data.bgz") samples <- c("sample1", "sample2") # you should see messages when running this yourself create_tabix_file(methy_calls, methy_tabix, samples) # you don't need to do this with real data # we have to use gzfile to tell R that we have a gzip compressed file methy_data <- read.table( gzfile(methy_tabix), col.names = methy_col_names(), nrows = 6) methy_data ## ----------------------------------------------------------------------------- # construct with a ModBamFiles object as the methylation data mbr <- ModBamResult( methy = ModBamFiles( samples = "sample1", system.file(package = "NanoMethViz", "peg3.bam") ), samples = data.frame( sample = "sample1", group = 1 ), exons = exon_tibble ) # use in the same way as you would a NanoMethResult object plot_gene(mbr, "Peg3", heatmap = TRUE) ## ----message = FALSE---------------------------------------------------------- nmr <- load_example_nanomethresult() bss <- methy_to_bsseq(nmr) bss ## ----------------------------------------------------------------------------- gene_regions <- exons_to_genes(NanoMethViz::exons(nmr)) edger_mat <- bsseq_to_edger(bss, gene_regions) edger_mat ## ----------------------------------------------------------------------------- anno <- rtracklayer::import(system.file(package = "NanoMethViz", "c_elegans.gtf.gz")) head(anno) ## ----------------------------------------------------------------------------- anno <- anno %>% as.data.frame() %>% dplyr::rename( chr = "seqnames", symbol = "gene_name" ) %>% dplyr::select("gene_id", "chr", "strand", "start", "end", "transcript_id", "symbol") head(anno) ## ----message = FALSE---------------------------------------------------------- nmr <- load_example_nanomethresult() plot_gene(nmr, "Peg3") ## ----------------------------------------------------------------------------- new_exons <- NanoMethViz::exons(nmr) %>% exons_to_genes() %>% mutate(transcript_id = gene_id) NanoMethViz::exons(nmr) <- new_exons plot_gene(nmr, "Peg3") ## ----message = FALSE---------------------------------------------------------- # convert to bsseq bss <- methy_to_bsseq(nmr) bss ## ----------------------------------------------------------------------------- # create gene annotation from exon annotation gene_anno <- exons_to_genes(NanoMethViz::exons(nmr)) # create log-methylation-ratio matrix lmr <- bsseq_to_log_methy_ratio(bss, regions = gene_anno) ## ----------------------------------------------------------------------------- plot_mds(lmr) + ggtitle("MDS Plot") plot_pca(lmr) + ggtitle("PCA Plot") ## ----------------------------------------------------------------------------- new_labels <- gsub("B6Cast_Prom_", "", colnames(lmr)) new_labels <- gsub("(\\d)_(.*)", "\\2 \\1", new_labels) groups <- gsub(" \\d", "", new_labels) plot_mds(lmr, labels = new_labels, groups = groups) + ggtitle("MDS Plot") + scale_colour_brewer(palette = "Set1") ## ----------------------------------------------------------------------------- plot_mds(lmr, labels = NULL, groups = groups) + ggtitle("MDS Plot") + scale_colour_brewer(palette = "Set1") ## ----eval = FALSE------------------------------------------------------------- # options("NanoMethViz.site_filter" = 5) ## ----eval = FALSE------------------------------------------------------------- # options("NanoMethViz.highlight_col" = "red")