---
title: "Simple model of ChIPseq"
author: Michael Stadler, michael.stadler@fmi.ch
date: "April 24, 2024"
output: html_document
toc: false
runtime: shiny
---

```{r encoding, echo=FALSE}
options(encoding="UTF-8")
```

A simple theoretical model can help to think about the process (described in Mikkelsen et al., Nature 2007 [PubMed](https://www.ncbi.nlm.nih.gov/pubmed/17603471)):

![ChIPseq_model](figures/11/chipseq.png)


* genome is divided into $N$ bins of size $b$
* fraction $f$ of these bins contain feature of interest  
* uniform sampling of reads: number of reads per bin $\sim$ Poisson$(\lambda)$
* immuno-precipitation enriches by factor $e$  
  (on average, bins with/without features have $eM$ and $M$ reads):  
  $\lambda_{background} = M$ and $\lambda_{enriched} = eM$  
* when sequencing a total of $R$ reads:  
  $M = \frac{R}{N(ef+(1-f))}$

Use the interactive Shiny app below (for more information see [shiny.rstudio.com](http://shiny.rstudio.com/)) to play with the model parameters and estimate the expected distribution of read counts in genome-tiling windows:

```{r shiny, echo=FALSE}
genomeSizes <- c("Yeast (12Mb)"=12e6, "Worm (100Mb)"=100e6, "Human (3.3Gb)"=3.3e9)
binSizes <- c(100,1000,10000,100000)
inputPanel(
  selectInput("genome", "Genome:",  choices=names(genomeSizes)),
  selectInput("b", "Binsize (b):",  choices=binSizes),
  sliderInput("f", "Fraction of bins with feature (f):", min=0.002, max=0.5, value=0.002, ticks=FALSE),
  selectInput("e", "IP fold-enrichment (e):",  choices=c(2,5,10,20)),
  selectInput("R", "Number of reads (R):",  choices=c(1e6,2e6,5e6,10e6,20e6,50e6,100e6)),
  #sliderInput("t", "Threshold:", min=0, max=200, value=12)
  numericInput("t", "Threshold:", min=0, max=1000, value=12)
  )


renderPlot({
  e <- as.numeric(input$e)
  f <- as.numeric(input$f)
  th <- as.numeric(input$t)
  N <- round(genomeSizes[input$genome] /as.numeric(input$b))
  M <- as.numeric(input$R)/(N*(e*f+(1-f)))
  TP <- f*N*ppois(th, lambda=e*M, lower.tail=FALSE)
  FN <- f*N*ppois(th, lambda=e*M, lower.tail=TRUE)
  FP <- (1-f)*N*ppois(th, lambda=M, lower.tail=FALSE)
  TN <- (1-f)*N*ppois(th, lambda=M, lower.tail=TRUE)
  FDR <- FP/(TP+FP)
  sens <- TP/(TP+FN)
  xs <- 0:max(10,round(2.2*e*M))
  mp <- barplot(rbind(dpois(xs, lambda=M),
                      dpois(xs, lambda=e*M),
                      f/(1-f)*dpois(xs, lambda=e*M)),
                names.arg=xs, legend.text=c("background","enriched","enriched (proportional)"),
                beside=TRUE, col=c("red","green3","lightgreen"), xlab="Reads per bin", ylab="Density of bins")
  thx <- colMeans(mp)[ifelse(th < xs[1] || th > xs[length(xs)], NA, findInterval(th, xs))] + 0.5*diff(colMeans(mp)[1:2])
  abline(v=thx, lty=2, col="black")
  text(x=thx+par('cxy')[1], y=par('usr')[4], adj=c(0,1.5), labels=sprintf("FDR = %.1f%%\nsensitivity = %.1f%%",100*FDR,100*sens))
})
#}, width="auto", height=400)
```


## References
