---
title: Discrete distributions for the analysis of Next Generation Sequencing (NGS)
  data
author: "Jacques van Helden"
date: 'First version: 2016-12-10; Last update: `r Sys.Date()`'
output:
  slidy_presentation:
    fig_caption: yes
    fig_height: 4
    fig_width: 6
    highlight: tango
    incremental: no
    keep_md: no
    self_contained: yes
    slide_level: 2
    smaller: yes
    theme: cerulean
    toc: yes
    widescreen: yes
  html_document:
    fig_caption: yes
    highlight: zenburn
    self_contained: yes
    theme: cerulean
    toc: yes
    toc_depth: 3
    toc_float: yes
  ioslides_presentation:
    colortheme: dolphin
    fig_caption: yes
    fig_height: 4
    fig_width: 6
    fonttheme: structurebold
    highlight: tango
    incremental: no
    keep_md: no
    slide_level: 2
    smaller: yes
    theme: cerulean
    toc: yes
    widescreen: yes
  pdf_document:
    fig_caption: yes
    highlight: zenburn
    toc: yes
    toc_depth: 3
  beamer_presentation:
    colortheme: dolphin
    fig_caption: yes
    fig_height: 4
    fig_width: 6
    fonttheme: structurebold
    highlight: tango
    incremental: no
    keep_tex: no
    slide_level: 2
    theme: Montpellier
    toc: yes
  word_document:
    fig_caption: yes
    toc: yes
    toc_depth: 3
font-import: http://fonts.googleapis.com/css?family=Risque
font-family: Garamond
transition: linear
---



```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, eval = TRUE, cache = FALSE, message = FALSE, warning = FALSE,
                      comment = "")

## Options to display pretty numbers
library(knitr)
knit_hooks$set(inline = function(x) {
  prettyNum(x, big.mark=" ")
})
options(scipen = 6, digits = 3)

knitr::asis_output("\\footnotesize")

```


# Introduction

## Discrete probabilities and NGS

The advent of Next Generation Sequencing (**NGS**) technologies revived the importance of discrete distributions of probabilities for biologists. 

This tutorial aims at providing a rapid overview of some discrete distributions commonly used to analyse NGS data, and highlight the relationship between them. 

## Overview

| Distribution |  Applications |
|-----------------|-----------------------------------------------|
| Geometric | Local read mapping without mismatch (read extension until first mismatch) |
| Binomial | Global read mapping with a given number of mismatches |
| Negative binomial | Local read mapping with $m$ mismatches (waiting time for $(m+1)^{th}$ mismatch); <br> Detection of differentially expressed genes from RNA-seq data |
| Poisson | ChIP-seq peak calling |
| Hypergeometric | Enrichment of a set of differentially expressed genes for functional classes |


****************************************************************
# Let us experiment first

## The Poisson distribution

The Poisson is a very simple and widely used discrete distribution. 

$$P(X=x) = \frac{e^{-\lambda}\lambda^x}{x!}$$

- represents the probability to observe $x$ successes when expecting $\lambda$ (say "lambda"). 
- expected mean (for a sample of infinite size): $\mu = \lambda$
- expected variance: $\sigma^2 = \lambda$

- **More info**: read the help for the Poisson distribution: `help(Poisson)`

## Exercise -- Poisson distribution


