CpG enrichment from a BAM file (MEDIPS-style)

Compute CpG enrichment metrics (relH and GoGe) from aligned reads in a BAM file. Reads are imported directly with Rsamtools, and Biostrings is used to interrogate the reference genome. Optionally exports a fragment-length density plot (PDF) and a serialized RDS with the GRanges of reads.

calculateCGEnrichment(
  file = NULL,
  BSgenome = NULL,
  exportPath = NULL,
  extend = 0,
  shift = 0,
  uniq = 0,
  chr.select = NULL,
  paired = TRUE
)

Arguments

file

character(1) Path to the BAM file. Default: NULL (must be supplied).

BSgenome

character(1) Name of a BSgenome package, e.g. "BSgenome.Hsapiens.NCBI.GRCh38". For GRCh38 and hg19, precomputed distributions are cached within mesa for speed; otherwise, calculateGenomicCGDistribution() is used. Default: NULL.

exportPath

character(1) or NULL Directory in which to write a fragment-length density PDF and an RDS file containing the read GRanges. If NULL, no files are written. Default: NULL.

extend

integer(1) Extension length for single-end reads (used only when paired = FALSE): reads shorter than extend + 1 bases are lengthened to extend + 1 in the 5'->3' direction, as MEDIPS did, and longer reads are left unchanged. Unused for paired reads. Default: 0.

shift

integer(1) Strand-aware offset applied to read positions, used only when paired = FALSE. Unused for paired reads, where the true fragment position is known. Default: 0.

uniq

numeric(1) How to handle duplicate fragments:

• 0 — keep every read.

• 1 — keep at most one read per genomic location and strand.

Any other value, and any logical, is an error. Default: 0.

chr.select

character() or NULL Subset of chromosomes to use. The BAM index is used to restrict the scan when one is present; otherwise the whole file is scanned and the selection applied afterwards. Default: NULL (all chromosomes).

paired

logical(1) Whether the BAM contains paired-end reads. When TRUE, properly paired fragments are reconstructed from the template length; when FALSE, single reads are used (optionally extended via extend). Default: TRUE.

Value

A data.frame with columns:

• file — input BAM path.

• relH — sample relH normalised by genome relH.

• GoGe — sample GoGe normalised by genome GoGe.

• nReads — total reads considered.

• nReadsWithoutPattern — reads lacking "CG" motif.

• n100bpReads — reads with length >= 100 bp.

• n100bpReadsWithoutPattern — reads >= 100 bp lacking "CG".

Details

Reads are scanned for "CG" dinucleotides. Counts are normalised against genome-wide expectations (relH, GoGe).

• For GRCh38 and hg19, cached genomic distributions are bundled with mesa for efficiency.

• For other genomes, calculateGenomicCGDistribution() is invoked.

See also

calculateCGEnrichmentGRanges(), calculateGenomicCGDistribution(), BSgenome

Examples

if (requireNamespace("MEDIPSData", quietly = TRUE) &&
    requireNamespace("BSgenome.Hsapiens.UCSC.hg19", quietly = TRUE)) {
    calculateCGEnrichment(
        file = system.file(
            "extdata",
            "hESCs.Input.chr22.bam",
            package = "MEDIPSData"
        ),
        BSgenome   = "BSgenome.Hsapiens.UCSC.hg19",
        exportPath = tempdir(),
        paired     = FALSE,
        chr.select = "chr22"
    )
}
#> No BAM index found for hESCs.Input.chr22.bam; scanning the whole file and selecting chr22 afterwards.
#> Saving 6.67 x 6.67 in image
#> Ignoring unknown labels:
#> • xlab : "Fragment Length"
#> • ylab : "Density"
#>                                                          file    relH     GoGe
#> 1 /__w/_temp/Library/MEDIPSData/extdata/hESCs.Input.chr22.bam 1.98234 1.050679
#>   nReads nReadsWithoutPattern n100bpReads n100bpReadsWithoutPattern
#> 1 111515                61353           0                         0