fdrEnrichedCounts {htSeqTools} | R Documentation

## Posterior probability that a certain number of repeats are higher than expected by chance.

### Description

Given a vector of number of repeats (e.g. there are 100 sequences appearing once, 50 sequences appearing twice etc.) the function computes the false discovery rate that each number of repeats is unusually high.

### Usage

    fdrEnrichedCounts(counts,use=1:10,components=0,mc.cores=1)

### Arguments

`counts` | vector with observed frequencies. The vector must have names. `tabDuplReads` function can be used for this purpose.
`use` | number of repeats to be used when estimating the null distribution. The number of repeats expected if no unusually high repeats are present. The first 10 are used by default.
`components` | number of negative binomials that will be used to fit the null distribution. The default value is 1. This value has to be between 0 and 4. If 0 is given the optimal number of negative biomials is chosen using the Bayesian information criterion (BIC)
`mc.cores` | number of cores to be used to compute calculations. This parameter will be passed bt to `mclappply`

### Details

The null distribution is a combination of n negative binomials where. n is assigned through the `components` parameter. If `components` is equal to 0 the optimal number of negative binomials is choosen using the Bayesian information criterion (BIC). The parameters of the null distribution are estimated from the number of observations with as many repeats as told in the `use` parameter. If use is 1:10 the null distribution will be estimated using repeats that appear 1 time, 2 times, ... or 10 times.

False discovery rate for usually high number of repeats is done following an empirical Bayes scheme similar to that in Efron et al. Let f0(x) be the null distribution, f(x) be the overall distribution and (1-pi0) the proportion of unusually high repeats. We assume the two component mixture f(x)= pi0 f0(x) + (1-pi0)f1(x). Essentially, f(x) is estimated from the data (imposing that f(x) must be monotone decreasing after its mode using `isoreg` from packabe `base`, to improve the estimate in the tails). Currently pi0 is set to 1, i.e. its maximum possible value, which provides an upper bound for the FDR. The estimated false discovery rate for enrichment is `1-pi0*(1-cumsum(f0(x)))/(1-cumsum(f(x)))`. A monotone regression (isoreg) is applied to remove small random fluctuations in the estimated FDR and to guarantee that it decreases with x.

### Value

`data.frame` with the following columns:

`pdfH0` | vector with pdf under the null hypothesis of no enrichment
`pdfOverall` | vector with pdf for mixture distribution
`fdrEnriched` | vector with false discovery rate that each count is significantly enriched

### References

Ji et al. An integrated software system for analyzing ChIP-chip and ChIP-seq data. Nature Biotechnology, 2008, 26, 1293-1300.

Efron et al. Empirical Bayes Analysis of a Microarray Experiment, Journal of the American Statistical Association, 2001, 96, 1151-1160.

### Examples

    #Generate 1000 sequences repeated once, on the average
    nrepeats <- c(rpois(10^4,1),rpois(10,10))
    nrepeats <- nrepeats[nrepeats>0]
    counts <- table(nrepeats)
    barplot(counts) -> xaxis #observe bimodality around 10
    fdrest <- fdrEnrichedCounts(counts,use=1:5,components=1)
    cutoff <- xaxis[which(fdrest$fdrEnriched<0.95)[1]]
    abline(v=cutoff,col=2)
    text(cutoff,counts[1]/2,'cut-off',col=2)
    head(fdrest)

[Package htSeqTools version 1.28.3 Index]