## ---- include = FALSE-----------------------------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  width = 100
)
options(width=100)
# devtools::load_all(".")   # delete later

## -------------------------------------------------------------------------------------------------
library(epialleleR)

capture.bam <- system.file("extdata", "capture.bam", package="epialleleR")
bam.data    <- preprocessBam(capture.bam)

## -------------------------------------------------------------------------------------------------
# data.table::data.table object for
# CpG VEF report
cg.vef.report <- generateCytosineReport(bam.data)
head(cg.vef.report[order(meth+unmeth, decreasing=TRUE)])

# CpG cytosine report
cg.report <- generateCytosineReport(bam.data, threshold.reads=FALSE)
head(cg.report[order(meth+unmeth, decreasing=TRUE)])

# CX cytosine report
cx.report <- generateCytosineReport(bam.data, threshold.reads=FALSE,
                                    report.context="CX")
head(cx.report[order(meth+unmeth, decreasing=TRUE)])

## -------------------------------------------------------------------------------------------------
# report for amplicon-based data
# matching is done by exact start or end positions plus/minus tolerance
amplicon.report <- generateAmpliconReport(
  bam=system.file("extdata", "amplicon010meth.bam", package="epialleleR"),
  bed=system.file("extdata", "amplicon.bed", package="epialleleR")
)
amplicon.report

# report for capture-based data
# matching is done by overlap
capture.report <- generateCaptureReport(
  bam=system.file("extdata", "capture.bam", package="epialleleR"),
  bed=system.file("extdata", "capture.bed", package="epialleleR")
)
head(capture.report)

# generateBedReport is a generic function for BED-guided reports
bed.report <- generateBedReport(
  bam=system.file("extdata", "capture.bam", package="epialleleR"),
  bed=system.file("extdata", "capture.bed", package="epialleleR"),
  bed.type="capture"
)
identical(capture.report, bed.report)

## ---- fig.width=10, fig.height=6, out.width="100%", out.height="100%"-----------------------------
# First, let's extract base methylation information for sequencing reads
# of 1:9 mix of methylated and non-methylated control DNA
patterns <- extractPatterns(
  bam=system.file("extdata", "amplicon010meth.bam", package="epialleleR"),
  bed=as("chr17:43125200-43125600","GRanges")
)

# that many read pairs overlap genomic region of interest
nrow(patterns)

# these are positions of bases within relevant read pairs
base.positions <- grep("^[0-9]+$", colnames(patterns), value=TRUE)

# let's make a summary table with counts for every pattern
patterns.summary <- patterns[, c(lapply(.SD, unique), .N),
                             by=.(pattern, beta), .SDcols=base.positions]

# that many unique methylation patterns were found
nrow(patterns.summary)

# let's melt and plot patterns
plot.data <- data.table::melt.data.table(patterns.summary,
  measure.vars=base.positions, variable.name="pos", value.name="base")

# categorical positions, all patterns sorted by beta, with counts on the right
if (require("ggplot2", quietly=TRUE) & require("ggstance", quietly=TRUE)) {
  ggplot(na.omit(plot.data),
         aes(x=pos, y=reorder(pattern,beta),
             group=pattern, color=factor(base))) +
    geom_line(color="grey", position=position_dodgev(height=0.5)) +
    geom_point(position=position_dodgev(height=0.5)) +
    scale_colour_grey(start=0.8, end=0) +
    theme_light() +
    theme(axis.text.x=element_text(angle=60, hjust=1, vjust=1),
          axis.text.y=element_blank()) +
    labs(x="position", y="pattern", title="epialleles", color="base") +
    scale_x_discrete(expand=c(0.05,0)) +
    geom_text(inherit.aes=FALSE, data=patterns.summary,
              mapping=aes(x="count", y=pattern, label=N), size=3)
}

# continuous positions, nonunique patterns according to their counts
if (require("ggplot2", quietly=TRUE) & require("ggstance", quietly=TRUE)) {
  ggplot(na.omit(plot.data)[N>1],
         aes(x=as.numeric(as.character(pos)), y=factor(N),
             group=pattern, color=factor(base))) +
    geom_line(color="grey", position=position_dodgev(height=0.5)) +
    geom_point(position=position_dodgev(height=0.5)) +
    scale_colour_grey(start=0.8, end=0) +
    theme_light() +
    labs(x="position", y="count", title="epialleles", color="base")
}

# now let's explore methylation patterns in RAD51C gene promoter using
# methylation capture data
capture.patterns <- extractPatterns(
  bam=system.file("extdata", "capture.bam", package="epialleleR"),
  bed=as("chr17:58691673-58693108", "GRanges"),
  verbose=FALSE
)
capture.positions <- grep("^[0-9]+$", colnames(capture.patterns), value=TRUE)
capture.summary <-
  capture.patterns[, c(lapply(.SD, unique), .N),
                   by=.(pattern, beta), .SDcols=capture.positions]
capture.data <- data.table::melt.data.table(capture.summary,
  measure.vars=capture.positions, variable.name="pos", value.name="base")
if (require("ggplot2", quietly=TRUE) & require("ggstance", quietly=TRUE)) {
  ggplot(na.omit(capture.data),
         aes(x=as.numeric(as.character(pos)), y=factor(N),
             group=pattern, color=factor(base))) +
    geom_line(color="grey", position=position_dodgev(height=0.9)) +
    geom_point(position=position_dodgev(height=0.9)) +
    scale_colour_grey(start=0.8, end=0) +
    theme_light() +
    labs(x="position", y="count", title="RAD51C", color="base")
}

