---
title: "QC Report"
date: "`r Sys.Date()`"
output:
  html_document:
    toc: true
    theme: united
params:
  threads: 1
  input_bgz: NULL
  input_vcf: NULL
  input_bam: NULL
  reference: NULL
  output_dir: "."
  region: NULL # By default there's no region
  tabix: "tabix" # Assume tabix is available in PATH
  bcftools: "bcftools" # Assume bcftools is available in PATH
  rastair: "rastair" # Assume rastair is available in PATH
  include_flags: 3
  exclude_flags: 3852
  read_len: NULL
  plot_vbias: TRUE
  plot_gc: TRUE
  is_WGBS: FALSE
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```
```{r message=FALSE, warning=FALSE, include=FALSE}
suppressPackageStartupMessages({
  library(data.table)
  library(ggplot2)
})

parse_positions <- function(mod_pos, unmod_pos, snp_pos, moddn_pos, unmoddn_pos) {
  m_pos <- as.integer(strsplit(as.character(mod_pos), ",", fixed=TRUE)[[1]])
  u_pos <- as.integer(strsplit(as.character(unmod_pos), ",", fixed=TRUE)[[1]])
  s_pos <- as.integer(strsplit(as.character(snp_pos), ",", fixed=TRUE)[[1]])
  dm_pos <- as.integer(strsplit(as.character(moddn_pos), ",", fixed=TRUE)[[1]])
  du_pos <- as.integer(strsplit(as.character(unmoddn_pos), ",", fixed=TRUE)[[1]])

  nm <- ifelse(any(is.na(m_pos)), 0, length(m_pos))
  nu <- ifelse(any(is.na(u_pos)), 0, length(u_pos))
  ns <- ifelse(any(is.na(s_pos)), 0, length(s_pos))
  ndm <- ifelse(any(is.na(dm_pos)), 0, length(dm_pos))
  ndu <- ifelse(any(is.na(du_pos)), 0, length(du_pos))
  return(data.table(pos_in_read=na.omit(c(m_pos, u_pos, s_pos, dm_pos, du_pos)), type=c(rep("M", nm), rep("U", nu), rep("S", ns), rep("m", ndm), rep("u", ndu))))
}


