diff --git a/DESCRIPTION b/DESCRIPTION index d1bb933..3334c2a 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: SingleMoleculeGenomicsIO Title: Tools for reading single molecule genomics data -Version: 0.1.8 +Version: 0.1.9 Authors@R: c( person("Michael", "Stadler", , "michael.stadler@fmi.ch", role = c("aut", "cre"), comment = c(ORCID = "0000-0002-2269-4934")), @@ -67,4 +67,4 @@ Config/testthat/edition: 3 Encoding: UTF-8 Roxygen: list(markdown = TRUE) SystemRequirements: GNU make -Config/roxygen2/version: 8.1.0 +Config/roxygen2/version: 8.1.1 diff --git a/NAMESPACE b/NAMESPACE index 68bfb1b..32b4787 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -36,6 +36,7 @@ export(subsetReads) importFrom(BSgenome,getSeq) importFrom(BiocGenerics, "%in%", + "strand<-", as.data.frame, as.list, cbind, diff --git a/NEWS.md b/NEWS.md index 1b3aaae..4d238c2 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,8 @@ +# SingleMoleculeGenomicsIO 0.1.9 + +* Allow `ignore.strand=FALSE` in `expandSEToBaseSpace` +* Add `flank` argument to `expandSEToBaseSpace` + # SingleMoleculeGenomicsIO 0.1.8 * Add `getBaseCoverageForBam` and `countAlignmentsInBam` functions diff --git a/R/expandSEToBaseSpace.R b/R/expandSEToBaseSpace.R index ac4a1a5..f1d7bf5 100644 --- a/R/expandSEToBaseSpace.R +++ b/R/expandSEToBaseSpace.R @@ -7,6 +7,10 @@ #' regions can be specified as character scalars (e.g. "chr1:1200-1300") #' that can be coerced into a \code{GRanges} object. If \code{NULL} (the #' default), \code{regions} is set to the range of the data in \code{se}. +#' @param flank Integer scalar. If \code{regions=NULL} and \code{flank} is +#' non-zero, \code{flank} positions will be added on both sides of each +#' range in the derived \code{regions}. This is useful to add a buffer of +#' NA values between observed values in non-overlapping regions. #' @param seqinfo \code{NULL} or a \code{\link[Seqinfo]{Seqinfo}} object #' containing information about the set of genomic sequences (chromosomes). #' Alternatively, a named numeric vector with genomic sequence names and @@ -16,7 +20,15 @@ #' expand to base space. Only these assays will be present in the returned #' object. #' @param ignore.strand A logical scalar defining whether to ignore the strand -#' of the \code{rowRanges(se)}. +#' of the \code{rowRanges(se)}. If \code{FALSE}, \code{se} is first +#' filtered (using \code{filterPositions(..., filters = "repeated.positions")}) +#' to resolve any positions represented on both strands. +#' @param assayNameAmbig A character scalar indicating the assay to use to +#' decide which row to retain if multiple rows represent the same genomic +#' position (on different strands). The row with the largest row sum in +#' this assay is retained. Only used if \code{ignore.strand} is \code{FALSE}. +#' @param verbose Logical scalar indicating whether to print out progress +#' messages. #' #' @importFrom SummarizedExperiment SummarizedExperiment colData assays #' @importFrom S4Vectors metadata @@ -24,12 +36,17 @@ #' @importFrom GenomicRanges GPos #' @importFrom utils modifyList #' @importFrom cli cli_abort -#' @importFrom BiocGenerics start end nrow ncol colnames unstrand +#' @importFrom BiocGenerics start end nrow ncol colnames unstrand strand strand<- #' @importFrom SparseArray NaArray nnawhich #' @importFrom Seqinfo seqnames seqinfo #' -#' @return A single-base resolution \code{\link[SummarizedExperiment]{RangedSummarizedExperiment}} -#' corresponding to \code{se}. +#' @return A single-base resolution +#' \code{\link[SummarizedExperiment]{RangedSummarizedExperiment}} +#' corresponding to \code{se}. If \code{ignore.strand} is \code{TRUE}, +#' the return object is unstranded. If \code{ignore.strand} is \code{FALSE}, +#' all positions present in the input object will retain the assigned +#' strand, while all other positions will be unstranded. +#' #' @author Charlotte Soneson, Michael Stadler #' #' @examples @@ -45,11 +62,16 @@ #' @export expandSEToBaseSpace <- function(se, regions = NULL, + flank = 0L, seqinfo = NULL, keepAssays = getReadLevelAssayNames(se), - ignore.strand = TRUE) { + ignore.strand = TRUE, + assayNameAmbig = "mod_prob", + verbose = FALSE) { # check arguments .assertVector(x = se, type = "RangedSummarizedExperiment") + .assertScalar(x = flank, type = "numeric") + flank <- as.integer(flank) .assertVector(x = keepAssays, type = "character", validValues = getReadLevelAssayNames(se)) @@ -80,35 +102,37 @@ expandSEToBaseSpace <- function(se, ranges = IRanges(start = nnaIdx$minpos, end = nnaIdx$maxpos), seqinfo = seqinfo(se))) - }))) + })) + flank) } .assertVector(x = regions, type = "GRanges") - .assertScalar(x = ignore.strand, type = "logical", validValues = TRUE) + .assertScalar(x = ignore.strand, type = "logical", validValues = c(TRUE, FALSE)) - # note: commenting out for now, in order to allow shrinking of the - # range covered by the SE as well - # if (countOverlaps(query = range(unstrand(rowRanges(se))), - # subject = region, type = "within") == 0L) { - # cli_abort("Not all positions in {.arg se} are within {.arg region}") - # } + # subset SE to only the requested regions + se <- subsetByOverlaps(se, regions, ignore.strand = TRUE) # create base-resolution rowRanges - if (ignore.strand) { - rr <- GPos(seqnames = rep(seqnames(regions), width(regions)), - pos = unlist(lapply(seq_along(regions), function(i) { - seq(start(regions)[i], end(regions)[i]) - })), - strand = "*", - seqinfo = seqinfo) - } else { - cli_abort("{.arg ignore.strand = FALSE} is not supported yet") # nocov + rr <- GPos(seqnames = rep(seqnames(regions), width(regions)), + pos = unlist(lapply(seq_along(regions), function(i) { + seq(start(regions)[i], end(regions)[i]) + })), + strand = "*", + seqinfo = seqinfo) + if (!ignore.strand) { + nrowpre <- nrow(se) + se <- filterPositions(se, filters = "repeated.positions", + assayNameNA = NULL, + assayNameAmbig = assayNameAmbig) + nrowpost <- nrow(se) + if (nrowpost != nrowpre) { + .message("{nrowpre - nrowpost} positions filtered out due to ambiguous strand assignment", noTimer = TRUE) + } + m <- match(rowRanges(se), rr, ignore.strand = TRUE) + i <- which(!is.na(m)) + strand(rr)[m[i]] <- strand(se)[i] } - # subset SE to only the requested regions - se <- subsetByOverlaps(se, regions, ignore.strand = TRUE) - # iterate over assays - newRowIdx <- match(unstrand(rowRanges(se)), rr) + newRowIdx <- match(rowRanges(se), rr, ignore.strand = ignore.strand) assayL <- lapply(assays(se, withDimnames = FALSE)[keepAssays], function(Df) { endoapply(Df, function(naArr) { naArrBases <- NaArray( diff --git a/man/expandSEToBaseSpace.Rd b/man/expandSEToBaseSpace.Rd index 129180e..a2b86a3 100644 --- a/man/expandSEToBaseSpace.Rd +++ b/man/expandSEToBaseSpace.Rd @@ -7,9 +7,12 @@ expandSEToBaseSpace( se, regions = NULL, + flank = 0L, seqinfo = NULL, keepAssays = getReadLevelAssayNames(se), - ignore.strand = TRUE + ignore.strand = TRUE, + assayNameAmbig = "mod_prob", + verbose = FALSE ) } \arguments{ @@ -22,6 +25,11 @@ regions can be specified as character scalars (e.g. "chr1:1200-1300") that can be coerced into a \code{GRanges} object. If \code{NULL} (the default), \code{regions} is set to the range of the data in \code{se}.} +\item{flank}{Integer scalar. If \code{regions=NULL} and \code{flank} is +non-zero, \code{flank} positions will be added on both sides of each +range in the derived \code{regions}. This is useful to add a buffer of +NA values between observed values in non-overlapping regions.} + \item{seqinfo}{\code{NULL} or a \code{\link[Seqinfo]{Seqinfo}} object containing information about the set of genomic sequences (chromosomes). Alternatively, a named numeric vector with genomic sequence names and @@ -33,11 +41,25 @@ expand to base space. Only these assays will be present in the returned object.} \item{ignore.strand}{A logical scalar defining whether to ignore the strand -of the \code{rowRanges(se)}.} +of the \code{rowRanges(se)}. If \code{FALSE}, \code{se} is first +filtered (using \code{filterPositions(..., filters = "repeated.positions")}) +to resolve any positions represented on both strands.} + +\item{assayNameAmbig}{A character scalar indicating the assay to use to +decide which row to retain if multiple rows represent the same genomic +position (on different strands). The row with the largest row sum in +this assay is retained. Only used if \code{ignore.strand} is \code{FALSE}.} + +\item{verbose}{Logical scalar indicating whether to print out progress +messages.} } \value{ -A single-base resolution \code{\link[SummarizedExperiment]{RangedSummarizedExperiment}} -corresponding to \code{se}. +A single-base resolution +\code{\link[SummarizedExperiment]{RangedSummarizedExperiment}} +corresponding to \code{se}. If \code{ignore.strand} is \code{TRUE}, +the return object is unstranded. If \code{ignore.strand} is \code{FALSE}, +all positions present in the input object will retain the assigned +strand, while all other positions will be unstranded. } \description{ Expand the rows of a \code{RangedSummarizedExperiment} to single base resolution. diff --git a/tests/testthat/test-expandSEToBaseSpace.R b/tests/testthat/test-expandSEToBaseSpace.R index 0c35ee2..3232463 100644 --- a/tests/testthat/test-expandSEToBaseSpace.R +++ b/tests/testthat/test-expandSEToBaseSpace.R @@ -16,18 +16,20 @@ test_that("expandSEToBaseSpace works", { se2 <- subsetReads(se, reads = unlist(lapply(assay(se1, "mod_prob"), colnames)), invert = TRUE) Seqinfo::seqlevels(SummarizedExperiment::rowRanges(se2)) <- "chr2" - suppressWarnings(gpos <- c(rowRanges(se1), rowRanges(se2))) + suppressWarnings(gpos <- c(SummarizedExperiment::rowRanges(se1), + SummarizedExperiment::rowRanges(se2))) modmat <- make_zero_col_DFrame(nrow = length(gpos)) for (nm in colnames(se)) { x1 <- assay(se1, "mod_prob")[[nm]] x2 <- assay(se2, "mod_prob")[[nm]] + colnames(x2) <- paste0("n-", colnames(x2)) namat <- NaArray(dim = c(length(gpos), ncol(x1) + ncol(x2)), dimnames = list(NULL, c(colnames(x1), colnames(x2))), type = "double") - i1 <- match(rowRanges(se1), gpos) + i1 <- match(SummarizedExperiment::rowRanges(se1), gpos) j1 <- seq.int(ncol(x1)) namat[cbind(rep(i1, length(j1)), rep(j1, each = length(i1)))] <- c(as.matrix(x1)) - i2 <- match(rowRanges(se2), gpos) + i2 <- match(SummarizedExperiment::rowRanges(se2), gpos) j2 <- (ncol(x1) + seq.int(ncol(x2))) namat[cbind(rep(i2, length(j2)), rep(j2, each = length(i2)))] <- c(as.matrix(x2)) modmat[[nm]] <- namat @@ -39,13 +41,43 @@ test_that("expandSEToBaseSpace works", { rowRanges = gpos ) + # repeat some positions + se1 <- se + se2 <- se1[2:5, ] + rr <- SummarizedExperiment::rowRanges(se2) + strand(rr) <- ifelse(strand(rr) == "+", "-", "+") + SummarizedExperiment::rowRanges(se2) <- rr + suppressWarnings(gpos <- c(SummarizedExperiment::rowRanges(se1), + SummarizedExperiment::rowRanges(se2))) + gpos <- gpos[order(start(gpos))] + modmat <- make_zero_col_DFrame(nrow = length(gpos)) + for (nm in colnames(se)) { + x1 <- assay(se1, "mod_prob")[[nm]] + x2 <- assay(se2, "mod_prob")[[nm]] / 2 + colnames(x2) <- paste0("n-", colnames(x2)) + namat <- NaArray(dim = c(length(gpos), ncol(x1) + ncol(x2)), + dimnames = list(NULL, c(colnames(x1), colnames(x2))), + type = "double") + i1 <- match(SummarizedExperiment::rowRanges(se1), gpos) + j1 <- seq.int(ncol(x1)) + namat[cbind(rep(i1, length(j1)), rep(j1, each = length(i1)))] <- c(as.matrix(x1)) + i2 <- match(SummarizedExperiment::rowRanges(se2), gpos) + j2 <- (ncol(x1) + seq.int(ncol(x2))) + namat[cbind(rep(i2, length(j2)), rep(j2, each = length(i2)))] <- c(as.matrix(x2)) + modmat[[nm]] <- namat + } + seDupPos <- SummarizedExperiment( + assays = list(mod_prob = modmat), + colData = colData(se), + metadata = metadata(se), + rowRanges = gpos + ) + # expected fails - # expect_error(expandSEToBaseSpace(se = se, regions = "chr1:6940100-6941000"), - # "Not all positions in .se. are within .region.") expect_error(expandSEToBaseSpace(se = se, keepAssays = "doesnt-exist"), "All values in .keepAssays. must be one of: mod_prob") - expect_error(expandSEToBaseSpace(se = se, ignore.strand = FALSE), - "All values in .ignore.strand. must be one of: TRUE") + expect_error(expandSEToBaseSpace(se = se, flank = c(1, 2)), + ".flank. must have length 1") # empty input seExp0 <- expandSEToBaseSpace(se = se[numeric(0), ]) @@ -55,8 +87,8 @@ test_that("expandSEToBaseSpace works", { expect_identical(lapply(SummarizedExperiment::assay(se, "mod_prob"), colnames), lapply(SummarizedExperiment::assay(seExp0, "mod_prob"), colnames)) - # expected return value - seExp1 <- expandSEToBaseSpace(se = se) + # expected return value - ignore.strand=TRUE + seExp1 <- expandSEToBaseSpace(se = se, ignore.strand = TRUE) expect_s4_class(seExp1, "RangedSummarizedExperiment") expect_identical(dim(seExp1), c(width(range(SummarizedExperiment::rowRanges(se), @@ -74,8 +106,8 @@ test_that("expandSEToBaseSpace works", { valsExp <- SparseArray::nnavals(assayExp) expect_identical(valsOrig, valsExp) - posOrig <- GenomicRanges::start(se)[SparseArray::nnawhich(assayOrig)] - posExp <- GenomicRanges::start(seExp1)[SparseArray::nnawhich(assayExp)] + posOrig <- GenomicRanges::start(se)[SparseArray::nnawhich(assayOrig, arr.ind = TRUE)[, 1]] + posExp <- GenomicRanges::start(seExp1)[SparseArray::nnawhich(assayExp, arr.ind = TRUE)[, 1]] expect_identical(posOrig, posExp) } expect_identical(as.character(unique(GenomicRanges::strand(seExp1))), "*") @@ -83,7 +115,7 @@ test_that("expandSEToBaseSpace works", { # expanded region reg2 <- GenomicRanges::GRanges("chr1", IRanges::IRanges(6939000, 6942000)) - seExp2 <- expandSEToBaseSpace(se = se, regions = reg2) + seExp2 <- expandSEToBaseSpace(se = se, regions = reg2, ignore.strand = TRUE) expect_s4_class(seExp2, "RangedSummarizedExperiment") expect_identical(dim(seExp2), c(width(reg2), ncol(se))) expect_identical(SummarizedExperiment::assayNames(seExp2), @@ -98,14 +130,16 @@ test_that("expandSEToBaseSpace works", { valsExp <- SparseArray::nnavals(assayExp) expect_identical(valsOrig, valsExp) - posOrig <- GenomicRanges::start(se)[SparseArray::nnawhich(assayOrig)] - posExp <- GenomicRanges::start(seExp2)[SparseArray::nnawhich(assayExp)] + posOrig <- GenomicRanges::start(se)[SparseArray::nnawhich(assayOrig, arr.ind = TRUE)[, 1]] + posExp <- GenomicRanges::start(seExp2)[SparseArray::nnawhich(assayExp, arr.ind = TRUE)[, 1]] expect_identical(posOrig, posExp) } expect_identical(as.character(unique(GenomicRanges::strand(seExp2))), "*") # multiple subregions - seTwoChromExp <- expandSEToBaseSpace(se = seTwoChrom, regions = "chr1:6940000-6941000") + seTwoChromExp <- expandSEToBaseSpace(se = seTwoChrom, + regions = "chr1:6940000-6941000", + ignore.strand = TRUE) expect_s4_class(seTwoChromExp, "RangedSummarizedExperiment") expect_identical(dim(seTwoChromExp), c(1001L, ncol(seTwoChrom))) expect_identical(SummarizedExperiment::assayNames(seTwoChromExp), @@ -120,8 +154,8 @@ test_that("expandSEToBaseSpace works", { valsExp <- SparseArray::nnavals(assayExp) expect_identical(valsOrig, valsExp) - posOrig <- GenomicRanges::start(seTwoChrom[seqnames(seTwoChrom) == "chr1", ])[SparseArray::nnawhich(assayOrig)] - posExp <- GenomicRanges::start(seTwoChromExp)[SparseArray::nnawhich(assayExp)] + posOrig <- GenomicRanges::start(seTwoChrom[seqnames(seTwoChrom) == "chr1", ])[SparseArray::nnawhich(assayOrig, arr.ind = TRUE)[, 1]] + posExp <- GenomicRanges::start(seTwoChromExp)[SparseArray::nnawhich(assayExp, arr.ind = TRUE)[, 1]] expect_identical(posOrig, posExp) } expect_identical(as.character(unique(GenomicRanges::strand(seTwoChromExp))), "*") @@ -129,7 +163,8 @@ test_that("expandSEToBaseSpace works", { # multiple chromosomes reg3 <- GRanges(seqnames = c("chr1", "chr2"), ranges = IRanges(start = c(6940500, 6940200), width = 300)) - seTwoChromExp <- expandSEToBaseSpace(se = seTwoChrom, regions = reg3) + seTwoChromExp <- expandSEToBaseSpace(se = seTwoChrom, regions = reg3, + ignore.strand = TRUE) expect_s4_class(seTwoChromExp, "RangedSummarizedExperiment") expect_identical(dim(seTwoChromExp), c(sum(width(reg3)), ncol(seTwoChrom))) expect_identical(SummarizedExperiment::assayNames(seTwoChromExp), @@ -144,9 +179,87 @@ test_that("expandSEToBaseSpace works", { valsExp <- SparseArray::nnavals(assayExp) expect_identical(valsOrig, valsExp) - posOrig <- GenomicRanges::start(subsetByOverlaps(seTwoChrom, reg3))[SparseArray::nnawhich(assayOrig)] - posExp <- GenomicRanges::start(seTwoChromExp)[SparseArray::nnawhich(assayExp)] + posOrig <- GenomicRanges::start(subsetByOverlaps(seTwoChrom, reg3))[SparseArray::nnawhich(assayOrig, arr.ind = TRUE)[, 1]] + posExp <- GenomicRanges::start(seTwoChromExp)[SparseArray::nnawhich(assayExp, arr.ind = TRUE)[, 1]] expect_identical(posOrig, posExp) } expect_identical(as.character(unique(GenomicRanges::strand(seTwoChromExp))), "*") + + # multiple chromosomes with flank != 0 + seTwoChromExp0 <- expandSEToBaseSpace(se = seTwoChrom, flank = 0L, + ignore.strand = TRUE) + seTwoChromExp10 <- expandSEToBaseSpace(se = seTwoChrom, flank = 10L, + ignore.strand = TRUE) + expect_identical(nrow(seTwoChromExp10), nrow(seTwoChromExp0) + 40L) + expect_identical(lapply(assay(seTwoChrom, "mod_prob"), colnames), + lapply(assay(seTwoChromExp10, "mod_prob"), colnames)) + for (s in colnames(seTwoChrom)) { + assayExp0 <- assay(seTwoChromExp0, "mod_prob")[[s]] + assayExp10 <- assay(seTwoChromExp10, "mod_prob")[[s]] + + valsExp0 <- SparseArray::nnavals(assayExp0) + valsExp10 <- SparseArray::nnavals(assayExp10) + expect_identical(valsExp0, valsExp10) + + posExp0 <- GenomicRanges::start(seTwoChromExp0)[SparseArray::nnawhich(assayExp0, arr.ind = TRUE)[, 1]] + posExp10 <- GenomicRanges::start(seTwoChromExp10)[SparseArray::nnawhich(assayExp10, arr.ind = TRUE)[, 1]] + expect_identical(posExp0, posExp10) + } + expect_identical(as.character(unique(GenomicRanges::strand(seTwoChromExp))), "*") + + # multiple chromosomes - ignore.strand=FALSE + reg3 <- GRanges(seqnames = c("chr1", "chr2"), + ranges = IRanges(start = c(6940500, 6940200), width = 300)) + seTwoChromExp <- expandSEToBaseSpace(se = seTwoChrom, regions = reg3, + ignore.strand = FALSE) + expect_s4_class(seTwoChromExp, "RangedSummarizedExperiment") + expect_identical(dim(seTwoChromExp), c(sum(width(reg3)), ncol(seTwoChrom))) + expect_identical(SummarizedExperiment::assayNames(seTwoChromExp), + getReadLevelAssayNames(seTwoChrom)) + expect_identical(lapply(assay(seTwoChrom, "mod_prob"), colnames), + lapply(assay(seTwoChromExp, "mod_prob"), colnames)) + for (s in colnames(seTwoChrom)) { + assayOrig <- assay(subsetByOverlaps(seTwoChrom, reg3), "mod_prob")[[s]] + assayExp <- assay(seTwoChromExp, "mod_prob")[[s]] + + valsOrig <- SparseArray::nnavals(assayOrig) + valsExp <- SparseArray::nnavals(assayExp) + expect_identical(valsOrig, valsExp) + + posOrig <- GenomicRanges::start(subsetByOverlaps(seTwoChrom, reg3))[SparseArray::nnawhich(assayOrig, arr.ind = TRUE)[, 1]] + posExp <- GenomicRanges::start(seTwoChromExp)[SparseArray::nnawhich(assayExp, arr.ind = TRUE)[, 1]] + expect_identical(posOrig, posExp) + } + expect_identical(GenomicRanges::strand(subsetByOverlaps(seTwoChrom, seTwoChromExp)), + GenomicRanges::strand(seTwoChromExp)[GenomicRanges::strand(seTwoChromExp) != "*"]) + + # duplicated positions - ignore.strand=FALSE + expect_message( + seDupPosExp <- expandSEToBaseSpace(se = seDupPos, + ignore.strand = FALSE, verbose = TRUE), + "4 positions filtered out") + expect_s4_class(seDupPosExp, "RangedSummarizedExperiment") + expect_identical(dim(seDupPosExp), + c(width(range(SummarizedExperiment::rowRanges(seDupPos), + ignore.strand = TRUE)), + ncol(seDupPos))) + expect_identical(SummarizedExperiment::assayNames(seDupPosExp), + getReadLevelAssayNames(seDupPos)) + expect_identical(lapply(SummarizedExperiment::assay(seDupPos, "mod_prob"), colnames), + lapply(SummarizedExperiment::assay(seDupPosExp, "mod_prob"), colnames)) + for (s in colnames(seDupPos)) { + assayOrig <- SummarizedExperiment::assay(seDupPos, "mod_prob")[[s]] + assayOrig <- assayOrig[, seq_len(ncol(assayOrig) / 2)] + assayExp <- SummarizedExperiment::assay(seDupPosExp, "mod_prob")[[s]] + + valsOrig <- SparseArray::nnavals(assayOrig) + valsExp <- SparseArray::nnavals(assayExp) + expect_identical(valsOrig, valsExp) + + posOrig <- GenomicRanges::start(seDupPos)[SparseArray::nnawhich(assayOrig, arr.ind = TRUE)[, 1]] + posExp <- GenomicRanges::start(seDupPosExp)[SparseArray::nnawhich(assayExp, arr.ind = TRUE)[, 1]] + expect_identical(posOrig, posExp) + } + expect_identical(GenomicRanges::strand(subsetByOverlaps(seDupPos, seDupPosExp)), + GenomicRanges::strand(seDupPosExp)[GenomicRanges::strand(seDupPosExp) != "*"]) })