## ---- fig.width=10, fig.height=6, out.width="100%", out.height="100%"-----------------------------
# First, let's visualise eCDFs for within- and out-of-context beta values
# for all four amplicons and unmatched reads. Note that within-the-context eCDF
# of 0.5 is very close to the expected 1-VEF value (0.1) for all amplicons
# produced from this 1:9 mix of methylated and non-methylated control DNA

# let's compute eCDF
amplicon.ecdfs <- generateBedEcdf(
  bam=system.file("extdata", "amplicon010meth.bam", package="epialleleR"),
  bed=system.file("extdata", "amplicon.bed", package="epialleleR"),
  bed.rows=NULL
)

# there are 5 items in amplicon.ecdfs, let's plot all of them
par(mfrow=c(1,length(amplicon.ecdfs)))

# cycle through items
for (x in 1:length(amplicon.ecdfs)) {
  # four of them have names corresponding to genomic regions of amplicon.bed
  # fifth - NA for all the reads that don't match to any of those regions
  main <- if (is.na(names(amplicon.ecdfs[x]))) "unmatched"
          else names(amplicon.ecdfs[x])
  
  # plotting eCDF for within-the-context per-read beta values (in red)
  plot(amplicon.ecdfs[[x]]$context, col="red", verticals=TRUE, do.points=FALSE,
       xlim=c(0,1), xlab="per-read beta value", ylab="cumulative density",
       main=main)
  
  # adding eCDF for out-of-context per-read beta values (in blue)
  plot(amplicon.ecdfs[[x]]$out.of.context, add=TRUE, col="blue",
       verticals=TRUE, do.points=FALSE)
}


# Second, let's compare eCDFs for within-the-context beta values for only one
# amplicon but all three sequenced samples: pure non-methylated DNA, 1:9 mix of
# methylated and non-methylated DNA, and fully methylated DNA

# our files
bam.files <- c("amplicon000meth.bam", "amplicon010meth.bam",
               "amplicon100meth.bam")

# let's plot all of them
par(mfrow=c(1,length(bam.files)))

# cycle through items
for (f in bam.files) {
  # let's compute eCDF
  amplicon.ecdfs <- generateBedEcdf(
    bam=system.file("extdata", f, package="epialleleR"),
    bed=system.file("extdata", "amplicon.bed", package="epialleleR"),
    # only the second amplicon
    bed.rows=2, verbose=FALSE
  )
  
  # plotting eCDF for within-the-context per-read beta values (in red)
  plot(amplicon.ecdfs[[1]]$context, col="red", verticals=TRUE, do.points=FALSE,
       xlim=c(0,1), xlab="per-read beta value", ylab="cumulative density",
       main=f)
  
   # adding eCDF for out-of-context per-read beta values (in blue)
  plot(amplicon.ecdfs[[1]]$out.of.context, add=TRUE, col="blue",
       verticals=TRUE, do.points=FALSE)
}

## -------------------------------------------------------------------------------------------------
# VCF report
vcf.report <- generateVcfReport(
  bam=system.file("extdata", "amplicon010meth.bam", package="epialleleR"),
  bed=system.file("extdata", "amplicon.bed", package="epialleleR"),
  vcf=system.file("extdata", "amplicon.vcf.gz", package="epialleleR"),
  # thresholds on alignment and base quality
  min.mapq=30, min.baseq=13,
  # when VCF seqlevels are different from BED and BAM it is possible
  # to convert them internally
  vcf.style="NCBI"
)

# NA values are shown for the C->T variants on the "+" and G->A on the "-"
# strands, because bisulfite conversion makes their counting impossible
head(vcf.report)

# let's sort the report by increasing Fisher's exact test's p-values.
# the p-values are given separately for reads that map to the "+"
head(vcf.report[order(`FEp-`, na.last=TRUE)])

# and to the "-" strand
head(vcf.report[order(`FEp+`, na.last=TRUE)])

# and finally, let's plot methylation patterns overlapping one of the most
# covered SNPs in the methylation capture test data set - rs573296191
# (chr17:61864584) in BRIP1 gene
brip1.patterns <- extractPatterns(
  bam=system.file("extdata", "capture.bam", package="epialleleR"),
  bed=as("chr17:61864583-61864585", "GRanges"),
  highlight.positions=61864584,
  verbose=FALSE
)
brip1.positions <- grep("^[0-9]+$", colnames(brip1.patterns), value=TRUE)
brip1.summary <-
  brip1.patterns[, c(lapply(.SD, unique), .N),
                 by=.(pattern, beta), .SDcols=brip1.positions]
brip1.data <- data.table::melt.data.table(brip1.summary,
  measure.vars=brip1.positions, variable.name="pos", value.name="base")
if (require("ggplot2", quietly=TRUE) & require("ggstance", quietly=TRUE)) {
  ggplot(na.omit(brip1.data),
         aes(x=as.numeric(as.character(pos)), y=factor(N),
             group=pattern, color=factor(base))) +
    geom_line(color="grey", position=position_dodgev(height=0.9)) +
    geom_point(position=position_dodgev(height=0.9)) +
    scale_colour_manual(values=c("blue","red","grey80","black")) +
    theme_light() +
    labs(x="position", y="count", title="BRIP1", color="base")
}

## ----session--------------------------------------------------------------------------------------
sessionInfo()