# For a subset of observations of read orientation on one strand,
# calculate the suggested cutoff points
calculate_cutoff <- function(mbias_table) {
  # Ensure we process things in orders
  mbias_table <- mbias_table[order(mbias_table$pos),]
  # get the mean beta of the middle 60% of the data
  n_obs <- nrow(mbias_table)
  if (n_obs < 3) {
    stop("Not enough data points to calculate cutoffs")
  }
  upper_lim <- min(ceiling(n_obs * 0.8), n_obs-1)
  lower_lim <- max(floor(n_obs * 0.2), 2)
  mean_middle <- mean(mbias_table[lower_lim:upper_lim,]$beta)
  upper_ci_max <- max(mbias_table[lower_lim:upper_lim,]$upper_ci)
  lower_ci_min <- min(mbias_table[lower_lim:upper_lim,]$lower_ci)

  # Search from the inside out for the first point left/right that is
  # either outside the 95% CI bounds of the middle, or the mean of the
  # middle is outside the CI of the point, or it's more than 0.05 away
  # from the mean of the middle
  check_match <- function(row, mean_middle, lower_ci_min, upper_ci_max) {
    if (is.na(mean_middle) || is.null(mean_middle) || is.na(row$beta) || is.null(row$beta) || is.na(lower_ci_min) || is.null(lower_ci_min) || is.na(upper_ci_max) || is.null(upper_ci_max)) {
      return(FALSE)
    }

    if (!row$beta %between% c(lower_ci_min, upper_ci_max)
        && abs(row$beta - mean_middle) > 0.05
        && !mean_middle %between% c(row$lower_ci, row$upper_ci)) {
      return(TRUE)
    }

    return(FALSE)
  }

  left_cutoff = 0
  right_cutoff = 0

  middle = ceiling(n_obs/2)
  for (pos in middle:1) {
    if (check_match(mbias_table[pos,], mean_middle, lower_ci_min, upper_ci_max)) {
      left_cutoff = pos
      break
    }
  }
  for (pos in 1:(n_obs-middle)) {
    if (check_match(mbias_table[middle+pos,], mean_middle, lower_ci_min, upper_ci_max)) {
      right_cutoff = n_obs - (middle+pos) + 1
      break
    }
  }
  return(data.frame(left=as.double(left_cutoff), right=as.double(right_cutoff)))
}
```

```{r transform_data, include=FALSE}
# if no bgz given, but a bam file is provided, generate bed file on the fly
if ((is.na(params$input_bgz) || is.null(params$input_bgz) || nchar(params$input_bgz)==0 || !file.exists(params$input_bgz)) && !is.na(params$input_bam) && !is.null(params$input_bam) && nchar(params$input_bam)>0 && file.exists(params$input_bam)) {
  if (is.na(params$reference) || is.null(params$reference) || nchar(params$reference)==0 || !file.exists(params$reference)) {
    stop("Cannot run from bam file without a reference")
  }
  
  if (!is.na(params$region) && nchar(params$region) > 0) {
    command=paste(params$rastair, "per-read --bed - -r", params$reference, "-l", params$region, "-f", params$include_flags, "-F", params$exclude_flags, "-@", params$threads, params$input_bam)
  } else {
    command=paste(params$rastair, "per-read --bed - -r", params$reference, "-f", params$include_flags, "-F", params$exclude_flags, "-@", params$threads, params$input_bam)
  }
} else if (!(is.na(params$input_bgz) || is.null(params$input_bgz) || nchar(params$input_bgz)==0 || !file.exists(params$input_bgz))) {
  if (!is.na(params$region) && nchar(params$region) > 0) {
    command=paste(params$tabix, "-h", params$input_bgz, params$region)
  } else {
    command=paste("zcat", params$input_bgz)
  } 
}  else {
  stop("Neither per-read bgz nor bam file given, cannot proceed")
}

data <- fread(cmd=command, header=TRUE, stringsAsFactors = FALSE)
setnames(data, 1, "chr")
data[, chr := as.factor(chr)]
data[, orientation := as.factor(orientation)]

if (!is.null(params$is_WGBS) && params$is_WGBS == TRUE) {
  data[, num_mod := num_cpg - num_mod] # flip mod/unmod
  setnames(data, c("mod_cpgs", "unmod_cpgs", "mod_denovos", "unmod_denovos"), c("unmod_cpgs", "mod_cpgs", "unmod_denovos", "mod_denovos"))
}

# remove reads without CpGs
data <- data[num_cpg > 0]

# filter by flags
data <- data[bitwAnd(flag, params$include_flags) == params$include_flags]
data <- data[bitwAnd(flag, params$exclude_flags) == 0]

if (is.null(params$read_len) || is.na(params$read_len) || !is.numeric(params$read_len)) {
  read_len <- max(data$read_length, na.rm=TRUE)
} else {
  read_len = params$read_len
}

data <- data[, .N, by=.(chr, start, end, read_id, mapq, orientation, insert_size, read_length, flag, mod_cpgs, unmod_cpgs, snp_cpgs, mod_denovos, unmod_denovos)] # unify duplicate reads

# Break comma-separated pos_in_read into a list
data <- data[, parse_positions(mod_cpgs, unmod_cpgs, snp_cpgs, mod_denovos, unmod_denovos), by=.(chr, start, end, read_id, mapq, orientation, insert_size, read_length, flag)]

mbias_table <- data[, .(mod=sum(type %in% c("M", "m")), unmod=sum(type %in% c("U","u"))), by=.(chr, pos=pos_in_read, read_pair=as.factor(fifelse(bitwAnd(flag, 64)>0, "First", "Second")), strand=as.factor(fifelse(bitwAnd(flag, 80)==80 | bitwAnd(flag, 160)==160, "OT", "OB")))]

