Skip to contents

Given a SummarizedExperiment with modified and total base counts in assays "Nmod" and "Nvalid", perform a pairwise statistical test for differential modification.

Usage

getDifferentiallyModifiedWindows(
  se,
  assayNameMod = "Nmod",
  assayNameValid = "Nvalid",
  groupCol = "group",
  verbose = FALSE
)

Arguments

se

SummarizedExperiment, for example returned by quantifyWindowsInRegion. It is expected to at least contain assays for modified and total counts (given by assayNameMod and assayNameValid) and a colData column that defines the groups (given by groupCol).

assayNameMod, assayNameValid

Character scalars that give the assay names in se containing the modified and total counts, respectively.

groupCol

Character scalar giving the column in colData(se) that defines the groups of samples to be compared.

verbose

Logical scalar. If TRUE, report on progress.

Value

The GRanges object constructed from the topTags output obtained for the statistical analysis, with an additional column named "dirNegLog10PValue", calculated as the sign of the logFC multiplied with the -log10(PValue).

Author

Panagiotis Papasaikas, Sebastien Smallwood, Charlotte Soneson, Michael Stadler

Examples

modbamfiles <- system.file("extdata",
                           c("6mA_1_10reads.bam", "6mA_1_10reads.bam",
                             "6mA_2_10reads.bam", "6mA_2_10reads.bam"),
                           package = "footprintR")
se <- quantifyWindowsInRegion(bamfiles = modbamfiles,
                              region = "chr1:6940000-6955000", modbase = "a",
                              BPPARAM = BiocParallel::SerialParam())
se$group <- c("group1", "group1", "group2", "group2")
gr <- getDifferentiallyModifiedWindows(se, groupCol = "group")
class(gr)
#> [1] "GRanges"
#> attr(,"package")
#> [1] "GenomicRanges"
head(gr)
#> GRanges object with 6 ranges and 6 metadata columns:
#>       seqnames          ranges strand |     logFC    logCPM        LR
#>          <Rle>       <IRanges>  <Rle> | <numeric> <numeric> <numeric>
#>   [1]     chr1 6940000-6940023      * |  6.696961   12.9861 22.480673
#>   [2]     chr1 6940012-6940035      * |  7.018526   13.4242 29.825436
#>   [3]     chr1 6940024-6940047      * |  5.675384   13.7867 13.472295
#>   [4]     chr1 6940036-6940059      * |  0.585622   13.7867  0.166427
#>   [5]     chr1 6940048-6940071      * |  1.214138   13.6466  1.284355
#>   [6]     chr1 6940060-6940083      * | -1.165075   13.5237  1.256727
#>            PValue         FDR dirNegLog10PValue
#>         <numeric>   <numeric>         <numeric>
#>   [1] 2.12269e-06 1.44343e-04          5.673114
#>   [2] 4.72749e-08 6.42938e-06          7.325370
#>   [3] 2.42112e-04 8.23181e-03          3.615984
#>   [4] 6.83307e-01 1.00000e+00          0.165384
#>   [5] 2.57091e-01 1.00000e+00          0.589913
#>   [6] 2.62272e-01 1.00000e+00         -0.581249
#>   -------
#>   seqinfo: 1 sequence from an unspecified genome; no seqlengths