---
title: "QC Report"
date: "`r Sys.Date()`"
output: 
  html_document:
    toc: true
    theme: united
params:
  input_bgz: NULL
  output_dir: "."
  region: NULL # By default there's no region
  tabix: "tabix" # Assume tabix is available in PATH
  include_flags: 3
  exclude_flags: 3852
  read_len: NULL
---

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

# 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}
# Load data
if (!is.na(params$region) && nchar(params$region) > 0) {
  data = fread(cmd=paste(params$tabix, "-h", params$input_bgz, params$region), header=TRUE, stringsAsFactors=FALSE)
} else {
  data = fread(params$input_bgz, header=TRUE, stringsAsFactors=FALSE)
}
setnames(data, 1, "chr")
data[, chr := as.factor(chr)]

# 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
}
  
split_pos <- function(mod, unmod, snp) {
  mod_pos <- as.numeric(tstrsplit(mod, ","))
  unmod_pos <- as.integer(tstrsplit(unmod, ","))
  snp_pos <- as.integer(tstrsplit(snp, ","))
  return(list(pos_in_read = c(mod_pos, unmod_pos, snp_pos), type=as.factor(c(rep("mod", length(mod_pos)), rep("unmod", length(unmod_pos)), rep("snp", length(snp_pos))))))
}
data <- data[, split_pos(mod_cpgs, unmod_cpgs, snp_cpgs), by=.(chr, start, end, read_id, mapq, orientation, insert_size, read_length, flag)]

mbias_table <- data[, .(mod=sum(type=="mod"), unmod=sum(type=="unmod")), 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}
merge_pos <- function(dt) {
  if (nrow(dt) < 2) {
    return(NULL)
  }
  read_len = dt$end-dt$start
  first_in_pair = dt[, bitwAnd(flag, 64)==64]
  if (!any(first_in_pair)) {
    return(NULL)
  }
  Orientation = head(dt[, fifelse(bitwAnd(flag, 80)==80 | bitwAnd(flag, 160)==160, "OT", "OB")], 1)
  pos_in_fragment = fifelse(first_in_pair, dt$pos_in_read, read_len-dt$pos_in_read)
  new_type <- fifelse(first_in_pair, dt$type, rev(dt$type))
  pos_in_fragment[!first_in_pair] = pos_in_fragment[!first_in_pair] + dt$insert_size[!first_in_pair]-read_len[!first_in_pair]
  duplicates <- pos_in_fragment[!first_in_pair] < max(read_len[first_in_pair])
  pos_in_fragment <- pos_in_fragment[!duplicates]
  new_type <- new_type[!duplicates]
  return(dt[, .(type=new_type, pos_in_fragment=pos_in_fragment, insert_size=rep(min(insert_size), length(new_type)), fragment_length=rep(max(dt$end)-min(dt$start), length(new_type)), mapq=rep(mean(mapq), length(new_type)), Orientation=rep(Orientation, length(new_type)))])
}
vdata <- data[, merge_pos(.SD), by=.(chr, read_id)]
```

```{r plot_vbias, echo=FALSE, fig.height=4, fig.width=7}
for (chrom in unique(data$chr)) {
  p <- ggplot(vdata[chr==chrom, .(beta=sum(type=="mod")/.N), by=.(`Pos in fragment`=pos_in_fragment, `Insert size`=insert_size, Orientation)], aes(x=`Pos in fragment`, y=`Insert size`, color=beta)) + scale_color_gradient(low="darkblue", high="red") + geom_point(size=0.001, alpha=0.2) + ggtitle(chrom) + facet_wrap(vars(Orientation))
  print(p)
}
```

