Quantify windows by calculating the difference between modification fractions for the positive and negative strand

Calculate sums of per-position Nmod and Nvalid counts for positions that overlap windows in gr, in a strand-specific manner, and calculate FracMod as the ratio between the two. Then calculate the difference between the FracMod for the + and - strand.

Usage

strandDiffFracMod(se, gr, pseudocount = 0)

Arguments

se: A RangedSummarizedExperiment object, as generated by readModBam(..., level = "summary") (i.e. with the assays "Nmod" and "Nvalid").
gr: A GRanges object defining the windows to quantify.
pseudocount: A numeric scalar with a pseudocount that will be added to the "Nmod" and "Nvalid" before calculating the modification fraction.

Value

A RangedSummarizedExperiment object with one row per window defined by gr and the window quantification contained in the assay "FracModDiff".

Author

Charlotte Soneson

Examples

library(GenomicRanges)
library(SummarizedExperiment)
modbamfiles <- system.file("extdata",
                           c("6mA_1_10reads.bam", "6mA_2_10reads.bam"),
                           package = "footprintR")
se <- readModBam(bamfiles = modbamfiles, regions = "chr1:6940000-6940500",
                 level = "summary", modbase = "a", trim = TRUE,
                 BPPARAM = BiocParallel::SerialParam())
s <- seq(6940000, 6940400, by = 100)
windowgr <- GRanges(seqnames = "chr1",
                    ranges = IRanges(start = s, width = 100, names = letters[1:5]))
seStrandDiff <- strandDiffFracMod(se = se, gr = windowgr)
seStrandDiff
#> class: RangedSummarizedExperiment 
#> dim: 5 2 
#> metadata(1): readLevelData
#> assays(6): Nmodpos Nmodneg ... dirNegLog10PValue FracModDiff
#> rownames(5): a b c d e
#> rowData names(0):
#> colnames(2): s1 s2
#> colData names(2): sample modbase
assayNames(seStrandDiff)
#> [1] "Nmodpos"           "Nmodneg"           "Nvalidpos"        
#> [4] "Nvalidneg"         "dirNegLog10PValue" "FracModDiff"      
assay(seStrandDiff, "Nvalidpos")
#>   s1 s2
#> a  0 23
#> b  0 26
#> c  0 27
#> d  0 28
#> e  0 18
assay(seStrandDiff, "Nvalidneg")
#>   s1 s2
#> a 99 32
#> b 91 37
#> c 51 25
#> d 82 40
#> e 76 23