calcPhasogram: Calculate phasograms (same strand alignment distances).

In fmicompbio/swissknife: Handy code shared in the FMI CompBio group

Calculate phasograms (same strand alignment distances).

Description

Calculate the frequencies of same strand alignment distances, for example from MNase-seq data to estimate nucleosome repeat length. Distance calculations are implemented in C++ (calcAndCountDist) for efficiency.

Usage

calcPhasogram(fname, regions = NULL, rmdup = TRUE, dmax = 3000L)

Arguments

fname	character vector with one or several bam files. If multiple files are given, distance counts from all will be summed.
regions	GRanges object. Only alignments falling into these regions will be used. If NULL (the default), all alignments are used.
rmdup	logical(1) indicating if duplicates should be removed. If TRUE (the default), only one of several alignments starting at the same coordinate is used.
dmax	numeric(1) specifying the maximal distance between same strand alignments to count.

Value

integer vector with dmax elements, with the element at position d giving the observed number of alignment pairs at that distance.

Author(s)

Michael Stadler

References

Phasograms were originally described in Valouev et al., Nature 2011 (doi:10.1038/nature10002). The implementation here differs in two ways from the original algorithms:

1. It does not implement removing of positions that have been seen less than n times (referred to as a n-pile subset in the paper).

2. It does allow to retain only alignments that fall into selected genomic intervals (regions argument).

See Also

estimateNRL to estimate the nucleosome repeat length from a phasogram, plotPhasogram to visualize an annotated phasogram, calcAndCountDist for low-level distance counting.

Examples

if (requireNamespace("GenomicAlignments", quietly = TRUE) &&
    requireNamespace("Rsamtools", quietly = TRUE)) {
    bamf <- system.file("extdata", "phasograms", "mnase_mm10.bam",
                        package = "swissknife")
    pg <- calcPhasogram(bamf)
    print(estimateNRL(pg, usePeaks = 1:4)[1:2])
    plotPhasogram(pg, usePeaks = 1:4, xlim = c(0,1000))
}

fmicompbio/swissknife documentation built on June 11, 2025, 4:17 p.m.