Add MEDIPS-style enrichment metrics to a qseaSet sample table

For each sample, compute MEDIPS-style CpG enrichment metrics (relH, GoGe, and read counts) from its BAM and append the results to the qseaSet sample metadata. Supports both pulldown (default) and input libraries.

addMedipsEnrichmentFactors(
  qseaSet,
  exportPath = NULL,
  nonEnrich = FALSE,
  extend = 0,
  shift = 0,
  uniq = 0,
  chr.select = NULL,
  paired = TRUE,
  file_name = "file_name",
  nCores = 1
)

Arguments

qseaSet

qseaSet Input object whose sample table contains BAM paths.

exportPath

character(1) or NULL Directory to export per-sample fragment-length PDFs and read-GRanges RDS files. If NULL, no files are written. Default: NULL.

nonEnrich

logical(1) If TRUE, treat samples as Input libraries (columns appended under @libraries$input_file). If FALSE, treat as Pulldown (columns appended under @libraries$file_name). Default: FALSE.

extend

integer(1) Passed to calculateCGEnrichment(). Extension length for unpaired reads: shorter reads are lengthened to extend + 1 bases, longer ones left unchanged. Default: 0.

shift

integer(1) Passed to calculateCGEnrichment(). Shift applied to unpaired read positions; unused for paired reads. Default: 0.

uniq

numeric(1) Passed to calculateCGEnrichment(). Duplicate handling: 0 keeps all reads, 1 keeps one per genomic location. Default: 0.

chr.select

character() or NULL Passed to calculateCGEnrichment(); subset of chromosomes to analyse. Default: NULL (all chromosomes).

paired

logical(1) Whether BAMs are paired-end. Paired reads are read as the fragment span between the first mate's start and its mate's end. Default: TRUE.

file_name

character(1) Column name in the sample table holding BAM paths when nonEnrich = FALSE. When nonEnrich = TRUE, the column input_file is used instead. Default: "file_name".

nCores

integer(1) Number of parallel cores for parallel::mclapply(). Set to 1 for serial. Default: 1.

Value

A qseaSet with new MEDIPS-style metrics appended:

• Sample/library metrics (added under the appropriate library frame): relH, GoGe, nReads, nReadsWithoutPattern, n100bpReads, n100bpReadsWithoutPattern.

• Provenance: any export artefacts (PDF/RDS) written to exportPath if provided.

Details

Internally calls calculateCGEnrichment() per sample to derive relH/GoGe and counts, then merges the results into the @libraries slot ($file_name for pulldown, $input_file for input when nonEnrich = TRUE). Parallelisation uses parallel::mclapply(); on non-Unix systems, nCores > 1 is ignored.

See also

calculateCGEnrichment(), calculateCGEnrichmentGRanges()

Examples

# A single-end MeDIP-seq input BAM on chr22 (hg19), from MEDIPSData
if (requireNamespace("MEDIPSData", quietly = TRUE) &&
    requireNamespace("BSgenome.Hsapiens.UCSC.hg19", quietly = TRUE)) {
    bam <- system.file(
        "extdata",
        "hESCs.Input.chr22.bam",
        package = "MEDIPSData"
    )

    # The enrichment factors are stored next to the coverage, so add
    # coverage first
    qs <- data.frame(
        sample_name = "hESC_Input_chr22",
        file_name = bam,
        group = "Input"
    ) %>%
        qsea::createQseaSet("BSgenome.Hsapiens.UCSC.hg19",
            chr.select = "chr22", window_size = 5000) %>%
        qsea::addCoverage(fragment_length = 300, uniquePos = FALSE,
            paired = FALSE) %>%
        addMedipsEnrichmentFactors(chr.select = "chr22", paired = FALSE)

    # Copy the factors into the sample table to inspect them
    qs %>%
        addLibraryInformation() %>%
        getSampleTable() %>%
        dplyr::select(sample_name, relH, GoGe, nReads, nReadsWithoutPattern)
}
#> ==== Creating qsea set ====
#> restricting analysis on chr22
#> Dividing selected chromosomes of BSgenome.Hsapiens.UCSC.hg19 in 5000nt windows
#> Reading bam alignment /__w/_temp/Library/MEDIPSData/extdata/hESCs.Input.chr22.bam
#> warning: reduced efficiency, no bam index found
#> Number of imported sequencing fragments: 111515
#> Calculating short read coverage for genome wide windows.
#> Number of selected sequencing fragments: 111515
#> Adding Medips Enrichment factors to 1 enriched samples, using 1 cores.
#> No BAM index found for hESCs.Input.chr22.bam; scanning the whole file and selecting chr22 afterwards.
#>                       sample_name    relH     GoGe nReads nReadsWithoutPattern
#> hESC_Input_chr22 hESC_Input_chr22 1.98234 1.050679 111515                61353