Add coverage to a qseaSet using proper pairs and/or high-quality R1 reads

Scan BAM files listed in a qseaSet and add per-region counts and library summaries using proper pairs (with MAPQ and insert-size filters) and, optionally, high-quality unpaired R1 reads extended to a fragment length. Parallelisation via BiocParallel is supported.

addBamCoveragePairedAndUnpaired(
  qs,
  fragmentLength = NULL,
  minMapQual = 30,
  maxInsertSize = 1000,
  minInsertSize = 100,
  minReferenceLength = 30,
  parallel = getMesaParallel(),
  properPairsOnly = FALSE
)

Arguments

qs

qseaSet Input object with valid file_name in its sample table and regions accessible via qsea::getRegions().

fragmentLength

integer(1) or NULL Length to extend unpaired R1 reads. If NULL, estimated as the median proper-pair width (when available). Default: NULL.

minMapQual

integer(1) Minimum MAPQ for either end (proper pairs) or R1 (unpaired). Default: 30.

maxInsertSize

integer(1) Maximum absolute insert size for proper pairs. Default: 1000.

minInsertSize

integer(1) Minimum absolute insert size for proper pairs. Default: 100.

minReferenceLength

integer(1) Minimum reference span for R1s (unpaired). Default: 30.

parallel

logical(1) If TRUE and a multi-core BPPARAM is registered, BAMs are processed in parallel (see BiocParallel). Default: getMesaParallel().

properPairsOnly

logical(1) If TRUE, only proper pairs (plus rescued near-proper pairs) are used. If FALSE, combine proper pairs with high-quality unpaired R1 reads. Default: FALSE.

Value

A qseaSet with updated slots:

• Counts: updated count matrix.

• Library: updated library table with summary metrics.

• Parameters: appended record of MAPQ/size thresholds.

Details

Per-region counts are computed from fragment midpoints to align with qsea’s counting scheme. When properPairsOnly = FALSE, unpaired R1 reads are extended from the 5' end by fragmentLength. If fragmentLength is not supplied, it is estimated from proper pairs in the same BAM.

Parallelisation

Parallel execution uses BiocParallel::bplapply() with the currently registered backend (see BiocParallel::register()). Pass parallel = TRUE and set a multi-core/cluster backend for speed.

See also

getBamCoveragePairedAndUnpairedR1() (internal low-level worker), qsea::getRegions(), qsea::getSampleTable(), BiocParallel::register()

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")

    qs <- data.frame(sample_name = "hESC", file_name = bam,
        group = "Input") %>%
        qsea::createQseaSet("BSgenome.Hsapiens.UCSC.hg19",
            chr.select = "chr22", window_size = 5000) %>%
        addBamCoveragePairedAndUnpaired(fragmentLength = 300)

    # Fragments counted over all windows, and the per-sample QC summary
    colSums(qsea::getCounts(qs))
    getSampleQCSummary(qs)
}
#> ==== Creating qsea set ====
#> restricting analysis on chr22
#> Dividing selected chromosomes of BSgenome.Hsapiens.UCSC.hg19 in 5000nt windows
#> 
#> Attaching package: ‘generics’
#> The following object is masked from ‘package:dplyr’:
#> 
#>     explain
#> The following objects are masked from ‘package:base’:
#> 
#>     as.difftime, as.factor, as.ordered, intersect, is.element, setdiff,
#>     setequal, union
#> 
#> Attaching package: ‘BiocGenerics’
#> The following object is masked from ‘package:dplyr’:
#> 
#>     combine
#> The following objects are masked from ‘package:stats’:
#> 
#>     IQR, mad, sd, var, xtabs
#> The following objects are masked from ‘package:base’:
#> 
#>     Filter, Find, Map, Position, Reduce, anyDuplicated, aperm, append,
#>     as.data.frame, basename, cbind, colnames, dirname, do.call,
#>     duplicated, eval, evalq, get, grep, grepl, is.unsorted, lapply,
#>     mapply, match, mget, order, paste, pmax, pmax.int, pmin, pmin.int,
#>     rank, rbind, rownames, sapply, saveRDS, table, tapply, unique,
#>     unsplit, which.max, which.min
#> 
#> Attaching package: ‘S4Vectors’
#> The following objects are masked from ‘package:dplyr’:
#> 
#>     first, rename
#> The following object is masked from ‘package:utils’:
#> 
#>     findMatches
#> The following objects are masked from ‘package:base’:
#> 
#>     I, expand.grid, unname
#> 
#> Attaching package: ‘IRanges’
#> The following objects are masked from ‘package:dplyr’:
#> 
#>     collapse, desc, slice
#> 
#> Attaching package: ‘Biostrings’
#> The following object is masked from ‘package:base’:
#> 
#>     strsplit
#> Scanning one file at a time. Use e.g. setMesaParallel(nCores = 4) to parallelise this process.
#> Reading hESCs.Input.chr22.bam, with minimum mapping quality 30.
#> No samtools fixmate MQ (mate quality) tags on the bam file, using R1 MAPQ only.
#> Using high quality R1 unpaired reads as well as pairs
#> Setting unpaired fragment lengths to 300.
#> Generated objects
#> Calculated coverage
#> # A tibble: 1 × 3
#>   sample_name valid_fragments  relH
#>   <chr>                 <int> <dbl>
#> 1 hESC                 111515  2.01