- open  [collective result table](https://docs.google.com/spreadsheets/d/1Kl_0ln0_dZycK17Nqyu44kw9R0dtVp5lflXRtN7pAhA/edit#gid=0)
- login with the email on which you were invited
- each student has been assigned a $\lambda$ comprized between 0.01 and 1000
- draw $rep=1000$ random numbers following a Poisson with this $\lambda$ value
- compute the mean and variance
- fill up the corresponding columns in the collective report

## Solution -- mean and variance of a Poisson random sampling

```{r poisson_sampling_solution}
lambda <- 3
rep <- 1000
x <- rpois(n=rep, lambda=lambda)
mean(x)
var(x)
```

## Replicating an experiment

- read the help for  `runif()` and  `replicate()` 
- make 1000 experiments consisting of the following steps:

    - select at random a $\lambda$ value between $0.5$ and $1000$
    - draw  $n=10$ random numbers following a Poisson with this $\lambda$
    - compute the mean and variance

- plot the relationship between mean and variance for the Poisson distribution

## Solution -- mean to variance relationship for the Poisson distribution

```{r poisson_mean_sampling_distrib}
# ?replicate
## Example of usage of the replicate function
sampling.means <- replicate(n = 10000, mean(rpois(n=10, lambda=3.5)))

hist(sampling.means, breaks=100)
```


```{r mean_var_sampling_distrib}
# This function returns the mean and variance of a 
# random sample drawn from a Poisson distribution
# 
# Parameters:
#   n sample size (number of elements to draw)
#   lambda expectation of the Poisson
rpois.mean.and.var <- function(n, lambda) {
  x <- rpois(n, lambda)
  return(data.frame(mean=mean(x), var=var(x), lambda=lambda, n=n))
}

rpois.mean.and.var(n=10, lambda=3.5)

## Generate a data frame with random Poisson sampling with increasing values of lambda
poisson.stats <- data.frame()
for (i in 1:1000) {
  poisson.stats <- rbind(poisson.stats, rpois.mean.and.var(n=1000, lambda=i))  
}

## Compute the coefficient of variation
poisson.stats$V <- poisson.stats$mean/poisson.stats$var

head(poisson.stats)

print(summary(poisson.stats))


## Summmarize the results with various plots
par(mfrow=c(2,2))
par(mar=c(5,4,1,1))
plot(poisson.stats$mean, poisson.stats$var, col=densCols(poisson.stats$mean, poisson.stats$var), panel.first = grid())
abline(a=0, b=1, col="brown", lwd=2)
plot(poisson.stats$mean, poisson.stats$var, col=densCols(poisson.stats$mean, poisson.stats$var),
    panel.first = grid(), log="xy")
abline(a=0, b=1, lwd=2)
plot(poisson.stats$mean, poisson.stats$V, col=densCols(poisson.stats$mean, poisson.stats$V),
     panel.first = grid())
abline(h=1, col="brown", lwd=2)

# plot(poisson.stats$mean, poisson.stats$V, col=densCols(poisson.stats$mean, poisson.stats$V), 
#      panel.first = grid(), log="xy")
# abline(h=1, col="brown", lwd=2)

hist(poisson.stats$V, breaks=100)

par(mar=c(5,4,4,2))
par(mfrow=c(1,1))
```

## Poisson mean vs variance

```{r plot_poisson_mean_var}
################################################################
# Quentin Ferre's solution with lapply
rep <- 10000
lambda <- runif(n=rep, min = 1, max=1000)
result <- lapply(lambda, rpois, n=10)

# Define a data frame with 2 columns indicating 
# the mean and variance of the random Poisson samples.
rpois.stats <- data.frame(
  mean=unlist(lapply(result, mean)), 
  var=unlist(lapply(result, var))*9/10
)

# Plot the relationship between mean and variance
plot(x=rpois.stats$mean, 
     y=rpois.stats$var, col="grey",
     main="Random poisson sampling",
     xlab="mean", ylab="variance")
grid()
abline(a=0, b=1, col="darkgreen", lwd=2)
```


## Poisson mean vs coefficient of variation

$$V = m / s^2$$

```{r plot_poisson_mean_varcoeff}
# Compute the coefficient of variation
rpois.stats$V <- rpois.stats$mean / rpois.stats$var

# Check the mean of the coefficient of variation
mean(rpois.stats$V)
median(rpois.stats$V)
sum(rpois.stats$mean >= rpois.stats$var)/nrow(rpois.stats)

## PROBLEM HERE: 
## THE VAR IS HIGHER THAN THE MEAN IN MORE CASES THAN EXPECTED
## THE MEAN VAR IS HIGHER THAN THE MEAN MEAN

# Plot the relationship between mean and coefficient of variation
plot(rpois.stats$mean, rpois.stats$V, col="grey",
     main="Random poisson sampling",
     xlab="mean", ylab="variance")
grid()
abline(h=1, col="darkgreen", lwd=2)
```



****************************************************************
# Perfect match probability

## Perfect match probability

We align a library of 50 million short reads of 25 base pairs onto a genome that comprises 23 chromosomes totalling 3 Gigabases. 

For the sake of simplicity, we assume that nucleotides are equiprobable and independently distributed in the genome.

What is the probability to observe the following events by chance?

a. A perfect match for a given read at a given genomic position.
b. A perfect match for a given read anywhere in the genome (searched on two strands).
c. A perfect match for any read of the library at any position of the genome.

d. How many matches do we expect by chance if the whole library is aligned onto the whole genome?

## Perfect match - parameters

Let us define the variables of our problem.
Since we assume equiprobable and independent nucleotides we can define $p$ as probability to observe a match by chance for a given nucleotide.

$$p = P(A) = P(C) = P(G) = P(T) = 0.25$$


```{r perfect_match_parameters}
k <- 25    # Read length
L <- 50e6  # Library size
C <- 23    # Number of chromosomes
G <- 3e9   # Genome size
p <- 1/4   # Matching probability for a nucleotide
```

**Exercise:** use these parameters to compute the matching probability for a read (*solution is on next slide*).

## Perfect match for a given read at a given genomic position

Since we assume independence, the joint probability (probability to match all the nucleotides) is the product of the individual matching probabilities for each nucleotide.


```{r perfect_match_solution}
# Matching probabilty for a given read 
# at a given genomic position
P.read <- p^k  
```

$$P_{\text{read}} = P(n_1 \land n_2 \land \ldots \land n_k) = p^k = `r p`^{`r k`} = `r signif(digits=2, p^k)`$$

This looks a rather small probability. However we need to take into account that this risk will be challenged many times: 

- the size of the genome (`r G`)
- the size of the sequencing library (`r L`)

## Number of genomic alignments

The read will be aligned to each genomic position, but we should keep in mind the following facts. 

1. For each chromosome, we will skip the last `r k-1` positions, since a `r k` bp read cannot be fully aligned there. 
2. We double the number of alignments since we try to map the read on two strands. 

$$N = 2 \sum_{i=1}^C (L_i - k + 1) = 2 \left( G - C (k - 1)  \right)$$

```{r number_of_alignments}
N <- 2 * (G - C * (k - 1))
```

In total, we will thus try to align each read on `r N` genomic positions. 


## Genome-wise matching probability for one read

We reason in 3 steps, by computing the following probabilities. 

| Formula | Rationale |
|---------------------|--------------------------------------------------------|
| $1 -P_{\text{read}} = 1 - p^k$ | no match at a given genomic position |
| $(1 -P_{\text{read}})^N$ | not a single match in the genome |
| $1 - (1 -P_{\text{read}})^N$ | at least one match in the genome |
| | |

```{r pm_one_read_any_position}
P.genomic <- 1 - (1 - P.read)^N
```

This gives $P_{\text{genomic}} = `r P.genomic`$.

## Library-wise probability

We can apply the same reasoning for the library-wise probability. 

| Formula | Rationale |
|-------------------------------|-------------------------------------------|
| $1 -P_{\text{genomic}} = (1 -P_{\text{read}})^N$ | no genomic match for a given read |
| $(1 -P_{\text{read}})^{N L}$ | not a single genomic match in the library |
| $1 - (1 -P_{\text{read}})^{N L}$ | at least one genomic match in the library |
| | |

```{r pm_any_read_any_position}
P.library <- 1 - (1 - P.read)^(N*L)
```

This gives $P_{\text{library}} = `r P.library`$, which should however not be literally interpreted as a certainty, but as a probability so close to $1$ that it cannot be distiguished from it. 

## Expected number of matches


The expected number of matches is the read matching probability mutliplie by the number of matching trials, i.e. $G \cdot L$ since each read will be matched against each genomic position.


$$E(X) = P_{read} \cdot N \cdot L$$

```{r expected_matches}
E <- P.read * N * L
```

In total, we expect `r E` perfect matches by chance for the whole library against the whole genome.

****************************************************************
# Geometric distribution: local alignment without mismatch

## Local alignment until the first mismatch

A local read-mapping algorithm starts by aligning the 5' base of a read, and extends the alignment until either the first mismatch or the end of the read. In the example below, the alignment stops after 11 nucleotides.

```
    ATGCG ACTAG CATAC GAGTG ACTAA
    11111 11111 10
... ATGCG ACTAG CGTTC GACTG ACTAA ...
```

What is the probability to obtain by chance:

a. an alignment of exactly $x=11$ nucleotides (11 matches followed by 1 mismatch)? 
b. an alignment of at least $x=11$ nucleotides (11 matches followed by anything)?

## Local alignment -- parameters


```{r local_no_mismatch}
p <- 0.25  # Matching probability for each nucleotide
x <- 11    # Number of matches before the first mismatch
P.x <- p^x * (1-p)
Pval.x <- p^x
```

$$P(X=`r x`) = p^x (1-p) = 0.25^{`r x`} 0.75 = `r P.x`$$

$$P(X \ge `r x`) = p^x = 0.25^{`r x`} = `r Pval.x`$$


## Geometric distribution


```{r geometric_densiy, echo=FALSE, fig.width=8, fig.height=4, fig.cap="**Geometric distribution.** "}
x <- 0:25  # Series of X values

# Geometric density function
P.x <- dgeom(x = x, prob = p)

# Geometric p-value distribution
Pval.x <- pgeom(q = x-1, prob = p, lower.tail = FALSE)

par(mfrow=c(1,2))
plot(x, P.x, type="h", lwd=2, col="#888888", ylim=c(0,1), las=1, panel.first=grid(), 
     xlab="Nucleotides before first mismatch", ylab="Probabilities")
lines(stepfun(x, c(1, Pval.x)), verticals=TRUE, col="darkgreen")
legend("topright", legend=c("Density: P(X=x)", "P-value: P(X>=x)"), col=c("#888888", "darkgreen"), lwd=2)

plot(x, P.x, type="p", lwd=2, col="#888888", las=1, 
     panel.first=grid(), log="y",
     xlab="Nucleotides before first mismatch", ylab="Probabilities")
lines(stepfun(x, c(1, Pval.x)), verticals=TRUE, col="darkgreen")
legend("topright", legend=c("Density: P(X=x)", "P-value: P(X>=x)"), col=c("#888888", "darkgreen"), lwd=2)

par(mfrow=c(1,1))

```


****************************************************************
# Binomial: global alignment with $m$ mismatches

## Global alignment with mismatches

What is the probability to observe a global alignment with at most $m=3$ mismatches for a given read of 25bp aligned on a particular genomic position?

This question can be formulated as a Bernoulli schema, where each nucleotide is a trial, which can result in either a success (nucleotide match between the read and the genome) or a failure (mismatch). We can label each position of the alignment with a Boolean value indicating whether it maches ($1$) or not ($0$), as examplified below. 

```
    ATGCG ACTAG CATAC GAGTG ACTAA
    11111 11111 10101 11011 11111
... ATGCG ACTAG CGTTC GACTG ACTAA ...
```

At each position, we have a probability of success $p=0.25$, and a probability of failure $q = 1-p = 0.75$.

## Probability to observe exactly $k$ matches

```{r}
n <- 25     # Number of trials, i.e. the length of the alignment
m <- 3      # Maximal number of accepted mismatches
k <- n -m   # Number of matches
p <- 1/4    # Matching probability for one nucleotide
```

Let us denote by $k$ the number of matching residues. 
The probability to observe $k$ successes in a Bernoulli schema with $n$ trials and 

$$P(X=k) = \mathcal{B}(k; n, p) = \binom{n}{k}p^k(1-p)^{n-k} = \frac{n!}{k!(n-k)!}p^k(1-p)^{n-k}$$

## Properties of the binomial distribution

- Mean = $n \cdot p$
- Variance = $n \cdot p \cdot (1 - p)$
- Shape: 

    - i-shaped when $p$ is close to 0, 
    - j-shaped when $p$ is close to 1,
    - bell-shaped for intermediate values of $p$.

## Binomial and perfect match

**Remark**: the perfect match probability seen above is a particular case of the binomial. 

$$P(X=n) = \frac{n!}{n!0!}p^n(1-p)^{n-n} = p^n$$

## Probability of hit with at most $m$ mismatches

We can sum the probabilities for all possible values of matches from $k = n -m$ ($m$ mismatches) to $k = n$ (no mismatch).

$$P(M \le m) = \sum_{k=n-m}^{n} \binom{n}{k}p^k(1-p)^{n-k}$$

## Binomial density

```{r echo=FALSE, fig.width=6, fig.height=4, fig.cap="**Binomial density function**. Alignemnts with at most $m$ mismatches are highlighted in blue. "}
x <- 0:n  # All possible numbers of matches
m <- 3    # Max number of mismatches
hits <- (n-m):n # Hits: at least n-m matches
plot(x, dbinom(x, n, p), type="h", lwd=2, col="grey", 
     xlab="X matching nucleotides", ylab="P(X=x)")
lines(hits, dbinom(hits, n, p), type="h", lwd=2, col="blue")
arrows(n-m, 0.02, n, 0.02, length = 0.1, lwd=2, 
       code = 3, angle=20, col="blue")
```

## Binomial P-value

```{r echo=FALSE, fig.width=6, fig.height=4, fig.cap="**Binomial p-value**. The ordinate indicates the probability to obtain at least $x$ matching nucleotides by chance. "}
x <- 0:n  # All possible numbers of matches
m <- 3    # Max number of mismatches
hits <- (n-m):n # Hits: at least n-m matches
y <- pbinom(x-1, n, p, lower=FALSE)
plot(stepfun(x, c(1, y)), verticals = TRUE,
#plot(x+0.5, y, type="S", lwd=2,
     col="darkgreen", las=1, main="Binomial P-value", 
     xlab="X", ylab="P(X >= x)", log="y",
     panel.first = grid())
pval <- pbinom(n-m-1, n, p, lower=FALSE)
arrows(n-m, pbinom(n-1, n, p, lower=FALSE), n-m, pval, 
       angle = 0, length = 0.1, lwd=2, col="blue")
arrows(n-m, pval, 0, pval, 
       angle = 20, length = 0.1, lwd=2, col="blue")
text(0, 10*pval, 
     labels = paste("P(X >= ", n-m, ") = ", 
                    signif(digits=2, pval), sep=""), pos = 4, col="blue")
```

## Simulated sequences

We can generate random sequences with equiprobable and independent residues from the nucleotide alphabet. 

$$\mathcal{A} = \{A, C, G, T\}$$

```{r random_sequences, echo=FALSE}
# Generate random sequences with equiprobable residues
alphabet <- c("A", "C", "G", "T") # Nucleotide alphabet
n <- 25    # Size of the reads
rep <- 100000 # number of repetitions
# Generate random "sequeences" 
# (actually a matrix of letters)
seq <- data.frame(matrix(nrow=rep, ncol=n, sample(alphabet, size=n * rep, replace = TRUE)))

# Print the 5 first sequences to check result
cat(c(apply(seq[1:5,], 1, paste, collapse=""), "..."), sep="\n")
```

## Match count distribution in simulated sequences

```{r count_matches, echo=FALSE, fig.width=6, fig.height=4, fig.cap="**Global alignment simulation**. A random read is aligned on random sequences. "}
## Count the number of matches between a random read 
## and random sequences of the same length

# Generate a random read
read <- sample(alphabet, size=n, replace = TRUE)

# Define a function counting the number of matching
# residues between two sequences
count.matches <- function(seq1, seq2) {
  return(sum(seq1==seq2))
}

# Compute alignment lengths
matches.per.seq <- apply(seq, 1, count.matches, read)

## Plot an histogram to compare the simulation results 
## and the theoretical distribution
par(mar=c(4,4,0.1,0.1))
hist(matches.per.seq, breaks=seq(from=-0.5, to=n+0.5, by=1), col="grey",
     las=1,
     main="",
     xlab="Number of matches",
     ylab="Sequences")
x <- 0:n
lines(x, rep*dbinom(x, size = n, prob = 0.25), type="h", col="blue", lwd=3)
legend("topright", 
       c(paste(sep="", "random sequences (r=", rep, ")"),
         paste(sep="", "Binom(n=",n,", p=", p,")")), 
       lwd=3,
       col=c("grey", "blue"))
```

## Exercise -- binomial Parameters

Each student will take a custom prior probability ($p$) among the following values: $\{0.001, 0.01, 0.1, 0.25, 0.5, 0.75, 0.9, 0.99, 0.999\}$.

1. Draw 10000 random numbers from a binomial distribution (`rbinom()`) with the custom $p$ and $25$ trials. 
2. Compute the expected mean and variance. 
3. Compute the classical descriptive statistics: mean, variance, standard deviation. 
4. Fill up the form on the [collective result table](https://docs.google.com/spreadsheets/d/1Kl_0ln0_dZycK17Nqyu44kw9R0dtVp5lflXRtN7pAhA/edit#gid=0)
5. Plot an histogram of the numbers drawn.
6. Overlay the theoretical distribution and check the consistency.

## Solution -- binomial 

```{r rbinom}
rand.rep <- 10000 # Random sample size
p <- 0.1          # Prior probability
n <- 25           # Number of trials for the binomial
exp.mean <- n*p   # Expected number of successes
exp.var <- n*p*(1-p)

# Generate random numbers
x <- rbinom(n = rand.rep, size = n, prob = p)

# Compute statistics
stats <- data.frame(p = p, n = n,exp.mean=exp.mean, mean=mean(x), 
  exp.var = exp.var, variance=var(x), sd=sd(x))
kable(stats, digits=4)
```

## Solution -- binomial plot

```{r fig.width=6, fig.height=4, fig.caption="**Random binomial sampling. **"}
hist(x, breaks=(0:26)-0.5, col="grey", main=paste("Binomial simulation, p =", p),
     xlab="Successes", ylab="Frequency", las=1)
lines(0:25, rand.rep*dbinom(x = 0:25, size = n, prob = p), col="blue", lwd=3, type="h")
```


****************************************************************
# Negative binomial: local alignment with at most $m$ mismatches


## Local alignment with mismatches: problem statement


```{r local_alignment_simul, echo=FALSE}
## Simulate local alignment with random sequences. 
## Align a given read with random sequences and measure the 
## length of the longest alignment with at most m mismatches. 

# Max accepted mismatches.
m <- 5
# Note: the alignment stops at m+1

## Define a function counting the length of the alignment 
## until the (m+1)th mismatch
ali.length <- function(seq1, seq2, max.mismatch) {
  mismatches <- which(seq1!=seq2)
  if (length(mismatches) > max.mismatch) {
    return(mismatches[max.mismatch+1]-1)
  } else {
    return (length(seq1))
  }
}

# Compute alignment lengths
ali.lengths <- apply(seq, 1, ali.length, read, m)
max.len <- max(ali.lengths)
```

A local alignment algorithm starts from the 5' end of a read, and stops either at the $x^{th}$ mismatch or when the end of the read is reached. What is the probability to obtain by chance an alignemnt of exactly $25$ nucleotides with exactly $m=`r m`$ mismatches?

This amounts to obtain exactly $k=`r n -m`$ matches and $m=`r m`$ mismatches (in any order), followed by a mismatch at the $(k+m+1)^{th}$ position.

We show here some examples of local alignments with at most `r m` mismatches. Note that the last residue can be either a match (uppercase) or a mismatch (lowercase).

```{r echo=FALSE}
## Show the first alignments
cat(read, sep="")

## Highlight matches in uppercases and mismatches in lowercases
match.cases <- function(seq1, seq2, max.mismatches=NA) {
  seq1 <- unlist(seq1)
  seq2 <- unlist(seq2)
  result <- toupper(seq1)
  mismatches <- seq1 != seq2
  result[mismatches] <- tolower(seq1[mismatches])
  if (!is.na(max.mismatches)) {
    mismatch.pos <- which(seq1 != seq2)
    if (length(mismatch.pos) > max.mismatches) {
      result[mismatch.pos[max.mismatches+1]:length(seq1)] <- ""
    }
  }
  return(paste(result, collapse=""))
}

result <- apply(seq[1:5,], 1, match.cases, read, m)
cat(c(result, "..."), sep="\n")

```


## Number of successes before the $r^{th}$ failure

The **negative binomial** distribution (also called **Pascal distribution**) indicates the probability of the number of successes ($k$) before the $r^{th}$ failure, in a Bernoulli schema with success probability $p$.

$$\mathcal{NB}(k|r, p) = \binom{k+r-1}{k}p^k(1-p)^r$$

This formula is a simple adaptation of the binomial, with the difference that we know that the last trial must be a failure. The binomial coefficient is thus reduced to choose the $k$ successes among the $n-1 = k+r-1$ trials preceding the $r^{th}$ failure. 

## Alternative formulation

It can also be adapted to indicate related probabilities. 

- Number of **failures** ($r$) before the $k^{th}$ **success**.

$$\mathcal{NB}(r|k, p) = \binom{k+r-1}{r}p^k(1-p)^r$$

- Number of **trials** ($n=k+r-1$) before the $r^{th}$ **failure**.

$$\mathcal{NB}(n|r, p) = \binom{n-1}{r-1}p^{n-r}(1-p)^r$$


## Properties of the negative binomial

## Negative binomial density

```{r negbin, echo=FALSE, fig.width=8, fig.height=4, fig.cap="Negative binomial. "}
k <- 0:25    # Number of successes

par(mfrow=c(2,2))
par(mar=c(4,4,1,1))
for (r in c(1, 2, 5, 10)) {
  P.k <- dnbinom(k, size = r, prob = 0.75)
  plot(k, P.k, type="h", panel.first = grid(), lwd=3, ylim=c(0,1), las=1, xlab="Number of successes (k)", ylab="P(K=k)", col="darkblue")
  legend("topright", legend=paste("Last failure: r = ", r))
}
par(mfrow=c(1,1))
```

## Local alignment with simulated sequences

```{r random_simul_local_alignment, fig.width=7, fig.height=4, echo=FALSE}
read.len <- 25  # Max alignment length
x <- m:read.len # All possible alignment lengths
r <- m+1        # Number of failures
n <- x+1        # Number of trials including the closing failure
p <- 0.25       # Matching probability

# Compute and draw the probability with the dnbinom function
p1 <- dnbinom(x-m, size = m+1, prob = 1-p) 

# Compute (approx) the probability manually and check
p2 <- choose(n -1, r-1) *p^(n-r)*(1-p)^r    
# (p1-p2) / (p1+p2) # Check the precision

## plot an histogram of the alignment lengths
hist(ali.lengths, breaks=seq(from=-0.5, to=max.len+1, by=1), col="grey",
     las=1,
     main=paste("Local alignemnt, at most", m, "mismatches"),
     xlab=paste(sep="", "Nucleotides before ", r, "th mismatch"),
     ylab="Number of sequences")
lines(x, rep*p1, type="h", col="blue", lwd=3)
# lines(x,rep*p2, type="h", col="red", lty="dashed")
legend("topright", 
       c(paste(sep="", "random sequences (r=", rep, ")"),
         paste(sep="", "NB(", r, ",", 1-p, ")")), 
       lwd=3,
       col=c("grey", "blue"))
```

## Exercise -- Negative binomial

Each student chooses a value for the maximal number of failures ($r$). 

1. Read carefully the help of the negative binomial functions: `help(NegBinomial)`
2. **Random sampling**: draw of $rep=100000$ random numbers from a negative binomial distribution (`rndbinom()`) to compute the distribution of the number of successes ($k$) before the $r^{th}$ failure. 
3. Compute the expected mean and variance of the negative binomial. 
4. Compute the mean and variance from your sampling distribution. 
5. Draw an histogram with the number of successes before the $r^{th}$ failure. 
6. Fill up the form on the [collective result table](https://docs.google.com/spreadsheets/d/1Kl_0ln0_dZycK17Nqyu44kw9R0dtVp5lflXRtN7pAhA/edit#gid=0)

## Solution to the exercise -- negative binomial

```{r}
r <- 6       # Number of failures
p <- 0.75     # Failure probability
rep <- 100000
k <- rnbinom(n = rep, size = r, prob = p)
max.k <- max(k)
exp.mean <- r*(1-p)/p
rand.mean <- mean(k)
exp.var <- r*(1-p)/p^2
rand.var <- var(k)
hist(k, breaks = -0.5:(max.k+0.5), col="grey", xlab="Number of successes (k)",
     las=1, ylab="", main="Random sampling from negative binomial")
abline(v=rand.mean, col="darkgreen", lwd=2)
abline(v=exp.mean, col="green", lty="dashed")
arrows(rand.mean, rep/20, rand.mean+sqrt(rand.var), rep/20, 
       angle=20, length = 0.1, col="purple", lwd=2)
text(x = rand.mean, y = rep/15, col="purple",
     labels = paste("sd =", signif(digits=2, sqrt(rand.var))), pos=4)
legend("topright", legend=c(
  paste("r =", r), 
  paste("mean =", signif(digits=4, rand.mean)), 
  paste("exp.mean =", signif(digits=4, exp.mean)), 
  paste("var =", signif(digits=4, rand.var)),
  paste("exp.var =", signif(digits=4, exp.var))
  ))
kable(data.frame(r=r, 
                 exp.mean=exp.mean, 
                 mean=rand.mean,
                 exp.var=exp.var,
                 var=rand.var), digits=4)
```

****************************************************************
# Negative binomial for over-dispersed counts


## To be treated in the afternoon !