## ----setup, include=FALSE----------------------------------------------------- knitr::opts_chunk$set(cache = 0, fig.width = 6, fig.height = 7) ## ----style, echo = FALSE, results = 'asis'------------------------------------ BiocStyle::markdown() ## ----eval=FALSE--------------------------------------------------------------- # if (!requireNamespace("BiocManager", quietly = TRUE)) # install.packages("BiocManager") # # BiocManager::install("DAMEfinder") ## ----------------------------------------------------------------------------- suppressPackageStartupMessages({ library(DAMEfinder) library(SummarizedExperiment) library(GenomicRanges) library(BSgenome.Hsapiens.UCSC.hg19) }) bam_files <- c(system.file("extdata", "NORM1_chr19_trim.bam", package = "DAMEfinder"), system.file("extdata", "CRC1_chr19_trim.bam", package = "DAMEfinder")) vcf_files <- c(system.file("extdata", "NORM1.chr19.trim.vcf", package = "DAMEfinder"), system.file("extdata", "CRC1.chr19.trim.vcf", package = "DAMEfinder")) sample_names <- c("NORM1", "CRC1") #Use another reference file for demonstration, and fix the seqnames genome <- BSgenome.Hsapiens.UCSC.hg19 seqnames(genome) <- gsub("chr","",seqnames(genome)) reference_file <- DNAStringSet(genome[[19]], use.names = TRUE) names(reference_file) <- 19 #Extract reads and extract methylation according to allele snp.list <- extract_bams(bam_files, vcf_files, sample_names, reference_file, coverage = 2) #CpG sites for first SNP in VCF file from sample NORM1 snp.list$NORM1[[1]] #CpG sites for first SNP in VCF file from sample CRC1 snp.list$CRC1[[1]] ## ----------------------------------------------------------------------------- derASM <- calc_derivedasm(snp.list) derASM assays(derASM) ## ----------------------------------------------------------------------------- x <- assay(derASM, "der.ASM") head(x) ## ----------------------------------------------------------------------------- data(extractbams_output) #The data loaded is an output from `split_bams()`, therefore we run #`calc_derivedasm` to get the SummarizedExperiment derASM <- calc_derivedasm(extractbams_output, cores = 1, verbose = FALSE) #We remove all CpG sites with any NA values, but not 0s filt <- rowSums(!is.na(assay(derASM, "der.ASM"))) == 8 derASM <- derASM[filt,] #set the design matrix grp <- factor(c(rep("CRC",4),rep("NORM",4)), levels = c("NORM", "CRC")) mod <- model.matrix(~grp) mod #Run default dames <- find_dames(derASM, mod) head(dames) #Run empirical method dames <- find_dames(derASM, mod, pvalAssign = "empirical") head(dames) ## ----------------------------------------------------------------------------- tuple_files <- c(system.file("extdata", "NORM1_chr19.qs.CG.2.tsv.gz", package = "DAMEfinder"), system.file("extdata", "CRC1_chr19.qs.CG.2.tsv.gz", package = "DAMEfinder")) sample_names <- c("NORM1", "CRC1") tuple_list <- read_tuples(tuple_files, sample_names) head(tuple_list$NORM1) ## ----------------------------------------------------------------------------- ASM_mat <- calc_asm(tuple_list) ASM_mat ## ----------------------------------------------------------------------------- #load package data data(readtuples_output) #run calc_asm and filter object ASMscore <- calc_asm(readtuples_output) filt <- rowSums(!is.na(assay(ASMscore, "asm"))) == 5 #filt to avoid warnings ASMscore <- ASMscore[filt,] #make design matrix (or specify a contrast) grp <- factor(c(rep("CRC",3),rep("NORM",2)), levels = c("NORM", "CRC")) mod <- model.matrix(~grp) #run default and increase maxGap to get longer, more sparse regions dames <- find_dames(ASMscore, mod, maxGap = 300) head(dames) #run alternative mode dames <- find_dames(ASMscore, mod, maxGap = 300, pvalAssign = "empirical") head(dames) ## ----dametrack---------------------------------------------------------------- #Here I will use the tuple-ASM SummExp colData(ASMscore)$group <- grp colData(ASMscore)$samples <- colnames(ASMscore) #Set a DAME as a GRanges. I choose a one from the tables we obtained above dame <- GRanges(19,IRanges(323736,324622)) dame_track(dame = dame, ASM = ASMscore) ## ----dt2---------------------------------------------------------------------- dame_track(dame = dame, ASM = ASMscore, window = 2) ## ----dt3---------------------------------------------------------------------- dame <- GRanges(19,IRanges(387966,387983)) grp <- factor(c(rep("CRC",4),rep("NORM",4)), levels = c("NORM", "CRC")) colData(derASM)$group <- grp dame_track(dame = dame, derASM = derASM) ## ----dt4---------------------------------------------------------------------- dame_track(dame = dame, derASM = derASM, plotSNP = TRUE) ## ----dt5---------------------------------------------------------------------- dame_track(dame = dame, derASM = derASM, ASM = ASMscore) ## ----dt6---------------------------------------------------------------------- dame_track_mean(dame = dame, derASM = derASM, ASM = ASMscore) ## ----fig1--------------------------------------------------------------------- #put SNP in GRanges (you can find the SNP with the dame_track function) snp <- GRanges(19, IRanges(267039, width = 1)) #always set the width if your #GRanges has 1 site snp bam.file <- system.file("extdata", "CRC1_chr19_trim.bam", package = "DAMEfinder") vcf.file <- system.file("extdata", "CRC1.chr19.trim.vcf", package = "DAMEfinder") methyl_circle_plot(snp = snp, vcfFile = vcf.file, bamFile = bam.file, refFile = reference_file) ## ----fig2--------------------------------------------------------------------- cpgsite <- GRanges(19, IRanges(266998, width = 1)) methyl_circle_plot(snp = snp, vcfFile = vcf.file, bamFile = bam.file, refFile = reference_file, cpgsite = cpgsite) ## ----fig3--------------------------------------------------------------------- cpgsite <- GRanges(19, IRanges(266998, width = 1)) methyl_circle_plotCpG(cpgsite = cpgsite, bamFile = bam.file, refFile = reference_file) ## ----fig4--------------------------------------------------------------------- #a random region dame <- GRanges(19, IRanges(266998,267100)) methyl_circle_plot(snp = snp, vcfFile = vcf.file, bamFile = bam.file, refFile = reference_file, dame = dame) ## ----fig5--------------------------------------------------------------------- grp <- factor(c(rep("CRC",3),rep("NORM",2)), levels = c("NORM", "CRC")) methyl_MDS_plot(ASMscore, group = grp) ## ----------------------------------------------------------------------------- utils::sessionInfo()