Merge possibly overlapping peaks

ENCODE has .broadPeak files for ChIP-seq experiments on histone
modifications that, oddly, contain peaks which are subsets of larger
peaks. With IRanges::reduce, overlapping and subsetted peaks are merged.
Support for .broadPeak and .narrowPeak files

Files downloaded directly from ENCODE tend to be .broadPeak or
.narrowPeak format. Both are of the BED6 + X variety. Only the first
three columns (chrom, start, end) are used for analysis.
New enrichment method: ChIP-Enrich Ratio

ChIP-Enrich Ratio is an additional enrichment method whose model is:
goterm ~ ratio + spline, where ratio refers to the percentage of a
gene locus which is covered by a ChIP-seq peak. It is meant for use
with histone modifications or transcription factors with broader
ChIP-seq peaks.

Three new functions are implemented for this method:
1. assign_peak_segments(...) determines the overlap between
the peaks and gene loci in base pairs. Note, a peak overlapping more
than one locus has it's base pairs proportionally assigned.
2. test_gam_ratio(...) builds the model and performs the
enrichment.
3. plot_gene_versus_coverage(…) plots the log10 length of the gene locus
against the percentage covered by a peak. It is a part of the QC plots.

Two functions are modified:
1. num_peaks_per_gene(assigned_peaks,locusdef,…,ratio) will function as
with method=='chipenrich' when ratio==F. When ratio==T an additional
column is added to the peaks per gene table indicating the ratio of
the gene locus covered by a peak.
2. chipenrich(…) is modified to use assign_peak_segments(…) and
num_peaks_per_gene(...,ratio=T) when method=='ce.ratio'.
Implemented chipenrich approximating enrichment method

*Added comments to some functions to give downstream/upstream context.
*Added test_approx to SUPPORTED_METHODS and METHOD_NAMES.
*num_peaks_per_gene(…) no longer includes calculations for ce.ratio
method. The functionality is broken out to calc_peak_gene_overlap(…).
Tested broken out versus not and obtained exact same results.
*Added calc_approx_weights(…) which determines the appropriate weights,
the probability of a peak, and the resid.dev based on the smoothing
spline. For use in the approximating method.
*Added test_approx(…). Still unclear on the meaning of nwp in the
prototype code. Also unsure of the correction case implementation.
*Also of note is that the chipenrich model is peak ~ goterm +
s(log10_length,bs='cr’), but for the approximating method peak is
replaced by peak*weight and is no longer binary as necessary for
logistic regression so the approximating model is goterm ~ peak*weight.
The variables are swapped.
Bug fixes

*For some reason, adding library(rms) at the top of main.R does not
load the package, but putting it in test_approx(…) works like a charm.
*Removed leftover message(i) statement in the set.gene loop.
*Fixed absence of num_peaks column in peak correction.
*Fixed longstanding bug where reading in a .narrowPeak gives a message
saying “Reading broadPeak…”. Silly.
Change p-value calculation in approximate method

1-pchisq(z^2,1) has a lower limit that using lower.tail=F doesn’t have.
Report unique peak ID in *_peaks.tab

Create a hash of the chrom, peak_start, peak_end and assign an ID to
the unique combination. The hash gets dropped in the chipenrich
function.

Also added a cap of 1 on the ratio of peak overlap.
Enable randomizations

There are two types of randomizations: rndall and rndlength.

rndall takes the peaks per gene table (after possible column additions
from different methods) and permutes the rows over all the genes.

rndlength does the same, but orders the gene loci by length, bins into
groups of 100, and permutes the peaks per gene table within those 100
bins to preserve the locus lengths of the genesets.

Randomizations are accessed by appending either “_rndall” or
“_rndlength” to any gene set in the genesets parameter of chipenrich(…).

Supporting functions randomize_ppg_all(…) and randomize_ppg_length(…)
are added, and the appropriate parsing of “_rndall” and “_rndlength”
are added within chipenrich(…).
Fixed bug in assigning unique Peak IDs
ce.ratio is now broadenrich

All references are changed, method name to use is ‘broadenrich’.
Binned locus length vs % coverage QC plot

More clear than the non-binned version previously implemented.
chipenrich, broadenrich, and chipapprox are multicore

We implemented multicore versions of the chipenrich, broadenrich, and
chipapprox methods using the multicore R package. Note, Windows does
not support the multicore R package, so these methods will run as
single core for Windows users.

Bug fix: A ‘continuity gene’ was introduced for genesets if all genes
in the geneset had a peak to help the logistic regression converge. The
‘continuity gene’ mistakenly persisted for all genesets thereafter,
rather than only for the required geneset, slightly altering enrichment
p-values. This is no longer the case.
multicore and plotting changes

To avoid having to add multicore to the NAMESPACE, which is not
supported in Windows, the library is loaded within the functions which
use it.

plot_gene_versus_coverage() and ..plot_gene_versus_coverage() are now
just plot_gene_coverage() and ..plot_gene_coverage().

When method = 'broadenrich' only output the ..plot_gene_coverage()
plot. The other plots are not relevant.
Documentation changes and bumped package version

*Changes to and creation of R help files. Updated the vignette for the
package.
*Altered the NAMESPACE to recognize plot_gene_coverage.
*Updated DESCRIPTION to include multicore as enhancer package. Bumped
up version number in keeping with Bioconductor versioning guidelines
(as far as I understand).
*News summarizes changes to the package.
*Updated install_chipenrich.R file. Tested on Mac, Windows, and Linux.
Merge with master

*Pulled changes made by Ryan since dev branch was created.
*Conflicts in main.R settled with preference to dev branch.
Specifically, did not include GO terms associated with genes in
peaks-per-gene table.
Broad-Enrich QC Plot axis label

% locus covered is inaccurate. Changed to proportion of locus covered...
Update to v1.1.2

All instances of merge_all are replaced with non-recursive Reduce. This
solves some recursion error problems for large gene sets (such as all
of GO, or MeSH).