compute_fragments_length: Computing the length of the fragments in the .bam file.
In FunChIP: Clustering and Alignment of ChIP-Seq peaks based on their shapes

Description Usage Arguments Details Value Author(s) Examples

Given a .bam file and GRanges object, it computes the positive and negative coverage for each GRanges element, estimates the distance between positive and negative peaks, and finally the fragment length d, i.e. the sum of the length of the reads and the distance between positive and negative peaks. See Details and the package vignette for the description of the method.

1	compute_fragments_length(object, bamf, min.d = 0, max.d = 200)

`object`	GRanges object of length N.
`bamf`	Path to the .bam file used to compute the coverage function. The associated .bam.bai index file must also be present.
`min.d`	integer. Minimum value for the distance between positive and negative peaks. Default is 0.
`max.d`	integer. Maximum value for the distance between positive and negative peaks. Default is 200.

Given a set of n= 1, … N regions, characterized by their positive and negative coverages, the function computes the distance between the positive peak f_{n+} and the negative peak, shifted by δ f_{n-}^{δ}:

D(f_{n+}, f_{n-}^{δ}) = \frac{\| f_{n+} - f_{n-}^{δ}\|^2_{L^2}}{\textrm{width}(\textrm{union}(f_{n+}, f_{n-}^{δ}))}

The function computes the d_{pn} minimizing the distance between postive and negative peaks

d_{pn} = \textrm{argmin}_{δ \in [\code{min.d}, \code{max.d}]} ∑_{n=1 }^N D(f_{n+}, f_{n-}^{δ})

The function returns both the plot of the global distance vs the fragment length d = d_{pn} + r, where r is the length of the reads, and the optimum value for d. r is also estimated from the .bam file as the average of the read lengths.

optimum value of the parameter d, to be used in the pileup_peak method.

Alice Parodi, Marco J. Morelli, Laura M. Sangalli, Piercesare Secchi, Simone Vantini

# load the data
# GRanges object

data(GR100)

# import the .bam file

bamf <- system.file("extdata", "test.bam", package="FunChIP",
                      mustWork=TRUE)
                      
# compute the estimated fragment length

d <- compute_fragments_length(GR[1:10], bamf, min.d = 0, max.d = 200)