# Flip the OB reads so that mbias is always reported 5' -> 3'
mbias_table[strand=="OT" & read_pair == "Second" | strand=="OB" & read_pair == "First", pos := read_len-pos]
mbias_table <- mbias_table[, .(beta=mod/(mod+unmod), lower_ci=ifelse(mod+unmod>0, binom.test(mod, mod+unmod, conf.level = 0.999)$conf.int[1], 0.0), upper_ci=ifelse(mod+unmod>0, binom.test(mod, mod+unmod, conf.level = 0.999)$conf.int[2], 0.0)), by=.(chr, pos, read_pair, strand)]
mbias_table <- mbias_table[!is.na(beta)]
setkey(mbias_table, chr, pos)
```

## M-bias plots

```{r plot_mbias, echo=FALSE, fig.height=5, fig.width=7}
for (chrom in unique(mbias_table$chr)) {
  cutoffs <- mbias_table[chr==chrom, calculate_cutoff(as.data.frame(.SD)), by=.(chr, read_pair, strand)]
  print(cutoffs)
  #write cutoffs to file if requested
  if (!is.na(params$output_dir) && !is.null(params$output_dir)) {
    out_table <- cutoffs[, .(l_ot = .SD[read_pair=="First", left],
                             r_ot = .SD[read_pair=="First", right],
                             l_ob = .SD[read_pair=="Second", left],
                             r_ob = .SD[read_pair=="Second", right]),
                         by=strand]
    setDF(out_table)
    write.table(out_table[,2:5], file.path(params$output_dir, paste0(chrom, "_cutoffs.txt")), row.names = FALSE, col.names = FALSE, quote=FALSE, sep=",")
  }

  g <- ggplot(mbias_table[chr==chrom], aes(x=pos, y=beta, color=read_pair)) +
        geom_line() +
        geom_ribbon(aes(ymin=lower_ci, ymax=upper_ci, fill=read_pair), alpha=0.2) +
        labs(color = "Read") +
        xlab("Position in read (5'->3')") +
        scale_color_discrete(labels=c('First in pair', 'Second in pair')) +
        scale_fill_discrete(guide="none") +
        theme_light()
  if(!is.null(cutoffs)) {
    g = g + geom_vline(aes(xintercept=left, color=read_pair), data=cutoffs)
    g = g + geom_vline(aes(xintercept=read_len-right, color=read_pair), data=cutoffs)
  }
  g = g + facet_wrap(vars(strand), nrow = 2, as.table = FALSE)
  g = g + ggtitle(chrom)

  print(g)
}
```

## V-bias plots

```{r transform_vbias, include=FALSE}
if (params$plot_vbias) {
  suppressPackageStartupMessages({
    library(Rsamtools)
    library(GenomicRanges)
    library(Biostrings)
    library(ggside)
  })
  
  # List all CpGs in the reference sequence for a particular region
  cpg_table_from_fasta <- function(region,
                                   fasta,
                                   build_index_if_missing = FALSE) {
    
    elements <- unlist(strsplit(region, "[:\\-]"))
    if (length(elements == 1)) {
      chr = elements[1]
      start_pos = 1
      end_pos = NULL
    } else if (length(elements) == 3) {
      chr = elements[1]
      start_pos = elements[2]
      end_pos = elements[3]
    } else {
      stop("Could not parse region argument")
    }
    # region format: "chr1:2000-10000" (1-based, inclusive)
    if (!is.null(end_pos) && start_pos > end_pos) stop("start must be <= end")
    if (!file.exists(fasta)) stop("FASTA file not found: ", fasta)
  
    fai <- paste0(fasta, ".fai")
    is_bgzip <- grepl("\\.(fa|fasta)\\.(gz|bgz)$", fasta, ignore.case = TRUE)
    gzi <- paste0(fasta, ".gzi")
  
    if (!file.exists(fai) || (is_bgzip && !file.exists(gzi))) {
      if (!build_index_if_missing) {
        stop("Missing FASTA index. Need ", fai,
             if (is_bgzip) paste0(" and ", gzi) else "",
             ". Set build_index_if_missing=TRUE to run indexFa().")
      }
      indexFa(fasta)
    }
  
    # Handle chr naming differences (e.g., chr1 vs 1)
    idx <- scanFaIndex(fasta)
    seq_levels <- as.character(seqnames(idx))
    if (!(chr %in% seq_levels)) {
      alt <- if (startsWith(chr, "chr")) sub("^chr", "", chr) else paste0("chr", chr)
      if (alt %in% seq_levels) {
        chr <- alt
      } else {
        stop("Chromosome '", chr, "' not found in FASTA index.")
      }
    }
    
    if (is.null(end_pos)) {
      end_pos = seqlengths(idx)[chr]
    }
    
    param <- GRanges(chr, IRanges::IRanges(start_pos, end_pos))
    dna <- scanFa(fasta, param = param)[[1]]
  
    cg_rel <- start(matchPattern("CG", dna, fixed = TRUE))
    if (length(cg_rel) == 0L) {
      return(data.table(chr = character(), pos = integer(), base = character(), cpg_id = integer()))
    }
  
    c_pos <- start_pos + cg_rel - 1L
    out <- data.table(
      chr = chr,
      pos = as.vector(rbind(c_pos, c_pos + 1L)),     # C then G
      base = rep(c("C", "G"), times = length(c_pos)),
      cpg_id = rep(seq_along(c_pos), each = 2L)
    )
  
    out[, chr := as.factor(chr)]
    out[, base := as.factor(base)]
    return(out)
  }
  
  # For a subset of observations of read orientation on one strand,
  # calculate the suggested cutoff points
  calculate_cutoff <- function(mbias_table) {
    # Ensure we process things in orders
    mbias_table <- mbias_table[order(mbias_table$pos),]
    # get the mean beta of the middle 60% of the data
    n_obs <- nrow(mbias_table)
    if (n_obs < 3) {
      stop("Not enough data points to calculate cutoffs")
    }
    upper_lim <- min(ceiling(n_obs * 0.8), n_obs-1)
    lower_lim <- max(floor(n_obs * 0.2), 2)
    mean_middle <- mean(mbias_table[lower_lim:upper_lim,]$beta)
    upper_ci_max <- max(mbias_table[lower_lim:upper_lim,]$upper_ci)
    lower_ci_min <- min(mbias_table[lower_lim:upper_lim,]$lower_ci)
  
    # Search from the inside out for the first point left/right that is
    # either outside the 95% CI bounds of the middle, or the mean of the
    # middle is outside the CI of the point, or it's more than 0.05 away
    # from the mean of the middle
    check_match <- function(row, mean_middle, lower_ci_min, upper_ci_max) {
      if (is.na(mean_middle) || is.null(mean_middle) || is.na(row$beta) || is.null(row$beta) || is.na(lower_ci_min) || is.null(lower_ci_min) || is.na(upper_ci_max) || is.null(upper_ci_max)) {
        return(FALSE)
      }
  
      if (!row$beta %between% c(lower_ci_min, upper_ci_max)
          && abs(row$beta - mean_middle) > 0.05
          && !mean_middle %between% c(row$lower_ci, row$upper_ci)) {
        return(TRUE)
      }
  
      return(FALSE)
    }
  
    left_cutoff = 0
    right_cutoff = 0
  
    middle = ceiling(n_obs/2)
    for (pos in middle:1) {
      if (check_match(mbias_table[pos,], mean_middle, lower_ci_min, upper_ci_max)) {
        left_cutoff = pos
        break
      }
    }
    for (pos in 1:(n_obs-middle)) {
      if (check_match(mbias_table[middle+pos,], mean_middle, lower_ci_min, upper_ci_max)) {
        right_cutoff = n_obs - (middle+pos) + 1
        break
      }
    }
    return(data.frame(left=as.double(left_cutoff), right=as.double(right_cutoff)))
  }
  
  if (is.na(params$reference) || nchar(params$reference)==0 || !file.exists(params$reference)) {
    stop("V-Bias plot required a reference fasta file")
  }
  
  cpgsInRegion <- lapply(unique(data$chr), function(chrom) {
    res <- cpg_table_from_fasta(region = as.character(chrom), fasta = params$reference)
  })
  cpgsInRegion <- rbindlist(cpgsInRegion)
  cpgsInRegion[, end:=pos]
  setkey(cpgsInRegion, chr, pos, end)
  
  readPositions <- data[, .N, by=.(chr, start, end, orientation)]
  readPositions[, N:=NULL]
  setkey(readPositions, chr, start, end)
  cpgsPerRead <- foverlaps(cpgsInRegion, readPositions, nomatch=NULL)
  cpgsPerRead <- cpgsPerRead[, .N, by=.(chr, start, end, orientation, cpg_id)][, .(ordinal = order(cpg_id), cpg_id), by=.(chr, start, end, orientation)]
  
  rm(readPositions)
  
  # annotate with ordinals for all reads per sample, too
  vdata <- data
  vdata[, ordinal := order(pos_in_read), by=.(read_id, chr, start, end, orientation)]
  
  setkey(vdata, chr, start, end, orientation, ordinal)
  setkey(cpgsPerRead, chr, start, end, orientation, ordinal)
  vdata <- vdata[cpgsPerRead[, .(chr, start, end, orientation, ordinal, cpg_id)], nomatch=NULL]
  rm(cpgsInRegion, cpgsPerRead)
  
  # merge reads into fragments
  vdata <- vdata[, .(start=min(start), end=max(end), pos_in_fragment=min(start+pos_in_read), type=fifelse(any(type=="S"), "S", fifelse(all(type %in% c("U","u")), "U", "M"))), by=.(chr, read_id, cpg_id, insert_size, Orientation=fifelse(bitwAnd(flag, 96)==96 | bitwAnd(flag, 144)==144, "OT", "OB"))]
  vdata[Orientation=="OT", pos_in_fragment := pos_in_fragment - min(start), by=.(chr, read_id)]
  vdata[Orientation=="OB", pos_in_fragment := max(end)-pos_in_fragment, by=.(chr, read_id)]
}
```

```{r plot_vbias, echo=FALSE, fig.height=4, fig.width=7, warning=FALSE}
if (params$plot_vbias) {
  for (chrom in unique(vdata$chr)) {
  
    # main aggregated data (your tiles)
    vdata_plot <- vdata[
      chr == chrom & pos_in_fragment < insert_size,
      .(beta = sum(type == "M") / .N,
        observations = .N),
      by = .(`Pos in fragment` = pos_in_fragment,
             `Insert size` = insert_size,
             Orientation)
    ]
  
    # collapse to counts per Insert size per Orientation (this is what you want plotted)
    side_counts <- vdata_plot[
      , .(obs = sum(observations)),
      by = .(`Insert size`, Orientation)
    ]
    side_counts[, obs_pdf:=obs/sum(obs), by=.(Orientation)]
    
    y_max <- as.numeric(quantile(vdata_plot[["Insert size"]], 0.99,  na.rm = TRUE))
    x_max <- as.numeric(quantile(vdata_plot[["Pos in fragment"]], 0.99, na.rm = TRUE))
  
    p <- ggplot(vdata_plot, aes(x = `Pos in fragment`, y = `Insert size`, fill = beta)) +
      geom_tile() +
  
      # left side panel: counts vs Insert size (one per facet via Orientation)
      geom_ysidebar(
        data = side_counts,
        aes(y = `Insert size`, x = obs_pdf),   # x in side panel = obs, y aligns with main y
        inherit.aes = FALSE,
        stat = "identity",
        orientation = "y"
      ) +
      ggside(y.pos = "left") +
      
      scale_ysidex_reverse() +
      
      scale_fill_gradient2(low = "black", mid = "orange", high = "purple", midpoint = 0.5) +
      scale_y_continuous(limits = c(0, 600), expand = c(0, 0)) +
      scale_x_continuous(limits = c(0, 600), expand = c(0, 0)) +
  
      ggtitle(chrom) +
      facet_wrap(vars(Orientation)) +
  
      theme_minimal() +
      theme(
        panel.grid = element_blank(),
        #axis.title.y = element_blank(),
  
        # make the side panel look like an annotation
        ggside.panel.grid = element_blank(),
        ggside.axis.title.x = element_blank(),
        ggside.axis.text.x  = element_blank(),  # removes the side-panel x tick labels
        ggside.axis.ticks.x = element_blank()   # (optional) removes the tick marks too
      )
  
    print(p)
  }
  rm(vdata, data)
  invisible({gc()})
}
```

## CpG/GC bias

```{r parse_per_pos, include=FALSE}
if (params$plot_gc) {
  suppressPackageStartupMessages({
    library(Rsamtools)
    library(Biostrings)
  })
  
  if ((is.na(params$input_vcf) || is.null(params$input_vcf) || !file.exists(params$input_vcf)) && !is.na(params$input_bam) && !is.null(params$input_bam) && file.exists(params$input_bam)) {
    if (!is.na(params$region) && nchar(params$region) > 0) {
      command <- paste(params$rastair,
    		"call -@", params$threads,
    		"-l", params$region,
    		"--vcf-info-fields SC5,CPG,CPGnovo --vcf-format-fields GT,DP,M5mC --vcf - -c",
    		"-r", params$reference,
    		params$input_bam,
    		"|", params$bcftools, "norm -m-any",
    		"|", params$bcftools, "query -f \"%CHROM\\t%POS\\t%REF\\t%ALT\\t%CPGnovo[\\t%TGT\\t%DP\\t%M5mC]\\n\"")
    } else {
      command <- paste(params$rastair,
    		"call -@", params$threads,
    		"--vcf-info-fields SC5,CPG,CPGnovo --vcf-format-fields GT,DP,M5mC --vcf - -c",
    		"-r", params$reference,
    		params$input_bam,
    		"|", params$bcftools, "norm -m-any",
    		"|", params$bcftools, "query -f \"%CHROM\\t%POS\\t%REF\\t%ALT\\t%CPGnovo[\\t%TGT\\t%DP\\t%M5mC]\\n\"")
    }
  } else if (!is.na(params$input_vcf) && !is.null(params$input_vcf) && file.exists(params$input_vcf)) {
    if (!is.na(params$region) && nchar(params$region) > 0) {
      command <- paste(params$bcftools, 
        "view", params$input_vcf, 
        params$region,
        "|", params$bcftools, "norm -m-any",
  		  "|", params$bcftools, "query -f \"%CHROM\\t%POS\\t%REF\\t%ALT\\t%CPGnovo[\\t%TGT\\t%DP\\t%M5mC]\\n\"")
    } else {
      command <- paste(params$bcftools, 
        "norm -m-any", params$input_vcf,
  		  "|", params$bcftools, "query -f \"%CHROM\\t%POS\\t%REF\\t%ALT\\t%CPGnovo[\\t%TGT\\t%DP\\t%M5mC]\\n\"")
    }
  } else {
    stop("Need either vcf or bam file to create gc bias plot")
  }
  
	all_calls <- fread(cmd=command,
		stringsAsFactors=FALSE,
		header=FALSE)
	setnames(all_calls, c("chr", "pos", "ref", "alt", "isDN", "genotype", "depth", "beta"))
	all_calls[grepl(",", beta, fixed=TRUE), beta := tstrsplit(beta, ",", fixed=TRUE)[[1]]]
	all_calls[, beta := as.numeric(beta)]
	all_calls[, `:=`(chr = as.factor(chr), ref=as.factor(ref), alt=as.factor(alt), genotype=as.factor(genotype))]
	
	if (!is.null(params$is_WGBS) && params$is_WGBS == TRUE) {
    all_calls[, beta := 1-beta]
  }
	site_windows <- all_calls[, .N, by=.(chr = as.character(chr), pos = as.integer(pos))]
  site_windows <- site_windows[order(chr, pos)]
  
  fa <- FaFile(params$reference)
  open(fa)
  on.exit(close(fa), add = TRUE)
  
  chrom_lengths <- seqlengths(fa)
  missing_chroms <- setdiff(unique(site_windows$chr), names(chrom_lengths))
  if (length(missing_chroms) > 0) {
    stop("Chromosome(s) missing from reference FASTA: ", paste(missing_chroms, collapse = ", "))
  }
  
  site_windows[, window_start := pmax(1L, pos - 50L)]
  site_windows[, window_end := pmin(as.integer(chrom_lengths[chr]), pos + 50L)]
  
  gr <- GenomicRanges::GRanges(
    seqnames = site_windows$chr,
    ranges = IRanges::IRanges(start = site_windows$window_start, end = site_windows$window_end)
  )
  
  seqs <- scanFa(fa, param = gr)
  base_counts <- alphabetFrequency(seqs, baseOnly = TRUE)[, c("A", "C", "G", "T"), drop = FALSE]
  
  site_windows[, gc_count := base_counts[, "G"] + base_counts[, "C"]]
  site_windows[, at_count := base_counts[, "A"] + base_counts[, "T"]]
  site_windows[, cpg_count := as.integer(vcountPattern("CG", seqs, fixed = TRUE))]
  site_windows[, window_width := window_end - window_start + 1L]
  rm(seqs)
  site_windows[, chr := factor(chr, levels = levels(all_calls$chr))]
  all_calls[site_windows, on = .(chr, pos), `:=`(
    gc_count = i.gc_count,
    at_count = i.at_count,
    cpg_count = i.cpg_count,
    window_start = i.window_start,
    window_end = i.window_end,
    window_width = i.window_width
  )]
}
```

GC bias:

```{r gc_bias, echo=FALSE, warning=FALSE}
if (params$plot_gc) {
  beta_per_pos <- all_calls[, .(beta=round(mean(beta),1), gc_frac=round(sum(gc_count)/sum(gc_count+at_count), 1), cpgs=round(mean(cpg_count), 1)), by=.(chr, pos)]
  for_plot <- beta_per_pos[, .N, by=.(beta, gc_frac)]
  ggplot(for_plot, aes(x=gc_frac, y=beta, fill=N)) +
  	geom_tile() +
  	scale_fill_gradient2(low = "black", mid = "orange", high = "purple", midpoint = max(for_plot$N)/2) +
  	xlab("%GC in 100bp window") +
    ylab("Methylation") +
    theme_classic()
}
```

CpG bias:

```{r cpg_bias, echo=FALSE, warning=FALSE}
if (params$plot_gc) {
  for_plot <- beta_per_pos[, .N, by=.(beta, cpgs)]
  ggplot(for_plot, aes(x=cpgs, y=beta, fill=N)) +
  	geom_tile() +
  	scale_fill_gradient2(low = "black", mid = "orange", high = "purple", midpoint = max(for_plot$N)/2) +
    xlab("CpGs in 100bp window") +
    ylab("Methylation") +
  	theme_classic()
}
```

Depth bias:

```{r depth_bias, echo=FALSE, warning=FALSE}
if (params$plot_gc) {
  depth_by_cpg <- all_calls[, .(depth=round(mean(depth)),  cpgs=round(mean(cpg_count))), by=.(chr, pos)]
  for_plot <- depth_by_cpg[, .N, by=.(cpgs, depth)]
  ggplot(for_plot, aes(x=cpgs, y=depth, fill=N)) +
  	geom_tile() +
  	scale_fill_gradient2(low = "black", mid = "orange", high="purple", midpoint = max(for_plot$N)/2) +
    xlab("CpGs in 100bp window") +
    ylab("Depth") +
    ylim(0, quantile(depth_by_cpg$depth, 0.99)) +
    theme_classic()
}
```
