R/COCOA.R
In databio/PCRSA: Coordinate Covariation Analysis
Defines functions
weightedAvePerRegion
averagePerRegion
BSBinAggregate
.makeSymmetric
makeSymmetric
getMetaRegionProfile
createCorFeatureMat
aggregateSignalGRList
.formatResults
aggregateSignal
Documented in
aggregateSignal
aggregateSignalGRList
getMetaRegionProfile
# PACKAGE DOCUMENTATION
#' Coordinate Covariation Analysis (COCOA)
#'
#'
#' COCOA is a method for understanding epigenetic variation among samples.
#' COCOA can be used with epigenetic data that includes
#' genomic coordinates and an epigenetic signal,
#' such as DNA methylation and chromatin accessibility
#' data.
#' To describe the method on a high level, COCOA quantifies
#' inter-sample variation with either a supervised or unsupervised
#' technique then uses a database of "region sets" to
#' annotate the variation among samples. A region set is a set of
#' genomic regions that share a biological annotation,
#' for instance transcription factor (TF) binding regions,
#' histone modification regions, or open chromatin regions.
#' COCOA can identify region sets that are associated with
#' epigenetic variation between samples and
#' increase understanding of variation in your data.
#
# In contrast to some other common techniques, COCOA offers both
# supervised (known groups/phenotype) and unsupervised (no known groups/
# phenotype) analyses. Also, COCOA focuses on continuous variation
# between samples instead of having cutoffs. Because of this, COCOA can
# be used as a complementary method alongside "differential" methods
# that find discrete differences between groups of samples and
# it can also be used in situations where there are no groups.
# COCOA can identify biologically meaningful
# sources of variation between samples and
#'
#' @docType package
#' @name COCOA
#' @author John Lawson
#' @author Nathan Sheffield
#'
#' @references \url{http://github.com/databio}
#' @importFrom ggplot2 ggplot aes facet_wrap geom_boxplot geom_jitter geom_line
#' theme_classic xlab ylab geom_hline ylim scale_color_discrete
#' scale_x_discrete scale_fill_brewer scale_color_manual
#' scale_color_brewer theme element_line element_text geom_point
#' @importFrom ComplexHeatmap Heatmap draw
#' @import BiocGenerics S4Vectors IRanges GenomicRanges simpleCache
#' @importFrom data.table ":=" setDT data.table setkey fread setnames
#' setcolorder rbindlist setattr setorder copy is.data.table
#' setorderv as.data.table
#' @importFrom Biobase sampleNames
#' @importFrom stats lm coefficients poly wilcox.test ecdf pgamma p.adjust
#' @importFrom methods is
#' @importFrom MIRA binRegion
#' @importFrom tidyr gather
#' @importFrom grid grid.newpage grid.grabExpr grid.draw popViewport
#' pushViewport unit viewport
#' @importFrom grDevices dev.off
#' @importFrom methods hasArg
#' @importFrom fitdistrplus fitdist
#' @importFrom simpleCache simpleCache
NULL
# @importFrom ppcor pcor.test
# now package lists GenomicRanges in "Depends" instead of "Imports" in
# DESCRIPTION, still import package with @import though
# @importFrom GenomicRanges GRanges GRangesList elementMetadata strand
# seqnames granges
# Because of some issues,
# (see here: http://stackoverflow.com/questions/9439256/)
# I have to register stuff used in data.table as non-standard evaluation,
# in order to pass some R CMD check NOTES.
if (getRversion() >= "2.15.1") {
utils::globalVariables(c(
".", "..calcCols", "bin", "binID", "chr", "id", "colsToAnnotate",
"coordinateDT",
"coverage", "Group", "pOlap", "regionGroupID", "regionID", "theme",
"meanRegionSize", "regionSetCoverage", "rowIndex", "rsIndex",
"rsRegionID", "totalRegionNumber", "signalCoverage", ".SD",
"sumProportionOverlap"))
}
#########################################################################
#' Score a region set using feature contribution scores
#'
#' First, this function identifies which epigenetic features
#' overlap the region set.
#' Then the region set is scored using the feature contribution scores
#' (`signal` input)
#' according to the `scoringMetric` parameter.
#'
#' @template signal
#' @template signalCoord
#' @template signalCoordType
#' @templateVar refGenomeWarning TRUE
#' @template regionSet
#' @template signalCol
#' @template scoringMetric
#' @template verbose
# Useful when using
# aggregateSignal with 'apply' to do many region sets at a time.
# @param wilcox.conf.int logical. Only applies when using "rankSum" scoring
# method. returns a 95% confidence interval from the Wilcoxon rank sum test
# instead of p value.
#' @template absVal
#' @template rsOL
#' @param pOlap Numeric vector. Only used if rsOL is given and scoringMetric
#' is "proportionWeightedMean". This vector should contain the proportion of
#' each regionSet region that is overlapped by a signalCoord region. The
#' order of pOlap should be the same as the overlaps in rsOL.
#' @template returnCovInfo
#' @template checkInput
#' @return a data.frame with various columns:
#' One column for each item of signalCol, representing scores for the region set for each signalCol.
#' Columns: signalCoverage (cytosine_coverage) shows the number of epigenetic features that overlapped with the region set.
#' regionSetCoverage: indicates the number of regions from regionSet that overlapped with any of the epigenetic features.
#' totalRegionNumber: represents the number of regions in regionSet.
#' meanRegionSize: gives the average size (in base pairs) of regions in regionSet.
#'
#' For "multiBase" data using the "proportionWeightedMean" scoring metric,
#' the output will also include a "sumProportionOverlap" column.
#' This column represents the sum of proportion overlaps between
#' each signalCoord region and overlapping regionSet region,
#' providing an additional measure of regionSet coverage.
#'
#' @examples
#' data("brcaATACCoord1")
#' data("brcaATACData1")
#' data("esr1_chr1")
#' featureContributionScores <- prcomp(t(brcaATACData1))$rotation
#' rsScores <- aggregateSignal(signal=featureContributionScores,
#' signalCoord=brcaATACCoord1,
#' regionSet=esr1_chr1,
#' signalCol=c("PC1", "PC2"),
#' scoringMetric="default")
#' @export
aggregateSignal <- function(signal,
signalCoord,
regionSet,
signalCol = c("PC1", "PC2"),
signalCoordType = "default",
scoringMetric = "default",
verbose = FALSE,
absVal=TRUE,
rsOL=NULL, pOlap=NULL,
returnCovInfo=TRUE,
.checkInput=FALSE) {
################### checking inputs #################################
# if it was already checked outside this function, don't need to re-check
if (.checkInput) {
########## check that inputs are the correct class
checkConvertInputClasses(signal=signal,
signalCoord=signalCoord,
regionSet=regionSet,
signalCol=signalCol,
rsOL=rsOL)
########## check that dimensions of inputs are consistent
# length of signal coord = nrow of signal
if (length(signalCoord) != nrow(signal)) {
stop(cleanws("The number of coordinates in
signalCoord (length(signalCoord)) does not equal the number of
rows in `signal`"))
}
######### check that appropriate columns are present
# signalCol are column names of signal
if (!all(signalCol %in% colnames(signal))) {
missingCols = signalCol[!(signalCol %in% colnames(signal))]
stop(cleanws(paste0("Some signalCol are not
columns of signal: ", missingCols)))
}
######## check that scoringMetric is appropriate
if (!(scoringMetric %in% getScoringMethods("both"))) {
stop(cleanws("scoringMetric was not recognized.
Check spelling and available options."))
}
###### check that signalCoordType is appropriate
if (!(signalCoordType %in% c("default", "singleBase", "multiBase"))) {
stop(cleanws("signalCoordType not recognized.
Check spelling/capitalization."))
}
#######
# what happens if there are NAs or Inf in `signal`?
# any NAs that overlap the regionSet will cause the score to be NA
if (is(signal, "data.table")) {
naRows = apply(X = signal[, signalCol, with=FALSE, drop=FALSE],
MARGIN = 1, FUN = function(x) any(is.na(x)))
} else {
naRows = apply(X = signal[, signalCol, drop=FALSE],
MARGIN = 1, FUN = function(x) any(is.na(x)))
}
if (any(naRows)) {
signal <- signal[!naRows, ]
signalCoord <- signalCoord[!naRows]
warning("Removing rows with NA from `signal`")
}
#################################################################
# detect signalCoordType
if (signalCoordType == "default") {
# when signalCoord is a GRanges object
if (any(start(signalCoord) != end(signalCoord))) {
signalCoordType <- "multiBase"
} else {
signalCoordType <- "singleBase"
}
}
# if "default" scoring method is given, choose based on signalCoordType
if (scoringMetric == "default") {
if (signalCoordType == "singleBase") {
scoringMetric <- "regionMean"
} else if (signalCoordType == "multiBase") {
scoringMetric <- "proportionWeightedMean"
} else {
stop(cleanws("signalCoordType not recognized.
Check spelling/capitalization."))
}
}
# make sure that scoringMetric is consistent with signalCoordType
if (signalCoordType == "singleBase") {
if (!(scoringMetric %in% getScoringMethods("singleBase"))) {
stop("The scoringMetric you selected is not available for
this data's signalCoordType")
}
} else if (signalCoordType == "multiBase") {
if (!(scoringMetric %in% getScoringMethods("multiBase"))) {
stop("The scoringMetric you selected is not available for
this data's signalCoordType")
}
}
}
# extreme positive or negative values both give important information
# take absolute value or not
if (absVal) {
signal <- abs(signal) # required for later code
}
# XX copies unnecessarily:reformat into data.table with chromosome location and weight
# make sure `signal` is the correct type for further steps
# (proportionWeightedMean method requires a matrix)
if (!is(signal, "data.table") && (scoringMetric != "proportionWeightedMean")) {
signal <- as.data.table(signal)
} else if (!is(signal, "matrix") && (scoringMetric == "proportionWeightedMean")) {
signal <- as.matrix(signal)
}
# restricting to signalCol so unnecessary computations
# are not done
if (is(signal, "data.table")) {
loadingDT <- signal[, signalCol, with=FALSE]
# # naming does not work if only using one PC so add this line for that case
# setnames(loadiangDT, signalCol)
} else {
loadingDT <- signal[, signalCol, drop=FALSE]
}
###########################################################################
# scoring
# would rounding speed up aggregation?, potentially make a sparse matrix
# if a lot of entries became 0
# works for both singleBase and multiBase
if (scoringMetric == "simpleMean") {
loadAgMain <- as.data.table(regionOLMean(signalDT = loadingDT,
signalGR = signalCoord,
regionSet = regionSet,
calcCols= signalCol,
metric = "mean",
rsOL = rsOL,
returnCovInfo=returnCovInfo))
results <- .formatResults(loadAgMain, scoringMetric = scoringMetric,
regionSet=regionSet, signalCol = signalCol,
returnCovInfo=returnCovInfo)
} else if (scoringMetric == "simpleMedian") {
# scoring singleBase and multiBase both with this function for
# simpleMedian
loadAgMain <- as.data.table(regionOLMean(signalDT = loadingDT,
signalGR = signalCoord,
regionSet = regionSet,
calcCols= signalCol,
metric = "median",
rsOL=rsOL,
returnCovInfo = returnCovInfo))
results <- .formatResults(loadAgMain, scoringMetric = scoringMetric,
regionSet=regionSet, signalCol = signalCol,
returnCovInfo=returnCovInfo)
} else if (signalCoordType == "singleBase") {
# do the actual aggregation
if (scoringMetric == "regionMean") {
# specify aggregation operation
# will be done separately for each PC specified
aggrCommand <- buildJ(signalCol,
rep("mean", length(signalCol)))
# previously used BSAggregate from RGenomeUtils but now using local,
# modified copy
loadAgMain <- BSAggregate(BSDT = loadingDT,
regionsGRL = regionSet,
BSCoord = signalCoord,
jExpr = aggrCommand,
byRegionGroup = TRUE,
splitFactor = NULL,
returnOLInfo = returnCovInfo,
rsOL=rsOL)
results <- .formatResults(loadAgMain,
scoringMetric = scoringMetric,
regionSet=regionSet, signalCol = signalCol,
returnCovInfo=returnCovInfo)
} else if (scoringMetric == "regionMedian") {
aggrCommand <- buildJ(signalCol,
rep("median", length(signalCol)))
loadAgMain <- BSAggregate(BSDT = loadingDT,
regionsGRL = regionSet,
BSCoord = signalCoord,
jExpr = aggrCommand,
byRegionGroup = TRUE,
splitFactor = NULL,
returnOLInfo = returnCovInfo)
results <- .formatResults(loadAgMain,
scoringMetric = scoringMetric,
regionSet=regionSet, signalCol = signalCol,
returnCovInfo=returnCovInfo)
}
# } else if (scoringMetric == "meanDiff") {
# # if (is.null(pcLoadAv)) {
# # # calculate (should already be absolute)
# # pcLoadAv <- apply(X = loadingDT[, signalCol, with=FALSE],
# # MARGIN = 2, FUN = mean)
# # }
# loadMetrics <- signalOLMetrics(dataDT=loadingDT, regionSet=regionSet,
# signalCol = signalCol,
# metrics=c("mean", "sd"),
# alsoNonOLMet=TRUE)
# if (is.null(loadMetrics)) {
# results <- as.data.table(t(rep(NA, length(signalCol))))
# setnames(results, signalCol)
# results[, signalCoverage := 0]
# results[, regionSetCoverage := 0]
# results[, totalRegionNumber := numOfRegions]
# results[, meanRegionSize := round(mean(width(regionSet)), 1)]
# } else {
# # calculate mean difference
# # pooled standard deviation
# sdPool <- sqrt((loadMetrics$sd_OL^2 + loadMetrics$sd_nonOL^2) / 2)
#
# # mean difference
# # error if signalCoverage > (1/2) * totalCpGs
# meanDiff <- (loadMetrics$mean_OL - loadMetrics$mean_nonOL) /
# (sdPool * sqrt((1 / loadMetrics$signalCoverage) - (1 / (totalCpGs - loadMetrics$signalCoverage))))
# results <- as.data.table(t(meanDiff))
# colnames(results) <- loadMetrics$testCol
#
# # add information about degree of overlap
# results <- cbind(results, loadMetrics[1, .SD, .SDcols = c("signalCoverage", "regionSetCoverage", "totalRegionNumber", "meanRegionSize")])
# }
#
# } else if (scoringMetric == "rankSum") {
#
# # if (wilcox.conf.int) {
# # # returns confidence interval
# # wRes <- rsWilcox(dataDT = loadingDT, regionSet=regionSet,
# # signalCol = signalCol,
# # conf.int = wilcox.conf.int)
# #
# # if (is.null(wRes)) {
# # results <- as.data.table(t(rep(NA, length(signalCol) * 2)))
# # setnames(results, paste0(rep(signalCol, each=2),
# # c("_low", "_high")))
# # results[, signalCoverage := 0]
# # results[, regionSetCoverage := 0]
# # results[, totalRegionNumber := numOfRegions]
# # results[, meanRegionSize := round(mean(width(regionSet)), 1)]
# # } else {
# # results <- as.data.table(wRes)
# # }
# #
# # } else {
# # returns p value
# # one sided test since I took the absolute value of the loadings
# wRes <- rsWilcox(dataDT = loadingDT, regionSet=regionSet,
# signalCol = signalCol, alternative="greater")
#
# if (is.null(wRes)) {
# results <- as.data.table(t(rep(NA, length(signalCol))))
# setnames(results, signalCol)
# results[, signalCoverage := 0]
# results[, regionSetCoverage := 0]
# results[, totalRegionNumber := numOfRegions]
# results[, meanRegionSize := round(mean(width(regionSet)), 1)]
# } else {
# results <- as.data.table(wRes)
# }
# #}
} else {
# signalCoordType == "multiBase"
# for ATAC-seq
if (scoringMetric == "proportionWeightedMean") {
loadAgMain <- regionOLWeightedMean(signalMat = loadingDT,
signalGR = signalCoord,
regionSet = regionSet,
calcCols= signalCol,
rsOL = rsOL,
pOlap = pOlap,
returnCovInfo=returnCovInfo)
results <- .formatResults(as.data.table(loadAgMain),
scoringMetric = scoringMetric,
regionSet=regionSet, signalCol = signalCol,
returnCovInfo=returnCovInfo)
}
}
# signalOLMetrics() # make sure it works with no overlap
if (verbose) {
message(":", appendLF=FALSE)
}
return(as.data.frame(results))
}
# format results of scoring functions in aggregateSignal()
.formatResults <- function(loadAgMain, scoringMetric,
regionSet, signalCol, returnCovInfo=TRUE) {
numOfRegions <- length(regionSet)
# if no cytosines from loadings were included in regionSet, result is NA
if (is.null(loadAgMain)) {
results <- as.data.table(t(rep(NA, length(signalCol))))
setnames(results, signalCol)
if (returnCovInfo) {
results[, signalCoverage := 0]
results[, regionSetCoverage := 0]
results[, totalRegionNumber := numOfRegions]
results[, meanRegionSize := round(mean(width(regionSet)), 1)]
# this column is only added by this scoring method
if (scoringMetric == "proportionWeightedMean") {
results[, sumProportionOverlap := 0]
}
}
} else {
# regionMean, regionMedian, simpleMean for region data
results <- loadAgMain[, .SD, .SDcols = signalCol]
if (returnCovInfo) {
results[, signalCoverage := loadAgMain[, .SD, .SDcols = "signalCoverage"]]
results[, regionSetCoverage := loadAgMain[, .SD, .SDcols = "regionSetCoverage"]]
results[, totalRegionNumber := numOfRegions]
results[, meanRegionSize := round(mean(width(regionSet)), 1)]
################################
# proportionWeightedMean
# this column is only added by this scoring method
if (scoringMetric == "proportionWeightedMean") {
results[, sumProportionOverlap := loadAgMain[, .SD, .SDcols = "sumProportionOverlap"]]
}
}
##############################
# simple mean for base pair data
}
return(results)
}
#' Score many region sets
#'
#' This function will give each region set a score for each target variable
#' given by `signalCol` based on
#' the `scoringMetric` parameter. Based on these scores, you can determine
#' which region sets out of a region set database (given by `GRList`)
#' are most associated with the target variables. See the vignette "Introduction
#' to Coordinate Covariation Analysis" for help interpreting your
#' results.
#'
#' @template signal
#' @template signalCoord
#' @template signalCoordType
#' @template GRList
#' @template signalCol
#' @template scoringMetric
#' @template verbose
# @param wilcox.conf.int logical. Only applies when using "rankSum" scoring
# method. returns a 95% confidence interval from the Wilcoxon rank sum test
# instead of p value.
#' @template absVal
#' @template olList
#' @template pOlapList
#' @template returnCovInfo
#' @return Data.frame of results, one row for each region set.
#' It has the following columns:
#' one column for each item of signalCol with names given by signalCol.
#' These columns have scores for the region set for each signalCol.
#' Other columns: signalCoverage (formerly cytosine_coverage) which
#' has number of epigenetic features that overlapped at all with regionSet,
#' regionSetCoverage which has number of regions from regionSet
#' that overlapped any of the epigenetic features,
#' totalRegionNumber that has number of regions in regionSet,
#' meanRegionSize that has average size in base pairs of regions in regionSet,
#' the average is based on all regions in regionSet and not just ones that overlap.
#'
#' For "multiBase" data, if the "proportionWeightedMean" scoring metric
#' is used, then the output will also have a "sumProportionOverlap" column.
#' During this scoring method, the proportion overlap between each signalCoord
#' region and overlapping regionSet region is calculated. This column is
#' the sum of all those proportion overlaps and is another way to quantify
#' coverage of regionSet in addition to regionSetCoverage.
#'
#'
#' @examples
#' data("brcaATACCoord1")
#' data("brcaATACData1")
#' data("esr1_chr1")
#' data("nrf1_chr1")
#' featureContributionScores <- prcomp(t(brcaATACData1))$rotation
#' GRList <- GRangesList(esr1_chr1, nrf1_chr1)
#' rsScores <- aggregateSignalGRList(signal=featureContributionScores,
#' signalCoord=brcaATACCoord1,
#' GRList= GRList,
#' signalCol=c("PC1", "PC2"),
#' scoringMetric="default")
#'
#' @export
aggregateSignalGRList <- function(signal,
signalCoord,
GRList,
signalCol = c("PC1", "PC2"),
signalCoordType = "default",
scoringMetric = "default",
verbose = TRUE,
absVal=TRUE, olList=NULL,
pOlapList=NULL, returnCovInfo=TRUE) {
################### checking inputs #################################
########## check that inputs are the correct class
checkConvertInputClasses(signal=signal,
signalCoord=signalCoord,
regionSet=NULL,
signalCol = signalCol,
GRList=GRList,
olList=olList)
########## check that dimensions of inputs are consistent
# length of signal coord = nrow of signal
if (length(signalCoord) != nrow(signal)) {
stop(cleanws("The number of coordinates in
signalCoord (length(signalCoord)) does not equal the number of
rows in `signal`"))
}
######### check that appropriate columns are present
# signalCol are column names of signal
if (!all(signalCol %in% colnames(signal))) {
missingCols = signalCol[!(signalCol %in% colnames(signal))]
stop(cleanws(paste0("Some signalCol are not
columns of signal: ", missingCols)))
}
######## check that scoringMetric is appropriate
if (!(scoringMetric %in% getScoringMethods("both"))) {
stop(cleanws("scoringMetric was not recognized.
Check spelling and available options."))
}
###### check that signalCoordType is appropriate
if (!(signalCoordType %in% c("default", "singleBase", "multiBase"))) {
stop(cleanws("signalCoordType not recognized.
Check spelling/capitalization."))
}
#######
# what happens if there are NAs or Inf in `signal`?
# any NAs that overlap the regionSet will cause the score to be NA
if (is(signal, "data.table")) {
naRows = apply(X = signal[, signalCol, with=FALSE, drop=FALSE],
MARGIN = 1, FUN = function(x) any(is.na(x)))
} else {
naRows = apply(X = signal[, signalCol, drop=FALSE],
MARGIN = 1, FUN = function(x) any(is.na(x)))
}
if (any(naRows)) {
signal <- signal[!naRows, ]
signalCoord <- signalCoord[!naRows]
warning("Removing rows with NA from `signal`")
}
#################################################################
# detect signalCoordType
if (signalCoordType == "default") {
# when signalCoord is a GRanges object
if (any(start(signalCoord) != end(signalCoord))) {
signalCoordType <- "multiBase"
} else {
signalCoordType <- "singleBase"
}
}
# if "default" scoring method is given, choose based on signalCoordType
if (scoringMetric == "default") {
if (signalCoordType == "singleBase") {
scoringMetric <- "regionMean"
} else if (signalCoordType == "multiBase") {
scoringMetric <- "proportionWeightedMean"
} else {
stop(cleanws("signalCoordType not recognized.
Check spelling/capitalization."))
}
}
# convert object class outside aggregateSignal to extra prevent copying
# (one scoring method needs `signal` as a matrix though)
if (!is(signal, "data.table") && (scoringMetric != "proportionWeightedMean")) {
signal <- as.data.table(signal)
} else if (!is(signal, "matrix") && (scoringMetric == "proportionWeightedMean")) {
signal <- as.matrix(signal)
}
# take absolute value outside aggregateSignal to prevent extra copying
if (absVal) {
if (is(signal, "data.table")) {
signal[, signalCol] <- abs(signal[, signalCol, with=FALSE])
} else {
signal[, signalCol] <- abs(signal[, signalCol])
}
absVal <- FALSE
}
if (is.null(olList)) {
# apply over the list of region sets
resultsList <- lapplyAlias(GRList,
function(x) aggregateSignal(
signal = signal,
signalCoord = signalCoord,
signalCoordType = signalCoordType,
regionSet = x,
signalCol = signalCol,
scoringMetric = scoringMetric,
verbose = verbose,
absVal = absVal, pOlap=NULL,
returnCovInfo=returnCovInfo))
# resultsList <- mapply(FUN = function(x) aggregateSignal(
# signal = signal,
# signalCoord = signalCoord,
# signalCoordType = signalCoordType,
# regionSet = x,
# signalCol = signalCol,
# scoringMetric = scoringMetric,
# verbose = verbose,
# absVal = absVal, pOlap=y,
# returnCovInfo=returnCovInfo), x=GRList, y=pOlapList)
#wilcox.conf.int = wilcox.conf.int,
} else {
# # apply over the list of region sets
# resultsList <- lapplyAlias(olList,
# function(x) aggregateSignal(
# signal = signal,
# signalCoord = signalCoord,
# signalCoordType = signalCoordType,
# regionSet = NULL,
# signalCol = signalCol,
# scoringMetric = scoringMetric,
# verbose = verbose,
# absVal = absVal,
# rsOL=x, pOlap=pOlapList, returnCovInfo=returnCovInfo))
# #wilcox.conf.int = wilcox.conf.int,
if (is.null(pOlapList)) {
# mapply does not work if pOlapList is null and other argument is not null
resultsList <- lapplyAlias(olList, FUN = function(x) aggregateSignal(
signal = signal,
signalCoord = signalCoord,
signalCoordType = signalCoordType,
regionSet = NULL,
signalCol = signalCol,
scoringMetric = scoringMetric,
verbose = verbose,
absVal = absVal,
rsOL=x,
returnCovInfo=returnCovInfo))
} else {
resultsList <- mapply(FUN = function(x, y) aggregateSignal(
signal = signal,
signalCoord = signalCoord,
signalCoordType = signalCoordType,
regionSet = NULL,
signalCol = signalCol,
scoringMetric = scoringMetric,
verbose = verbose,
absVal = absVal,
rsOL=x, pOlap=y,
returnCovInfo=returnCovInfo),
x=olList, y=pOlapList)
}
}
resultsDT <- do.call(rbind, resultsList)
# # add names if they are present
# if (!is.null(names(GRList))) {
# # resultsDT[, rsName := names(GRList)]
# row.names(resultsDT) <- names(GRList)
# }
# Convert the first column to numeric for proper sorting
resultsDT[, 1] <- as.numeric(resultsDT[, 1])
# rsName <- row.names(resultsDT)
resultsDF <- as.data.frame(resultsDT) #, row.names = rsName)
# Sort the results based on the first column
resultsDF <- resultsDF[order(-resultsDF[, 1]), ]
return(resultsDF)
}
# @param dataMat columns of dataMat should be samples/patients, rows should be genomic signal
# (each row corresponds to one genomic coordinate/range)
# @param featureMat Rows should be samples, columns should be "features"
# (whatever you want to get correlation with: eg PC scores),
# all columns in featureMat will be used (subset when passing to function
# in order to not use all columns)
# @param centerDataMat logical object. Should rows in dataMat be centered based on
# their means? (subtracting row mean from each row)
# @param centerFeatureMat logical. Should columns in featureMat be centered based
# on their means? (subtract column mean from each column)
# @param testType character object. Can be "cor" (Pearson correlation),
# "spearmanCor (Spearman correlation)
# "pcor" (partial correlation), "cov" (covariance (Pearson)),
# @param covariate
#
# If a row in dataMat has 0 stand. deviation, correlation will be set to 0
# instead of NA as would be done by cor()
#
# @return returns a matrix where rows are the genomic signal (eg a CpG or region) and
# columns are the columns of featureMat
# @examples dataMat = matrix(rnorm(50), 5, 10)
# featureMat = matrix(rnorm(20), 10, 2)
createCorFeatureMat <- function(dataMat, featureMat,
centerDataMat=TRUE, centerFeatureMat = TRUE,
testType="cor", covariate=NULL) {
featureMat <- as.matrix(featureMat)
featureNames <- colnames(featureMat)
nFeatures <- ncol(featureMat)
nDataDims <- nrow(dataMat)
if (centerDataMat) {
cpgMeans <- rowMeans(dataMat, na.rm = TRUE)
# centering before calculating correlation
dataMat <- apply(X = dataMat, MARGIN = 2, function(x) x - cpgMeans)
}
if (centerFeatureMat) {
featureMeans <- colMeans(featureMat, na.rm = TRUE)
# centering before calculating correlation(also, t() converts to matrix)
featureMat <- t(apply(X = t(featureMat), MARGIN = 2, function(x) x - featureMeans))
if (dim(featureMat)[1] == 1) {
featureMat <- t(featureMat)
}
}
# avoid this copy and/or delay transpose until after calculating correlation?
dataMat <- as.data.frame(t(dataMat))
if (testType == "cor") {
# create feature correlation matrix with PCs (rows: features/CpGs, columns:PCs)
# how much do features correlate with each PC?
# put epigenetic data first in cor()
featurePCCor <- cor(dataMat, featureMat, use="pairwise.complete.obs", method="pearson")
} else if (testType == "spearmanCor") {
# xtfrm(x) ranking
featurePCCor <- cor(dataMat, featureMat, use="pairwise.complete.obs", method="spearman")
# } else if (testType == "pcor") {
# # partial correlation (account for covariates), ppcor package
#
# featurePCCor <- apply(X = featureMat, MARGIN = 2, function(y) apply(X = dataMat, 2,
# FUN = function(x) pcor.test(x = x, y=y,
# z=covariate,
# method="pearson")$estimate))
#
} else if (testType == "cov") {
featurePCCor <- cov(dataMat, featureMat, use="pairwise.complete.obs")
} else {
stop("invalid testType")
}
# if standard deviation of the data was zero, NA will be produced
# set to 0 because no standard deviation means no correlation with attribute of interest
featurePCCor[is.na(featurePCCor)] <- 0
colnames(featurePCCor) <- featureNames
return(featurePCCor)
# corLoadRatio <- signal[, signalCol] / featurePCCor
# hist(corLoadRatio[, "PC10"])
}
#' Create a "meta-region" profile
#'
#' This profile can show enrichment
#' of genomic signals with high feature contribution scores
#' in the region set but not in the
#' surrounding genome, suggesting that variation is linked specifically
#' to that region set.
#'
#' All regions in a given region set
#' are combined into a single aggregate profile. Regions in `regionSet`
#' should be expanded on each side to include a wider area of the genome around
#' the regions of interest (see example and vignettes).
#' To make the profile, first we optionally take
#' the absolute value of `signal` (`absVal` parameter).
#' Then each expanded regionSet region is
#' split into `binNum` bins. The corresponding
#' bins from each region (e.g. all bin1's, all bin2's, etc.) are grouped.
#' All overlapping values from `signal` are aggregated in each bin group
#' according to the `aggrMethod` parameter to get a meta-region profile.
#' Since DNA strand information is not considered,
#' the profile is averaged symmetrically around the center.
#' A peak in the middle of this profile suggests that variability is specific
#' to the region set of interest and is not a product of the surrounding genome.
#' A region set can still be significant even if it does not have a peak.
#' For example, some histone modification region sets may be in large genomic blocks
#' and not show a peak, despite having variation across samples.
#'
#' @template signal
#' @template signalCoord
#' @template regionSet
#' @template signalCol
#' @template signalCoordType
#' @param binNum Number of bins to split each region into when
#' making the aggregate profile. More bins will
#' give a higher resolution but perhaps more noisy profile.
#' @template verbose
#' @templateVar usesAggrMethod TRUE
#' @template scoringMetric
#' @template absVal
#' @return A data.frame with the binned meta-region profile,
#' one row per bin. columns: binID and one column for each target variable
#' in signalCol. The function will return NULL if there
#' is no overlap between signalCoord and any of the bin groups that come
#' from regionSet (e.g. none of the bin1's overlapped signalCoord,
#' NULL returned).
#'
#' @examples
#' data("brcaATACCoord1")
#' data("brcaATACData1")
#' data("esr1_chr1")
#' featureContributionScores <- prcomp(t(brcaATACData1))$rotation
#' esr1_chr1_expanded <- resize(esr1_chr1, 12000, fix="center")
#' mrProfile <- getMetaRegionProfile(signal=featureContributionScores,
#' signalCoord=brcaATACCoord1,
#' regionSet=esr1_chr1_expanded,
#' signalCol=c("PC1", "PC2"),
#' binNum=21)
#' @export
getMetaRegionProfile <- function(signal, signalCoord, regionSet,
signalCol = c("PC1", "PC2"),
signalCoordType = "default",
binNum = 21,
verbose=TRUE,
aggrMethod = "default", absVal=TRUE) {
################### checking inputs #################################
########## check that inputs are the correct class, convert
checkConvertInputClasses(signal=signal,
signalCoord=signalCoord,
regionSet=regionSet,
signalCol = signalCol)
########## check that dimensions of inputs are consistent
# length of signal coord = nrow of signal
if (length(signalCoord) != nrow(signal)) {
stop(cleanws("The number of coordinates in
signalCoord (length(signalCoord)) does not equal the number of
rows in `signal`"))
}
######### check that appropriate columns are present
# signalCol are column names of `signal`
if (!all(signalCol %in% colnames(signal))) {
missingCols = signalCol[!(signalCol %in% colnames(signal))]
stop(cleanws(paste0("Some signalCol are not
columns of signal: ", missingCols)))
}
######## check that aggregation method is appropriate
if (!(aggrMethod %in% getScoringMethods("metaRegionProfile"))) {
stop(cleanws("scoringMetric was not recognized.
Check spelling and available options."))
}
###### check that signalCoordType is appropriate
if (!(signalCoordType %in% c("default", "singleBase", "multiBase"))) {
stop(cleanws("signalCoordType not recognized.
Check spelling/capitalization."))
}
#######
# what happens if there are NAs or Inf in `signal`?
#################################################################
# detect signalCoordType
if (signalCoordType == "default") {
# when signalCoord is a GRanges object
if (any(start(signalCoord) != end(signalCoord))) {
signalCoordType <- "multiBase"
} else {
signalCoordType <- "singleBase"
}
}
# if "default" aggregation method is given, choose based on signalCoordType
if (aggrMethod == "default") {
if (signalCoordType == "singleBase") {
aggrMethod <- "regionMean"
} else if (signalCoordType == "multiBase") {
aggrMethod <- "proportionWeightedMean"
} else {
stop(cleanws("signalCoordType not recognized.
Check spelling/capitalization."))
}
}
##################################################################
# take absolute value or not
if (absVal) {
loadingDT <- as.data.table(abs(signal))
} else {
loadingDT <- as.data.table(signal)
}
GRDT <- grToDt(regionSet)
loadProf <- BSBinAggregate(BSDT = loadingDT,
rangeDT = GRDT,
binCount = binNum,
BSCoord = signalCoord,
byRegionGroup = TRUE,
splitFactor = NULL,
signalCol = signalCol,
aggrMethod = aggrMethod)
# if loadProf is NULL, return NULL from function, otherwise make symmetrical
# it will be NULL when there was no overlap between data and any of the bins
if (!is.null(loadProf)) {
loadProf <- makeSymmetric(loadProf)
loadProf[, regionGroupID := seq_len(binNum)][]
setnames(loadProf, old = "regionGroupID", new="binID")
loadProf <- as.data.frame(loadProf)
} else {
warning("Insufficient overlap between regionSet and signalCoord")
}
return(loadProf)
}
makeSymmetric <- function(prof) {
symProf <- apply(prof, 2, .makeSymmetric)
symProf <- as.data.table(symProf)
return(symProf)
}
.makeSymmetric <- function(vec) {
symVec <- (vec + rev(vec)) / 2
return(symVec)
}
# Produced originally for binning Ewing RRBS data across various region sets
#
# @param BSDT A data.table. For COCOA, a data.table of loading values
# with the PCs to be annotated. One column for the signal of each
# target variable
# and also has columns with the coordinates for the epigenetic features:
# chr (chromosome) and start column (also possibly end)
# @param rangeDT A data.table with the sets of regions to be binned,
# with columns named start, end
# @param binCount Number of bins across the region
# @param BSCoord GRanges. Coordinates for BSDT. If NULL, then "chr" and "start"
# columns must be in BSDT.
# @param byRegionGroup Pass along to binCount (see ?binCount)
# @template signalCol
# @param verbose A "logical" object. Whether progress
# of the function should be shown, one
# bar indicates the region set is completed.
# useful when using BSBinAggregate with 'apply' to do many
# region sets at a time.
# @param aggrMethod see ?getMetaRegionProfile()
BSBinAggregate <- function(BSDT, rangeDT, binCount,
BSCoord=NULL,
byRegionGroup = TRUE,
splitFactor = NULL,
signalCol,
verbose = FALSE,
aggrMethod) {
if (!is(rangeDT, "data.table")) {
stop("rangeDT must be a data.table")
}
if (is.null(BSCoord)) {
BSCoord <- BSdtToGRanges(list(BSDT))
}
seqnamesColName <- "seqnames" # default column name
if (! "seqnames" %in% colnames(rangeDT)) {
if ("chr" %in% colnames(rangeDT)) {
# message("seqnames column name set to: chr")
seqnamesColName <- "chr"
} else {
# Got neither.
stop("rangeDT must have a seqnames column")
}
}
# message("Binning...")
binnedDT <- rangeDT[, MIRA::binRegion(start, end,
binCount, get(seqnamesColName))]
# output is a list of GRanges objects, does not play well with vapply
# one GRanges object for each bin, containing a segment of each original rangeDT region
binnedGR <- sapply(split(binnedDT, binnedDT$binID), dtToGr)
# message("Aggregating...")
if (aggrMethod == "proportionWeightedMean") {
binMeansList <- lapply(X = binnedGR,
FUN = function(x) regionOLWeightedMean(signalMat = BSDT,
# signalGR = dtToGr(BSDT[, .(chr, start, end)]),
signalGR = BSCoord,
regionSet = x,
calcCols = signalCol))
# any bins that had no overlap with data will be NULL
# if any bins had no data, return NULL
if (any(vapply(X = binMeansList, FUN = is.null, FUN.VALUE = TRUE))) {
return(NULL)
}
# rbindlist of all NULL's will still return an (empty) data.table
# otherwise any NULL items will just be skipped and other rows will
# be concatenated
binnedBSDT <- rbindlist(binMeansList)
regionGroupID = 1:length(binMeansList)
# regionGroupID = regionGroupID[!vapply(X = binMeansList, FUN = is.null, FUN.VALUE = TRUE)]
binnedBSDT[, regionGroupID := regionGroupID]
} else if (aggrMethod == "simpleMean") {
binMeansList <- lapply(X = binnedGR,
FUN = function(x) regionOLMean(signalDT = BSDT,
# signalGR = dtToGr(BSDT[, .(chr, start, end)]),
signalGR = BSCoord,
regionSet = x,
calcCols = signalCol))
# any bins that had no overlap with data will be NULL
# if any bins had no data, return NULL
if (any(vapply(X = binMeansList, FUN = is.null, FUN.VALUE = TRUE))) {
return(NULL)
}
binnedBSDT <- rbindlist(binMeansList)
regionGroupID = 1:length(binMeansList)
# regionGroupID = regionGroupID[!vapply(X = binMeansList, FUN = is.null, FUN.VALUE = TRUE)]
binnedBSDT[, regionGroupID := regionGroupID]
} else if (aggrMethod == "regionMean") { # aggrMethod == "regionMean"
# what is output if a region set has no overlap?
binnedBSDT <- BSAggregate(BSDT,
regionsGRL=GRangesList(binnedGR),
BSCoord = BSCoord,
jExpr=buildJ(signalCol,
rep("mean", length(signalCol))),
byRegionGroup = byRegionGroup,
splitFactor = splitFactor)
# if any bins had no data
if (nrow(binnedBSDT) < binCount) {
return(NULL)
}
} else if (aggrMethod == "regionMedian") {
# what is output if a region set has no overlap?
binnedBSDT <- BSAggregate(BSDT,
regionsGRL=GRangesList(binnedGR),
BSCoord = BSCoord,
jExpr=buildJ(signalCol,
rep("median", length(signalCol))),
byRegionGroup = byRegionGroup,
splitFactor = splitFactor)
# if any bins had no data
if (nrow(binnedBSDT) < binCount) {
return(NULL)
}
}
# RGenomeUtils::BSAggregate
# # If we aren't aggregating by bin, then don't restrict to min reads!
# if (byRegionGroup) {
# binnedBSDT <- binnedBSDT[readCount > minReads,]
# }
if (verbose) {
message(".", appendLF=FALSE)
}
return(binnedBSDT)
}
# modification of BSAggregate to just return mean per region
#
# @template signal
# @template signalCoord
# @template regionSet
# Must be from the same reference genome
# as the coordinates for the actual data/samples (signalCoord).
# @template signalCol
# @param returnQuantile "logical" object. If FALSE, return region averages. If TRUE,
# for each region, return the quantile of that region's average value
# based on the distribution of individual genomic signal/feature values
# @template absVal
# @return a data.table with region coordinates and average loading
# values for each region. Has columns chr, start, end, and a column for each
# target variable in signalCol.
# Regions are not in order along the rows of the data.table.
#
# @example averagePerRegion(BSDT = BSDT, regionsGRL,
# jCommand = MIRA:::buildJ(cols = "methylProp", "mean"))
# Devel note: I could add a column for how many cytosines are in each region
averagePerRegion <- function(signal,
signalCoord,
regionSet,
signalCol = c("PC1", "PC2"),
returnQuantile = FALSE,
absVal=TRUE) {
################### checking inputs #################################
########## check that inputs are the correct class and converts
checkConvertInputClasses(signal=signal,
signalCoord=signalCoord,
regionSet=regionSet,
signalCol = signalCol)
########## check that dimensions of inputs are consistent
# length of signal coord = nrow of signal
if (length(signalCoord) != nrow(signal)) {
stop(cleanws("The number of coordinates in
signalCoord (length(signalCoord)) does not equal the number of
rows in `signal`"))
}
######### check that appropriate columns are present
# signalCol are column names of signal
if (!all(signalCol %in% colnames(signal))) {
missingCols = signalCol[!(signalCol %in% colnames(signal))]
stop(cleanws(paste0("Some signalCol are not
columns of signal: ", missingCols)))
}
#######
# what happens if there are NAs or Inf in `signal`?
# any NAs that overlap the regionSet will cause the score to be NA
if (is(signal, "data.table")) {
naRows = apply(X = signal[, signalCol, with=FALSE, drop=FALSE],
MARGIN = 1, FUN = function(x) any(is.na(x)))
} else {
naRows = apply(X = signal[, signalCol, drop=FALSE],
MARGIN = 1, FUN = function(x) any(is.na(x)))
}
if (any(naRows)) {
signal <- signal[!naRows, ]
signalCoord <- signalCoord[!naRows]
warning("Removing rows with NA from `signal`")
}
#################################################################
#################################################################
# determine whether coordinates are single base or a range
# if (!("end" %in% colnames(coordinateDT))) {
# dataCoordType <- "singleBase"
# } else {
if (any(start(signalCoord) != end(signalCoord))) {
dataCoordType <- "multiBase"
} else {
dataCoordType <- "singleBase"
}
# }
# take absolute value or not
if (absVal) {
signalDT <- as.data.table(abs(signal))
} else {
signalDT <- as.data.table(signal)
}
# use different function for single base data and for region data
if (dataCoordType == "singleBase") {
# linking coordinates to loading values, has columns chr start, signalCol
BSDT <- signalDT[, .SD, .SDcols = signalCol]
jExpr <- buildJ(signalCol, rep("mean", length(signalCol)))
avPerRegion <- BSAggregate(BSDT = BSDT,
regionsGRL = regionSet,
BSCoord = signalCoord,
excludeGR = NULL,
regionsGRL.length = NULL,
splitFactor = NULL,
keepCols = NULL,
sumCols = NULL,
jExpr = jExpr,
byRegionGroup = FALSE,
keep.na = FALSE,
returnSD = FALSE,
returnOLInfo = FALSE,
meanPerRegion = TRUE,
returnQuantile = returnQuantile)
} else if (dataCoordType == "multiBase") {
avPerRegion <- weightedAvePerRegion(signalDT = signalDT,
signalCoord=signalCoord,
regionSet=regionSet,
calcCols = signalCol,
returnQuantile = returnQuantile)
} else {
stop("dataCoordType not recognized.")
}
return(avPerRegion)
}
# signalDT
# average by region for region-based data (eg ATAC-seq)
# @return data.table. chr, start, end columns. One column for each calcCols
# that has the average value in each region for that col
weightedAvePerRegion <- function(signalDT,
signalCoord,
regionSet,
calcCols = c("PC1", "PC2"),
returnQuantile = FALSE) {
if (is(signalCoord, "data.frame")) {
signalCoord <- dtToGr(signalCoord)
}
hits <- findOverlaps(query = signalCoord, subject = regionSet)
# if no overlap, return NULL
if (length(hits) == 0) {
return(NULL)
}
olap <- pintersect(signalCoord[queryHits(hits)],
regionSet[subjectHits(hits)])
polap <- width(olap) / width(regionSet[subjectHits(hits)])
# get total proportion overlap per region
# aggregate polap by region
pOlapDT <- data.table(signalDT[queryHits(hits), calcCols, with=FALSE],
rsRegionID = subjectHits(hits),
pOlap = polap)
pOlapByRegionDT <- pOlapDT[, .(regionPOlap = sum(pOlap)), by=rsRegionID]
# specify aggregation operation
# will be done separately for each PC specified
aggrCommand <- paste("list(", paste(paste0(calcCols, "=", "sum", "(",
calcCols, " * pOlap)"), collapse = ", "), ")")
weightedSumByRegionDT <- pOlapDT[, eval(parse(text=aggrCommand)), by=rsRegionID]
# weightedSumByRegionDT and pOlapByRegionDT should be in the same order
regionInd <- weightedSumByRegionDT$rsRegionID
olCoord <- grToDt(regionSet)[regionInd, .(chr, start, end)]
# divide by total proportion overlap to get mean value
jCommand <- paste("list(", paste(paste0(calcCols, "=", calcCols, " / regionPOlap"), collapse = ", "), ")")
meanPerRegion <- cbind(pOlapByRegionDT, weightedSumByRegionDT)
meanPerRegion <- meanPerRegion[, eval(parse(text=jCommand))]
meanPerRegion <- cbind(olCoord, meanPerRegion)
if (returnQuantile) {
for (i in seq_along(calcCols)) {
# perhaps this could be more efficient with mapply
# ecdf example: ecdf(distributionData)(getPercentileOfThis)
meanPerRegion[, c(calcCols[i]) := ecdf(signalDT[[calcCols[i]]])(meanPerRegion[[calcCols[i]]])]
# I tried set() to improve performance but it took about the same time
}
}
# meanPerRegion <- pOlapDT[, .(regionMean = sum(score * (pOlap/sum(pOlap))), by=rsRegionID]
return(meanPerRegion)
}
#' Get regions that are most associated with target variable
#'
#' Get a GRanges with top regions from the region set based on
#' average feature contribution scores
#' for the regions or the quantile of the region's average
#' feature contribution score based on the
#' distribution of all feature contribution scores for the target variable.
#' Returns average feature contribution score or quantile as GRanges metadata.
#'
#' @template signal
#' @template signalCoord
#' @template regionSet
#' @template signalCol
#' @param cutoff Numeric. Only regions with at least this value will be
#' returned (either above this average `signal` value or above this quantile
#' if returnQuantile=TRUE).
#' @param returnQuantile Logical. If FALSE, return region averages. If TRUE,
#' for each region, return the quantile of that region's average value
#' based on the distribution of individual feature values in `signal` for
#' that `signalCol`.
#' @return A GRanges object with region coordinates for regions with
#' scores/quantiles above "cutoff" for any target variable in signalCol.
#' The scores/quantiles
#' for signalCol are given as metadata in the GRanges.
# Are regions in order along the rows of the data.table?
#
#' @examples
#' data("brcaATACCoord1")
#' data("brcaATACData1")
#' data("esr1_chr1")
#' featureContributionScores <- prcomp(t(brcaATACData1))$rotation
#' topRegions <- getTopRegions(signal=featureContributionScores,
#' signalCoord=brcaATACCoord1,
#' regionSet=esr1_chr1,
#' returnQuantile = TRUE)
#' @export
getTopRegions <- function(signal,
signalCoord,
regionSet,
signalCol = c("PC1", "PC2"),
cutoff = 0.8,
returnQuantile=TRUE) {
regionLoadDT <- averagePerRegion(signal=signal,
signalCoord=signalCoord, regionSet=regionSet,
signalCol = signalCol,
returnQuantile = returnQuantile)[]
keepInd <- regionLoadDT[, signalCol, with=FALSE] >= cutoff
# keep region if it is above cutoff in any of the PCs in signalCol
keepInd <- apply(X = keepInd, MARGIN = 1, FUN = any)
highGR <- dtToGr(regionLoadDT[keepInd, ])
values(highGR) <- as.data.frame(regionLoadDT[keepInd, signalCol, with=FALSE])
return(highGR)
}
# different scoring metrics
# remniscent of LOLA:
# support (number of regions, for us regions that have at least 1 cytosine and can be scored)
# mean loading value (or could do ratio of peak to surrounding regions)
# signal to noise ratio, how big is peak compared to noise of surrounding area
# with SNR, even a small peak could have a high SNR
# @param loadingProf
# loadingProfileSNR = function() {
# # magnitude of the peak divided by standard deviation of
# # noise (signal in surrounding areas)
# }
# support is just number of regions from each input region set that overlap at all with
# the cytosines that have loading values
# findOverlaps(regionSet, dtToGr(coordinateDT))
###########################################################################
# BSaggregate -- Aggregate a BSDT across regions or region groups,
# for multiple samples at a time.
# This function is as BScombineByRegion, but can handle not only multiple
# samples in BSDT, but also simultaneously multiple region sets by passing
# a regionsGRL (GRangesList object).
# you can use jExpr to do other functions.
# Given a bisulfite data table as input, with an identifier column for
# different samples; plus a GRanges objects with regions to aggregate.
#
# @param BSDT The bisulfite data.table (output from one of the parsing
# functions for methylation calls) that you wish to aggregate. It can
# be a combined table, with individual samples identified by column passed
# to splitFactor.
# @param regionsGRL Regions across which you want to aggregate. Should be
# from a single region set. eg GRangesList(regionSet)
# @param BSCoord GRanges. Coordinates for BSDT. If NULL, then "chr" and "start"
# columns must be in BSDT.
# @param excludeGR A GRanges object with regions you want to
# exclude from the aggregation function. These regions will be eliminated
# from the input table and not counted.
# @param jExpr You can pass a custom command in the j slot to data.table
# specifying which columns to aggregate, and which functions to use. You
# can use buildJ() to build a jExpr argument easily.
# @param byRegionGroup You can aggregate by regionID or by regionGroupID;
# this reflects the regionsGRL that you pass; by default, BSAggregate will
# aggregate each region individually -- scores will then be contiguous, and
# the output is 1 row per region.
# Turn on this flag to aggregate across all region groups, making the result
# uncontiguous, and resulting in 1 row per *region group*.
# @param returnSD Whether the standard deviation of the columns of interest
# should be returned. Standard deviation for rows that overlap with region set
# and also for rows that do not overlap with region set. Not currently functional.
# @param returnOLInfo if true, include number of overlapping cpgs and number
# of overlapping regions in the result
# @param meanPerRegion Will return the mean value in each region in a data.table
# instead of averaging all regions into a single value, returnOLInfo does
# nothing if meanPerRegion=TRUE (function exits before that code, which
# expects a different data structure)
# @param returnQuantile Only used if meanPerRegion=TRUE, instead of mean
# return the quantile/percentile of the mean of each region
# in relation to the distribution of original values in BSDT
# @template rsOL
#
BSAggregate <- function(BSDT, regionsGRL, BSCoord=NULL, excludeGR=NULL,
regionsGRL.length = NULL, splitFactor=NULL,
keepCols=NULL, sumCols=NULL,
jExpr=NULL, byRegionGroup=FALSE,
keep.na=FALSE, returnSD=FALSE,
returnOLInfo=FALSE, meanPerRegion=FALSE,
returnQuantile=FALSE, rsOL=NULL) {
# Assert that regionsGRL is a GRL.
# If regionsGRL is given as a GRanges, we convert to GRL
if(is(regionsGRL, "GRanges")) {
regionsGRL <- GRangesList(regionsGRL)
} else if (!is(regionsGRL, "GRangesList") && is.null(rsOL)) {
stop("regionsGRL is not a GRanges or GRangesList object")
}
# will cause error if BSDT is only a data.frame
if (is(BSDT, "data.frame") & !is(BSDT, "data.table")) {
BSDT <- as.data.table(BSDT)
}
if (!is(BSDT, "data.table")) {
stop("BSDT must be a data.table")
}
if(! is.null(excludeGR)) {
BSDT <- BSFilter(BSDT, minReads=0, excludeGR)
}
if (returnQuantile) {
# keep all data so quantiles can be calculated
# later code will change BSDT
origBSDT <- data.table::copy(BSDT)
}
# TODO: BSdtToGRanges needs to include end coordinate!!!
if (is.null(BSCoord)) {
bsgr <- BSdtToGRanges(list(BSDT))
} else {
bsgr <- BSCoord
}
additionalColNames <- setdiff(colnames(BSDT),
c("chr","start", "end",
"hitCount","readCount", splitFactor))
# specify that "character" outcome is expected from mode by
# supplying "a" as last vapply argument (any character object would work)
colModes <- vapply(BSDT, mode, "a")
if (is.null(sumCols)) {
sumCols <- setdiff(colnames(BSDT),c("chr", "start", "end",
"strand", splitFactor, keepCols))
# Restrict to numeric columns.
sumCols <- intersect(sumCols,
names(colModes[which(colModes == "numeric")]))
}
if (is.null(rsOL)) {
# It's required to do a findoverlaps on each region individually,
# Not on a GRL, because of the way overlaps with GRLs work. So,
# we must convert the GRL to a GR, but we must keep track of which
# regions came from which group.
regionsGR <- unlist(regionsGRL)
if(is.null(regionsGRL.length)) {
# if (length(regionsGRL) > 100) {
# message("BSAggregate: Calculating sizes. You can speed this up by supplying a regionsGRL.length vector...", appendLF=FALSE)
# }
# specify that output should be numeric with vapply
# (any number would work instead of 1)
regionsGRL.length <- vapply(regionsGRL, length, 1)
# message("Done counting regionsGRL lengths.")
}
# Build a table to keep track of which regions belong to which group
# BIOC note: sapply returns a list where each item is of different length
# therefore, I'm not using vapply
# one regionID per region of "regionsGR", withinGroupID: for when regionsGRL
# has multiple GRanges, Gives the index of each region within its corresponding
# GRanges object. regionGroupID: for when regionsGRL has multiple GRanges,
# indicates which regionsGRL GRanges object the region has come from (
# index of the GRanges object in regionsGRL)
region2group <- data.table(
regionID=seq_along(regionsGR),
chr=as.vector(seqnames(regionsGR)),
start=as.vector(start(regionsGR)),
end=as.vector(end(regionsGR)),
# withinGroupID= as.vector(unlist(sapply(regionsGRL.length, seq))),
regionGroupID=rep(seq_along(regionsGRL), regionsGRL.length)) # repeats each number regionGRL.length times
setkey(region2group, regionID)
#TODO: if overlap method is single, but "end" is present, find center
# and set that to start!
# if (all(start(bsgr[[1]]) != end(bsgr[[length(unique(queryHits(hits)))1]]))) {
if (all(start(bsgr) != end(bsgr))) {
stop("BSDT start and end coordinates are not the same. Choose a different aggrMethod.")
} else {
# fo <- findOverlaps(query = bsgr[[1]], subject = regionsGR)
fo <- findOverlaps(query = bsgr, subject = regionsGR)
}
if (length(subjectHits(fo)) < 1) {
warning("Insufficient overlap between signalCoord and the region set.")
return(NULL)
}
} else {
# Build a table to keep track of which regions belong to which group
# BIOC note: sapply returns a list where each item is of different length
# therefore, I'm not using vapply
# region2group <- data.table(
# regionID=seq_len(regionsGR),
# withinGroupID= as.vector(unlist(sapply(regionsGRL.length, seq))),
# regionGroupID=rep(seq_along(regionsGRL), regionsGRL.length))
if (length(subjectHits(rsOL)) < 1) {
warning("Insufficient overlap between signalCoord and the region set.")
return(NULL)
}
# only works for when one region set is given in regionsGRL (ie does
# not work for metaregion profiles)
region2group <- data.table(
regionID=seq_len(max(subjectHits(rsOL))),
regionGroupID=rep(1, max(subjectHits(rsOL)))) # assumes only 1 region set in regionsGRL
setkey(region2group, regionID)
fo <- rsOL
}
### use info from findOverlaps to see how many individual
# cytosines (or input regions) and how many regions (in region sets)
# overlap with one another
if (returnOLInfo) {
signalCoverage <- length(unique(queryHits(fo)))
regionSetCoverage <- length(unique(subjectHits(fo)))
}
# if("end" %in% colnames(BSDT)){
# setkey(BSDT, chr, start, end)
# } else{
# setkey(BSDT, chr, start)
# }
# Gut check:
# stopifnot(all(elementMetadata(bsgr[[1]])$readCount == BSDT$readCount))
# message("Setting regionIDs...")
BSDT <- BSDT[queryHits(fo),] #restrict the table to CpGs (or input region) in any region set.
BSDT[,regionID:=subjectHits(fo)] #record which region they overlapped.
#BSDT[queryHits(fo),regionID:=subjectHits(fo)]
#if (!keep.na) {
# BSDT = BSDT[queryHits(fo),]
#}
if (is.null(jExpr)) {
cols=c(sumCols, keepCols)
funcs <- c(rep("sum", length(sumCols)), rep("unique", length(keepCols)))
jExpr <- buildJ(cols, funcs)
}
# message("jExpr: ", jExpr)
# Define aggregation column. aggregate by region or by region group?
if (byRegionGroup) {
agCol <- "regionGroupID"
} else {
agCol <- "regionID" # Default
}
# Build the by string
if (is.null(splitFactor)) {
byString <- paste0("list(regionID)")
} else {
byString <- paste0("list(", paste("regionID",
paste0(splitFactor, ""),
collapse=", ", sep=", "), ")")
}
# Now actually do the aggregate: (aggregate within each region)
# message("Combining...")
bsCombined <- BSDT[,eval(parse(text=jExpr)), by=eval(parse(text=byString))]
setkey(bsCombined, regionID)
if (meanPerRegion) {
setkey(region2group, regionID)
avPerRegion <- merge(bsCombined, region2group)
avPerRegion[, c("regionID",
# "withinGroupID",
"regionGroupID") := NULL]
if (returnQuantile) {
for (i in seq_along(sumCols)) {
# perhaps this could be more efficient with mapply
# ecdf example: ecdf(distributionData)(getPercentileOfThis)
avPerRegion[, c(sumCols[i]) := ecdf(origBSDT[[sumCols[i]]])(avPerRegion[[sumCols[i]]])]
# I tried set() to improve performance but it took about the same time
}
}
return(avPerRegion)
}
# Now aggregate across groups.
# I do this in 2 steps to avoid assigning regions to groups,
# which takes awhile. I think this preserves memory and is faster.
# Define aggregation column. aggregate by region or by region group?
# byRegionGroup=TRUE means to aggregate within each individiual GRanges from
# regions.GRL
if (byRegionGroup) {
# must set allow=TRUE here in case there are multiple IDs (splitCol)
bsCombined[region2group, regionGroupID:=regionGroupID, allow=TRUE]
if (! is.null(splitFactor) ) {
byStringGroup <- paste0("list(", paste("regionGroupID",
paste0(splitFactor,
collapse=", "),
sep=", "), ")")
} else {
byStringGroup <- "list(regionGroupID)"
}
# aggregate by regionGroupID (region set/GRanges in regionsGRL)
bsCombined=bsCombined[,eval(parse(text=jExpr)),
by=eval(parse(text=byStringGroup))]
if (returnOLInfo) {
bsCombined[, signalCoverage := signalCoverage]
bsCombined[, regionSetCoverage := regionSetCoverage]
}
return(bsCombined)
} else {
e <- region2group[bsCombined,]
setkey(e, regionID)
return(e)
}
# WARNING: There are now 2^2 ways to aggregate, sum vs mean
# at each level: across regions, then across region sets. This
# doesn't give you a choice at this point.
}
# Given a BSDT (bisulfite data.table), remove any entries that overlap
# regions given in the excludeGR argument and/or filter out sites
# that have lower than a minimum number of reads.
# @param BSDT Bisulfite data.table to filter
# @param minReads Require at least this level of coverage at a cpg.
# @param excludeGR GRanges object with regions to filter.
#
# @return The BSDT with appropriate regions removed.
BSFilter <- function(BSDT, minReads = 10, excludeGR = NULL) {
# First, filter for minimum reads.
if (minReads > 0) {
BSDT <- BSDT[coverage >= minReads, ]
}
if (NROW(BSDT) == 0) { return(data.table(NULL)) }
# Now, filter entries overlapping a region in excludeGR.
if (!is.null(excludeGR)) {
gr <- dtToGr(BSDT)
fo <- findOverlaps(gr, excludeGR)
qh <- unique(queryHits(fo))
length(qh)
nrow(BSDT)
BSDT <- BSDT[-qh, ]
}
return(BSDT)
}
#' Get indices for top scored region sets
#'
#' For each target variable, get index of original region sets
#' but ordered by rsScores
#' ranking for each target variable.
#' The original index refers to that region set's position
#' in the `GRList` param given to `aggregateSignalGRList` which is
#' also that region set's
#' row index in the COCOA output. The first number in a given column
#' of this function's output will be the
#' original index of the region set ranked first for that target variable.
#' The second row for a
#' column will be the original index of the region set that ranked second
#' for that target variable, etc. You can use this function to make it easier
#' when you want to select the top region sets for further analysis or
#' just for sorting the results. Region set scores are sorted in decreasing
#' or increasing order according to the `decreasing` parameter.
#'
#' @template rsScores
#' @templateVar isRSRankingIndex TRUE
#' @template signalCol
#' @param decreasing Logical. Whether to sort rsScores in decreasing
#' or increasing order.
#' @param newColName Character. The names of the columns of the output data.frame.
#' The order should correspond to the order of the
#' input columns given by signalCol.
#' @return A data.frame with one column for each `signalCol`.
#' Column names are given by `signalCol` or `newColName` (if used).
#' Each column has been
#' sorted by score for region sets for that target variable
#' (order given by `decreasing`
#' param).
#' Original indices for region sets that were used to create rsScores
#' are given. Region sets with a score of NA are counted as having the
#' lowest scores and indices for these region sets will be at the bottom of the
#' returned data.frame (na.last=TRUE in sorting)
#' @examples data("rsScores")
#' rsRankInd = rsRankingIndex(rsScores=rsScores,
#' signalCol=c("PC1", "PC2"))
#' # region sets sorted by score for PC1
#' rsScores[rsRankInd$PC1, ]
#' # region sets sorted by score for PC2
#' rsScores[rsRankInd$PC2, ]
#'
#' @export
#'
rsRankingIndex <- function(rsScores, signalCol,
decreasing=TRUE, newColName = signalCol) {
if (!(is(rsScores, "data.frame") || is(rsScores, "matrix"))) {
stop("rsScores should be a data.frame. Check object class.")
}
rsScores <- as.data.table(rsScores)
if (!(is(signalCol, "character") | is(signalCol, "list"))) {
stop("signalCol should be a character vector or list of character vectors.")
}
if (is(signalCol, "list")) {
if (!all(vapply(X = signalCol, FUN = class, FUN.VALUE = "a") %in% "character")) {
stop("Items of signalCol should be character vectors.")
}
if (!length(unique(vapply(signalCol, FUN = length, FUN.VALUE = 2)))) {
stop("Items of signalCol should be the same length as each other.")
}
if (!all(unlist(signalCol) %in% colnames(rsScores))) {
stop("Some column names in signalCol are not present in rsScores.")
}
# the first item in the list is taken for newColName
if (is(newColName, "list")) {
newColName <- signalCol[[1]]
}
} else { # signalCol is character
if (!all(signalCol %in% colnames(rsScores))) {
stop("Some column names in signalCol are not present in rsScores.")
}
}
dtOrder <- rep(-99L, length(decreasing))
# how to sort scores
# -1 for decreasing order of scores
dtOrder[decreasing] <- -1L
# +1 for increasing order of scores
dtOrder[!decreasing] <- 1L
dtOrder <- rep(-99L, length(decreasing))
# how to sort scores
# -1 for decreasing order of scores
dtOrder[decreasing] <- -1L
# +1 for increasing order of scores
dtOrder[!decreasing] <- 1L
# so by references changes will not be a problem
rsScores <- copy(rsScores)
rsScores[, rsIndex := seq_len(nrow(rsScores))]
if (is(signalCol, "list")) {
if (length(newColName) != length(signalCol[[1]])) {
stop("newColName is not the same length as columns given in signalCol.")
}
rsEnSortedInd <- subset(rsScores, select= signalCol[[1]])
setnames(rsEnSortedInd, newColName)
colNameMat <- do.call(rbind, signalCol)
# then scores by each PC and make a column with the original index for sorted region sets
# this object will be used to pull out region sets that were top hits for each PC
for (i in seq_along(signalCol[[1]])) {
theseOrderCols <- colNameMat[, i]
setorderv(rsScores, cols = theseOrderCols, order=dtOrder, na.last=TRUE)
rsEnSortedInd[, newColName[i] := rsScores[, rsIndex]]
}
} else if (is(signalCol, "character")) {
if (length(newColName) != length(signalCol)) {
stop("newColName is not the same length as columns given in signalCol.")
}
signalCol <- signalCol[signalCol %in% colnames(rsScores)]
rsEnSortedInd <- subset(rsScores, select= signalCol)
setnames(rsEnSortedInd, newColName)
# then scores by each PC and make a column with the original index for sorted region sets
# this object will be used to pull out region sets that were top hits for each PC
for (i in seq_along(signalCol)) {
setorderv(rsScores, cols = signalCol[i], order=dtOrder, na.last=TRUE)
rsEnSortedInd[, newColName[i] := rsScores[, rsIndex]]
}
} else {
stop("signalCol should be a character vector or list of character vectors.")
}
# reset order
# setorderv(rsScores, cols = "rsIndex", order=1L)
return(as.data.frame(rsEnSortedInd))
}
#################### Metric functions ########################################
# scores, metrics, or statistical tests
# Instead of averaging within regions first as BSAggregate does,
# this function does a simple average and standard deviation
# for all CpGs that overlap
# with regions of a region set, also does average and
# standard deviation for non overlapping CpGs. Created to
# get metrics of loading values for each PC.
#
# Faster if given total average for each column of interest
#
# @param dataDT a data.table with
# columns to get metrics of eg (PC1, PC2). All columns
# will be considered
# columns to get the metrics from so no unnecessary columns should be
# included.
# @param dataGR GRanges. Coordinates for dataDT.
# @template regionSet
# Metrics will be calculated on
# only coordinates within this region set (and optionally separately
# on those outside this region set with alsoNonOLMet parameter)
# @param signalCol the columns to calculate the metrics on. The
# metrics will be calculated on each one of these columns separately.
# @param metrics character vector with the name of a function or functions
# to calculate on selected cytosines. Function should only require one
# input which should be the values of the cytosines.
# @param alsoNonOLMet also include same metrics
# for non overlapping CpGs (still also returns metrics for overlapping CpGs)
# param columnMeans Not functional/is deprecated. The idea was to use
# the mean of the column to speed up calculations (eg when calculating
# mean of overlapping CpGs, use that info and column mean to get
# mean for non overlapping CpGs without manually calculating it)
#
signalOLMetrics <- function(dataDT,
dataGR,
regionSet,
signalCol = colnames(dataDT)[!(colnames(dataDT) %in% c("chr", "start", "end"))],
metrics=c("mean", "sd"),
alsoNonOLMet=TRUE, rsOL=NULL) {
# convert DT to GR for finding overlaps
# dataGR <- BSdtToGRanges(list(dataDT))[[1]]
if (is.null(rsOL)) {
OL <- findOverlaps(query = dataGR, subject = regionSet)
# region set info
totalRegionNumber <- length(regionSet)
meanRegionSize <- round(mean(width(regionSet)), 1)
} else {
OL <- rsOL
}
# if no overlap, exit
if (length(OL) == 0) {
return(NULL)
}
# get indices for overlapping and non overlapping CpGs
olCpG <- queryHits(OL)
# get info on degree of overlap
# number of CpGs that overlap
signalCoverage <- length(unique(olCpG))
# number of regions that overlap
regionSetCoverage <- length(unique(subjectHits(OL)))
# gets metrics for all columns except chr, start, end
jExpr <- buildJ(cols=rep(signalCol, each=length(metrics)),
funcs=rep(metrics, length(signalCol)),
newColNames = paste0(rep(signalCol,
each=length(metrics)),
"_", metrics))
# getting the metrics
olMetrics <- as.data.frame(dataDT[olCpG, eval(parse(text=jExpr))])
# calculate average of nonOLCpGs based on columnMean if given
# if (!is.null())
#
# formatting so there is one row per PC/testCol
# output is a matrix with ncol = length(metrics)
# for vapply, FUN.VALUE should have length equal to a single output of FUN
olResults <- vapply(X = metrics,
FUN = function(x) as.numeric(olMetrics[, grepl(pattern = x, colnames(olMetrics))]),
as.numeric(seq_along(signalCol)))
olResults <- as.data.table(olResults)
setnames(olResults, old = colnames(olResults), new = paste0(colnames(olResults), "_OL"))
if (alsoNonOLMet) {
nonOLCpG <- (seq_len(nrow(dataDT)))[-olCpG]
# if no OL for region set, don't calculate for non region set
# TODO make conditional on region set having any overlap
nonOLMetrics <- as.data.frame(dataDT[nonOLCpG, eval(parse(text=jExpr))])
nonOLResults <- vapply(X = metrics,
FUN = function(x) as.numeric(nonOLMetrics[, grepl(pattern = x, colnames(nonOLMetrics))]),
as.numeric(seq_along(signalCol)))
nonOLResults <- as.data.table(nonOLResults)
setnames(nonOLResults, old = colnames(nonOLResults), new = paste0(colnames(nonOLResults), "_nonOL"))
if (is.null(rsOL)) {
metricDT <- cbind(data.table(testCol=signalCol),
olResults,
nonOLResults,
data.table(signalCoverage,
regionSetCoverage,
totalRegionNumber,
meanRegionSize))
} else {
metricDT <- cbind(data.table(testCol=signalCol),
olResults,
nonOLResults,
data.table(signalCoverage,
regionSetCoverage))
}
} else {
if (is.null(rsOL)) {
metricDT <- cbind(data.table(testCol=signalCol),
olResults,
data.table(signalCoverage,
regionSetCoverage,
totalRegionNumber,
meanRegionSize))
} else {
metricDT <- cbind(data.table(testCol=signalCol),
olResults,
data.table(signalCoverage,
regionSetCoverage))
}
}
return(metricDT)
}
# Wilcoxon rank sum test for a region set
# @param dataDT a data.table with chr, start, end columns as well
# as columns to get metrics of eg (PC1, PC2). All columns
# except chr, start, and end will be considered
# columns to get the metrics from so no unnecessary columns should be
# included.
# @template regionSet
# @param signalCol the columns of interest. You will do ranksum test separately
# on each of these columns (given test only uses info in one column)
# @param ... Additional parameters of wilcox.test function. See ?wilcox.test.
# For instance specify alternative hypothesis: alternative = "greater".
# @return A vector with a p value for each column other than chr, start or end.
# @examples data("brcaLoadings1")
# data("brcaMCoord1")
# data("nrf1_chr1")
# dataDT = as.data.table(cbind(brcaMCoord1, brcaLoadings1))
# rsWilcox(dataDT = dataDT, regionSet = nrf1_chr1, conf.int=TRUE)
rsWilcox <- function(dataDT,
regionSet,
signalCol = colnames(dataDT)[!(colnames(dataDT) %in% c("chr", "start", "end"))],
conf.int=FALSE,
...) {
# region set info
totalRegionNumber <- length(regionSet)
meanRegionSize <- round(mean(width(regionSet)), 1)
# convert DT to GR for finding overlaps
dataGR <- BSdtToGRanges(list(dataDT))[[1]]
OL <- findOverlaps(query = regionSet, subject = dataGR)
# if no overlap, exit
if (length(OL) == 0) {
return(NULL)
}
# get indices for overlapping and non overlapping CpGs
olCpG <- unique(subjectHits(OL))
nonOLCpG <- (seq_len(nrow(dataDT)))[-olCpG]
# get info on degree of overlap
# number of CpGs that overlap
signalCoverage <- length(unique(olCpG))
# number of regions that overlap
regionSetCoverage <- length(unique(queryHits(OL)))
# each confidence interval has length of 2: [low, high]
if (conf.int) {
# calculate Wilcoxon rank sum test for each column
# additional parameters given with ...
# confIntervals will be [low1, high1, low2, high2, etc.]
confIntervals <- as.numeric(vapply(X = signalCol, FUN = function(x) wilcox.test(x = as.numeric(as.matrix(dataDT[olCpG, x, with=FALSE])),
y = as.numeric(as.matrix(dataDT[nonOLCpG, x, with=FALSE])),
conf.int = conf.int, ...)$conf.int, c(1, 1)))
names(confIntervals) <- paste0(rep(signalCol, each=2), c("_low", "_high"))
wRes <- data.frame(t(confIntervals),
signalCoverage,
regionSetCoverage,
totalRegionNumber,
meanRegionSize)
} else {
# calculate Wilcoxon rank sum test for each column
# additional parameters given with ...
pVals <- vapply(X = signalCol, FUN = function(x) wilcox.test(x = as.numeric(as.matrix(dataDT[olCpG, x, with=FALSE])),
y = as.numeric(as.matrix(dataDT[nonOLCpG, x, with=FALSE])), ...)$p.value, 1)
wRes <- data.frame(t(pVals),
signalCoverage,
regionSetCoverage,
totalRegionNumber,
meanRegionSize)
}
return(wRes)
}
# I optimized upstream so that a matrix would be given to this function
# if this function is rewritten and no longer requires a matrix input,
# then in order to prevent unnecessary object copying,
# rewrite upstream code that converts signalDT to matrix class
# @param signalMat Data to be aggregated (e.g. raw data: ATAC-seq,
# region based DNA methylation or loading values)
# @param signalGR GRanges object with coordinates for signalMat
# @template regionSet
# The region set to score.
# @param calcCols character object. Column names. A weighted sum will be done
# for each of these columns (columns should be numeric).
# @template rsOL
# @param pOlap see "?aggregateSignal"
# @template returnCovInfo
# @value Returns data.frame with columns 'calcCols', signalCoverage col has
# number of signalGR regions that overlapped with any regionSet regions,
# regionSetCoverage has the sum of all proportion overlaps of regions from
# signalGR with regionSet (regionSet region is denominator)
# containing weighted mean for each col.
# Returns NULL if there is no overlap between signalGR and regionSet
regionOLWeightedMean <- function(signalMat, signalGR,
regionSet, calcCols, rsOL=NULL,
pOlap=NULL, returnCovInfo=TRUE) {
if (!is(signalMat, "matrix")) {
signalMat <- as.matrix(signalMat)
}
if (is.null(rsOL)) {
hits <- findOverlaps(query = signalGR, subject = regionSet)
} else {
hits <- rsOL
}
# if no overlap, return NULL
if (length(hits) == 0) {
return(NULL)
}
if (is.null(pOlap)) {
olap <- pintersect(signalGR[queryHits(hits)],
regionSet[subjectHits(hits)])
pOlap <- width(olap) / width(regionSet[subjectHits(hits)])
}
# some rows may be duplicated if a signalMat region overlapped multiple
# regions from signalGR but that is ok
# signalMat <- signalMat[queryHits(hits), ] # done in next step to prevent extra copying
# weight the signalMat values by the proportion overlap (weighted average)
weightedSum <- t(pOlap) %*% signalMat[queryHits(hits), calcCols]
# weighted average
denom <- sum(pOlap)
weightedAve <- as.data.frame(weightedSum / denom)
colnames(weightedAve) <- calcCols
# add columns for coverage info
if (returnCovInfo) {
weightedAve$signalCoverage = length(unique(queryHits(hits)))
weightedAve$regionSetCoverage = length(unique(subjectHits(hits)))
weightedAve$sumProportionOverlap = denom
}
return(weightedAve)
}
# @param signalDT Data to be aggregated (e.g. raw data: ATAC-seq,
# region based DNA methylation or loading values)
# @param signalGR GRanges object with coordinates for signalDT
# @template regionSet
# The region set to score.
# @param calcCols character object. Column names. A mean will be calculated for
# each of these columns (columns should be numeric).
# @param metric character. "mean" or "median"
# @template rsOL
# @template returnCovInfo
# @value Returns data.frame with columns 'calcCols', signalCoverage col has
# number of signalGR regions that overlapped with any regionSet regions,
# regionSetCoverage has the number of regions from
# signalGR that overlapped with regionSet
# Returns NULL if there is no overlap between signalGR and regionSet
regionOLMean <- function(signalDT, signalGR, regionSet,
calcCols, metric="mean", rsOL=NULL, returnCovInfo=TRUE) {
if (is.null(rsOL)) {
hits <- findOverlaps(query = signalGR, subject = regionSet)
} else {
hits <- rsOL
}
# if no overlap, return NULL
if (length(hits) == 0) {
return(NULL)
}
# some rows may be duplicated if a signalDT region overlapped multiple
# regions from signalGR but that is ok
if (metric == "mean") {
# mean of the overlapping signalDT values
signalAve <- as.data.frame(t(colMeans(signalDT[queryHits(hits),..calcCols])))
} else if (metric == "median") {
# median of the overlapping signalDT values
signalAve <- as.data.frame(t(apply(X = signalDT[queryHits(hits),..calcCols], 2, median)))
} else {
stop("Error in regionOLMean function. Invalid metric specified.")
}
if (returnCovInfo) {
# add columns for coverage info
signalAve$signalCoverage <- length(unique(queryHits(hits)))
signalAve$regionSetCoverage <- length(unique(subjectHits(hits)))
}
return(signalAve)
}
databio/PCRSA documentation built on Sept. 2, 2023, 11:24 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions weightedAvePerRegion averagePerRegion BSBinAggregate .makeSymmetric makeSymmetric getMetaRegionProfile createCorFeatureMat aggregateSignalGRList .formatResults aggregateSignal
Documented in aggregateSignal aggregateSignalGRList getMetaRegionProfile
# PACKAGE DOCUMENTATION
#' Coordinate Covariation Analysis (COCOA)
#'
#'
#' COCOA is a method for understanding epigenetic variation among samples.
#' COCOA can be used with epigenetic data that includes
#' genomic coordinates and an epigenetic signal,
#' such as DNA methylation and chromatin accessibility
#' data.
#' To describe the method on a high level, COCOA quantifies
#' inter-sample variation with either a supervised or unsupervised
#' technique then uses a database of "region sets" to
#' annotate the variation among samples. A region set is a set of
#' genomic regions that share a biological annotation,
#' for instance transcription factor (TF) binding regions,
#' histone modification regions, or open chromatin regions.
#' COCOA can identify region sets that are associated with
#' epigenetic variation between samples and
#' increase understanding of variation in your data.
#
# In contrast to some other common techniques, COCOA offers both
# supervised (known groups/phenotype) and unsupervised (no known groups/
# phenotype) analyses. Also, COCOA focuses on continuous variation
# between samples instead of having cutoffs. Because of this, COCOA can
# be used as a complementary method alongside "differential" methods
# that find discrete differences between groups of samples and
# it can also be used in situations where there are no groups.
# COCOA can identify biologically meaningful
# sources of variation between samples and
#'
#' @docType package
#' @name COCOA
#' @author John Lawson
#' @author Nathan Sheffield
#'
#' @references \url{http://github.com/databio}
#' @importFrom ggplot2 ggplot aes facet_wrap geom_boxplot geom_jitter geom_line
#' theme_classic xlab ylab geom_hline ylim scale_color_discrete
#' scale_x_discrete scale_fill_brewer scale_color_manual
#' scale_color_brewer theme element_line element_text geom_point
#' @importFrom ComplexHeatmap Heatmap draw
#' @import BiocGenerics S4Vectors IRanges GenomicRanges simpleCache
#' @importFrom data.table ":=" setDT data.table setkey fread setnames
#' setcolorder rbindlist setattr setorder copy is.data.table
#' setorderv as.data.table
#' @importFrom Biobase sampleNames
#' @importFrom stats lm coefficients poly wilcox.test ecdf pgamma p.adjust
#' @importFrom methods is
#' @importFrom MIRA binRegion
#' @importFrom tidyr gather
#' @importFrom grid grid.newpage grid.grabExpr grid.draw popViewport
#' pushViewport unit viewport
#' @importFrom grDevices dev.off
#' @importFrom methods hasArg
#' @importFrom fitdistrplus fitdist
#' @importFrom simpleCache simpleCache
NULL
# @importFrom ppcor pcor.test
# now package lists GenomicRanges in "Depends" instead of "Imports" in
# DESCRIPTION, still import package with @import though
# @importFrom GenomicRanges GRanges GRangesList elementMetadata strand
# seqnames granges
# Because of some issues,
# (see here: http://stackoverflow.com/questions/9439256/)
# I have to register stuff used in data.table as non-standard evaluation,
# in order to pass some R CMD check NOTES.
if (getRversion() >= "2.15.1") {
utils::globalVariables(c(
".", "..calcCols", "bin", "binID", "chr", "id", "colsToAnnotate",
"coordinateDT",
"coverage", "Group", "pOlap", "regionGroupID", "regionID", "theme",
"meanRegionSize", "regionSetCoverage", "rowIndex", "rsIndex",
"rsRegionID", "totalRegionNumber", "signalCoverage", ".SD",
"sumProportionOverlap"))
}
#########################################################################
#' Score a region set using feature contribution scores
#'
#' First, this function identifies which epigenetic features
#' overlap the region set.
#' Then the region set is scored using the feature contribution scores
#' (`signal` input)
#' according to the `scoringMetric` parameter.
#'
#' @template signal
#' @template signalCoord
#' @template signalCoordType
#' @templateVar refGenomeWarning TRUE
#' @template regionSet
#' @template signalCol
#' @template scoringMetric
#' @template verbose
# Useful when using
# aggregateSignal with 'apply' to do many region sets at a time.
# @param wilcox.conf.int logical. Only applies when using "rankSum" scoring
# method. returns a 95% confidence interval from the Wilcoxon rank sum test
# instead of p value.
#' @template absVal
#' @template rsOL
#' @param pOlap Numeric vector. Only used if rsOL is given and scoringMetric
#' is "proportionWeightedMean". This vector should contain the proportion of
#' each regionSet region that is overlapped by a signalCoord region. The
#' order of pOlap should be the same as the overlaps in rsOL.
#' @template returnCovInfo
#' @template checkInput
#' @return a data.frame with various columns:
#' One column for each item of signalCol, representing scores for the region set for each signalCol.
#' Columns: signalCoverage (cytosine_coverage) shows the number of epigenetic features that overlapped with the region set.
#' regionSetCoverage: indicates the number of regions from regionSet that overlapped with any of the epigenetic features.
#' totalRegionNumber: represents the number of regions in regionSet.
#' meanRegionSize: gives the average size (in base pairs) of regions in regionSet.
#'
#' For "multiBase" data using the "proportionWeightedMean" scoring metric,
#' the output will also include a "sumProportionOverlap" column.
#' This column represents the sum of proportion overlaps between
#' each signalCoord region and overlapping regionSet region,
#' providing an additional measure of regionSet coverage.
#'
#' @examples
#' data("brcaATACCoord1")
#' data("brcaATACData1")
#' data("esr1_chr1")
#' featureContributionScores <- prcomp(t(brcaATACData1))$rotation
#' rsScores <- aggregateSignal(signal=featureContributionScores,
#' signalCoord=brcaATACCoord1,
#' regionSet=esr1_chr1,
#' signalCol=c("PC1", "PC2"),
#' scoringMetric="default")
#' @export
aggregateSignal <- function(signal,
signalCoord,
regionSet,
signalCol = c("PC1", "PC2"),
signalCoordType = "default",
scoringMetric = "default",
verbose = FALSE,
absVal=TRUE,
rsOL=NULL, pOlap=NULL,
returnCovInfo=TRUE,
.checkInput=FALSE) {
################### checking inputs #################################
# if it was already checked outside this function, don't need to re-check
if (.checkInput) {
########## check that inputs are the correct class
checkConvertInputClasses(signal=signal,
signalCoord=signalCoord,
regionSet=regionSet,
signalCol=signalCol,
rsOL=rsOL)
########## check that dimensions of inputs are consistent
# length of signal coord = nrow of signal
if (length(signalCoord) != nrow(signal)) {
stop(cleanws("The number of coordinates in
signalCoord (length(signalCoord)) does not equal the number of
rows in `signal`"))
}
######### check that appropriate columns are present
# signalCol are column names of signal
if (!all(signalCol %in% colnames(signal))) {
missingCols = signalCol[!(signalCol %in% colnames(signal))]
stop(cleanws(paste0("Some signalCol are not
columns of signal: ", missingCols)))
}
######## check that scoringMetric is appropriate
if (!(scoringMetric %in% getScoringMethods("both"))) {
stop(cleanws("scoringMetric was not recognized.
Check spelling and available options."))
}
###### check that signalCoordType is appropriate
if (!(signalCoordType %in% c("default", "singleBase", "multiBase"))) {
stop(cleanws("signalCoordType not recognized.
Check spelling/capitalization."))
}
#######
# what happens if there are NAs or Inf in `signal`?
# any NAs that overlap the regionSet will cause the score to be NA
if (is(signal, "data.table")) {
naRows = apply(X = signal[, signalCol, with=FALSE, drop=FALSE],
MARGIN = 1, FUN = function(x) any(is.na(x)))
} else {
naRows = apply(X = signal[, signalCol, drop=FALSE],
MARGIN = 1, FUN = function(x) any(is.na(x)))
}
if (any(naRows)) {
signal <- signal[!naRows, ]
signalCoord <- signalCoord[!naRows]
warning("Removing rows with NA from `signal`")
}
#################################################################
# detect signalCoordType
if (signalCoordType == "default") {
# when signalCoord is a GRanges object
if (any(start(signalCoord) != end(signalCoord))) {
signalCoordType <- "multiBase"
} else {
signalCoordType <- "singleBase"
}
}
# if "default" scoring method is given, choose based on signalCoordType
if (scoringMetric == "default") {
if (signalCoordType == "singleBase") {
scoringMetric <- "regionMean"
} else if (signalCoordType == "multiBase") {
scoringMetric <- "proportionWeightedMean"
} else {
stop(cleanws("signalCoordType not recognized.
Check spelling/capitalization."))
}
}
# make sure that scoringMetric is consistent with signalCoordType
if (signalCoordType == "singleBase") {
if (!(scoringMetric %in% getScoringMethods("singleBase"))) {
stop("The scoringMetric you selected is not available for
this data's signalCoordType")
}
} else if (signalCoordType == "multiBase") {
if (!(scoringMetric %in% getScoringMethods("multiBase"))) {
stop("The scoringMetric you selected is not available for
this data's signalCoordType")
}
}
}
# extreme positive or negative values both give important information
# take absolute value or not
if (absVal) {
signal <- abs(signal) # required for later code
}
# XX copies unnecessarily:reformat into data.table with chromosome location and weight
# make sure `signal` is the correct type for further steps
# (proportionWeightedMean method requires a matrix)
if (!is(signal, "data.table") && (scoringMetric != "proportionWeightedMean")) {
signal <- as.data.table(signal)
} else if (!is(signal, "matrix") && (scoringMetric == "proportionWeightedMean")) {
signal <- as.matrix(signal)
}
# restricting to signalCol so unnecessary computations
# are not done
if (is(signal, "data.table")) {
loadingDT <- signal[, signalCol, with=FALSE]
# # naming does not work if only using one PC so add this line for that case
# setnames(loadiangDT, signalCol)
} else {
loadingDT <- signal[, signalCol, drop=FALSE]
}
###########################################################################
# scoring
# would rounding speed up aggregation?, potentially make a sparse matrix
# if a lot of entries became 0
# works for both singleBase and multiBase
if (scoringMetric == "simpleMean") {
loadAgMain <- as.data.table(regionOLMean(signalDT = loadingDT,
signalGR = signalCoord,
regionSet = regionSet,
calcCols= signalCol,
metric = "mean",
rsOL = rsOL,
returnCovInfo=returnCovInfo))
results <- .formatResults(loadAgMain, scoringMetric = scoringMetric,
regionSet=regionSet, signalCol = signalCol,
returnCovInfo=returnCovInfo)
} else if (scoringMetric == "simpleMedian") {
# scoring singleBase and multiBase both with this function for
# simpleMedian
loadAgMain <- as.data.table(regionOLMean(signalDT = loadingDT,
signalGR = signalCoord,
regionSet = regionSet,
calcCols= signalCol,
metric = "median",
rsOL=rsOL,
returnCovInfo = returnCovInfo))
results <- .formatResults(loadAgMain, scoringMetric = scoringMetric,
regionSet=regionSet, signalCol = signalCol,
returnCovInfo=returnCovInfo)
} else if (signalCoordType == "singleBase") {
# do the actual aggregation
if (scoringMetric == "regionMean") {
# specify aggregation operation
# will be done separately for each PC specified
aggrCommand <- buildJ(signalCol,
rep("mean", length(signalCol)))
# previously used BSAggregate from RGenomeUtils but now using local,
# modified copy
loadAgMain <- BSAggregate(BSDT = loadingDT,
regionsGRL = regionSet,
BSCoord = signalCoord,
jExpr = aggrCommand,
byRegionGroup = TRUE,
splitFactor = NULL,
returnOLInfo = returnCovInfo,
rsOL=rsOL)
results <- .formatResults(loadAgMain,
scoringMetric = scoringMetric,
regionSet=regionSet, signalCol = signalCol,
returnCovInfo=returnCovInfo)
} else if (scoringMetric == "regionMedian") {
aggrCommand <- buildJ(signalCol,
rep("median", length(signalCol)))
loadAgMain <- BSAggregate(BSDT = loadingDT,
regionsGRL = regionSet,
BSCoord = signalCoord,
jExpr = aggrCommand,
byRegionGroup = TRUE,
splitFactor = NULL,
returnOLInfo = returnCovInfo)
results <- .formatResults(loadAgMain,
scoringMetric = scoringMetric,
regionSet=regionSet, signalCol = signalCol,
returnCovInfo=returnCovInfo)
}
# } else if (scoringMetric == "meanDiff") {
# # if (is.null(pcLoadAv)) {
# # # calculate (should already be absolute)
# # pcLoadAv <- apply(X = loadingDT[, signalCol, with=FALSE],
# # MARGIN = 2, FUN = mean)
# # }
# loadMetrics <- signalOLMetrics(dataDT=loadingDT, regionSet=regionSet,
# signalCol = signalCol,
# metrics=c("mean", "sd"),
# alsoNonOLMet=TRUE)
# if (is.null(loadMetrics)) {
# results <- as.data.table(t(rep(NA, length(signalCol))))
# setnames(results, signalCol)
# results[, signalCoverage := 0]
# results[, regionSetCoverage := 0]
# results[, totalRegionNumber := numOfRegions]
# results[, meanRegionSize := round(mean(width(regionSet)), 1)]
# } else {
# # calculate mean difference
# # pooled standard deviation
# sdPool <- sqrt((loadMetrics$sd_OL^2 + loadMetrics$sd_nonOL^2) / 2)
#
# # mean difference
# # error if signalCoverage > (1/2) * totalCpGs
# meanDiff <- (loadMetrics$mean_OL - loadMetrics$mean_nonOL) /
# (sdPool * sqrt((1 / loadMetrics$signalCoverage) - (1 / (totalCpGs - loadMetrics$signalCoverage))))
# results <- as.data.table(t(meanDiff))
# colnames(results) <- loadMetrics$testCol
#
# # add information about degree of overlap
# results <- cbind(results, loadMetrics[1, .SD, .SDcols = c("signalCoverage", "regionSetCoverage", "totalRegionNumber", "meanRegionSize")])
# }
#
# } else if (scoringMetric == "rankSum") {
#
# # if (wilcox.conf.int) {
# # # returns confidence interval
# # wRes <- rsWilcox(dataDT = loadingDT, regionSet=regionSet,
# # signalCol = signalCol,
# # conf.int = wilcox.conf.int)
# #
# # if (is.null(wRes)) {
# # results <- as.data.table(t(rep(NA, length(signalCol) * 2)))
# # setnames(results, paste0(rep(signalCol, each=2),
# # c("_low", "_high")))
# # results[, signalCoverage := 0]
# # results[, regionSetCoverage := 0]
# # results[, totalRegionNumber := numOfRegions]
# # results[, meanRegionSize := round(mean(width(regionSet)), 1)]
# # } else {
# # results <- as.data.table(wRes)
# # }
# #
# # } else {
# # returns p value
# # one sided test since I took the absolute value of the loadings
# wRes <- rsWilcox(dataDT = loadingDT, regionSet=regionSet,
# signalCol = signalCol, alternative="greater")
#
# if (is.null(wRes)) {
# results <- as.data.table(t(rep(NA, length(signalCol))))
# setnames(results, signalCol)
# results[, signalCoverage := 0]
# results[, regionSetCoverage := 0]
# results[, totalRegionNumber := numOfRegions]
# results[, meanRegionSize := round(mean(width(regionSet)), 1)]
# } else {
# results <- as.data.table(wRes)
# }
# #}
} else {
# signalCoordType == "multiBase"
# for ATAC-seq
if (scoringMetric == "proportionWeightedMean") {
loadAgMain <- regionOLWeightedMean(signalMat = loadingDT,
signalGR = signalCoord,
regionSet = regionSet,
calcCols= signalCol,
rsOL = rsOL,
pOlap = pOlap,
returnCovInfo=returnCovInfo)
results <- .formatResults(as.data.table(loadAgMain),
scoringMetric = scoringMetric,
regionSet=regionSet, signalCol = signalCol,
returnCovInfo=returnCovInfo)
}
}
# signalOLMetrics() # make sure it works with no overlap
if (verbose) {
message(":", appendLF=FALSE)
}
return(as.data.frame(results))
}
# format results of scoring functions in aggregateSignal()
.formatResults <- function(loadAgMain, scoringMetric,
regionSet, signalCol, returnCovInfo=TRUE) {
numOfRegions <- length(regionSet)
# if no cytosines from loadings were included in regionSet, result is NA
if (is.null(loadAgMain)) {
results <- as.data.table(t(rep(NA, length(signalCol))))
setnames(results, signalCol)
if (returnCovInfo) {
results[, signalCoverage := 0]
results[, regionSetCoverage := 0]
results[, totalRegionNumber := numOfRegions]
results[, meanRegionSize := round(mean(width(regionSet)), 1)]
# this column is only added by this scoring method
if (scoringMetric == "proportionWeightedMean") {
results[, sumProportionOverlap := 0]
}
}
} else {
# regionMean, regionMedian, simpleMean for region data
results <- loadAgMain[, .SD, .SDcols = signalCol]
if (returnCovInfo) {
results[, signalCoverage := loadAgMain[, .SD, .SDcols = "signalCoverage"]]
results[, regionSetCoverage := loadAgMain[, .SD, .SDcols = "regionSetCoverage"]]
results[, totalRegionNumber := numOfRegions]
results[, meanRegionSize := round(mean(width(regionSet)), 1)]
################################
# proportionWeightedMean
# this column is only added by this scoring method
if (scoringMetric == "proportionWeightedMean") {
results[, sumProportionOverlap := loadAgMain[, .SD, .SDcols = "sumProportionOverlap"]]
}
}
##############################
# simple mean for base pair data
}
return(results)
}
#' Score many region sets
#'
#' This function will give each region set a score for each target variable
#' given by `signalCol` based on
#' the `scoringMetric` parameter. Based on these scores, you can determine
#' which region sets out of a region set database (given by `GRList`)
#' are most associated with the target variables. See the vignette "Introduction
#' to Coordinate Covariation Analysis" for help interpreting your
#' results.
#'
#' @template signal
#' @template signalCoord
#' @template signalCoordType
#' @template GRList
#' @template signalCol
#' @template scoringMetric
#' @template verbose
# @param wilcox.conf.int logical. Only applies when using "rankSum" scoring
# method. returns a 95% confidence interval from the Wilcoxon rank sum test
# instead of p value.
#' @template absVal
#' @template olList
#' @template pOlapList
#' @template returnCovInfo
#' @return Data.frame of results, one row for each region set.
#' It has the following columns:
#' one column for each item of signalCol with names given by signalCol.
#' These columns have scores for the region set for each signalCol.
#' Other columns: signalCoverage (formerly cytosine_coverage) which
#' has number of epigenetic features that overlapped at all with regionSet,
#' regionSetCoverage which has number of regions from regionSet
#' that overlapped any of the epigenetic features,
#' totalRegionNumber that has number of regions in regionSet,
#' meanRegionSize that has average size in base pairs of regions in regionSet,
#' the average is based on all regions in regionSet and not just ones that overlap.
#'
#' For "multiBase" data, if the "proportionWeightedMean" scoring metric
#' is used, then the output will also have a "sumProportionOverlap" column.
#' During this scoring method, the proportion overlap between each signalCoord
#' region and overlapping regionSet region is calculated. This column is
#' the sum of all those proportion overlaps and is another way to quantify
#' coverage of regionSet in addition to regionSetCoverage.
#'
#'
#' @examples
#' data("brcaATACCoord1")
#' data("brcaATACData1")
#' data("esr1_chr1")
#' data("nrf1_chr1")
#' featureContributionScores <- prcomp(t(brcaATACData1))$rotation
#' GRList <- GRangesList(esr1_chr1, nrf1_chr1)
#' rsScores <- aggregateSignalGRList(signal=featureContributionScores,
#' signalCoord=brcaATACCoord1,
#' GRList= GRList,
#' signalCol=c("PC1", "PC2"),
#' scoringMetric="default")
#'
#' @export
aggregateSignalGRList <- function(signal,
signalCoord,
GRList,
signalCol = c("PC1", "PC2"),
signalCoordType = "default",
scoringMetric = "default",
verbose = TRUE,
absVal=TRUE, olList=NULL,
pOlapList=NULL, returnCovInfo=TRUE) {
################### checking inputs #################################
########## check that inputs are the correct class
checkConvertInputClasses(signal=signal,
signalCoord=signalCoord,
regionSet=NULL,
signalCol = signalCol,
GRList=GRList,
olList=olList)
########## check that dimensions of inputs are consistent
# length of signal coord = nrow of signal
if (length(signalCoord) != nrow(signal)) {
stop(cleanws("The number of coordinates in
signalCoord (length(signalCoord)) does not equal the number of
rows in `signal`"))
}
######### check that appropriate columns are present
# signalCol are column names of signal
if (!all(signalCol %in% colnames(signal))) {
missingCols = signalCol[!(signalCol %in% colnames(signal))]
stop(cleanws(paste0("Some signalCol are not
columns of signal: ", missingCols)))
}
######## check that scoringMetric is appropriate
if (!(scoringMetric %in% getScoringMethods("both"))) {
stop(cleanws("scoringMetric was not recognized.
Check spelling and available options."))
}
###### check that signalCoordType is appropriate
if (!(signalCoordType %in% c("default", "singleBase", "multiBase"))) {
stop(cleanws("signalCoordType not recognized.
Check spelling/capitalization."))
}
#######
# what happens if there are NAs or Inf in `signal`?
# any NAs that overlap the regionSet will cause the score to be NA
if (is(signal, "data.table")) {
naRows = apply(X = signal[, signalCol, with=FALSE, drop=FALSE],
MARGIN = 1, FUN = function(x) any(is.na(x)))
} else {
naRows = apply(X = signal[, signalCol, drop=FALSE],
MARGIN = 1, FUN = function(x) any(is.na(x)))
}
if (any(naRows)) {
signal <- signal[!naRows, ]
signalCoord <- signalCoord[!naRows]
warning("Removing rows with NA from `signal`")
}
#################################################################
# detect signalCoordType
if (signalCoordType == "default") {
# when signalCoord is a GRanges object
if (any(start(signalCoord) != end(signalCoord))) {
signalCoordType <- "multiBase"
} else {
signalCoordType <- "singleBase"
}
}
# if "default" scoring method is given, choose based on signalCoordType
if (scoringMetric == "default") {
if (signalCoordType == "singleBase") {
scoringMetric <- "regionMean"
} else if (signalCoordType == "multiBase") {
scoringMetric <- "proportionWeightedMean"
} else {
stop(cleanws("signalCoordType not recognized.
Check spelling/capitalization."))
}
}
# convert object class outside aggregateSignal to extra prevent copying
# (one scoring method needs `signal` as a matrix though)
if (!is(signal, "data.table") && (scoringMetric != "proportionWeightedMean")) {
signal <- as.data.table(signal)
} else if (!is(signal, "matrix") && (scoringMetric == "proportionWeightedMean")) {
signal <- as.matrix(signal)
}
# take absolute value outside aggregateSignal to prevent extra copying
if (absVal) {
if (is(signal, "data.table")) {
signal[, signalCol] <- abs(signal[, signalCol, with=FALSE])
} else {
signal[, signalCol] <- abs(signal[, signalCol])
}
absVal <- FALSE
}
if (is.null(olList)) {
# apply over the list of region sets
resultsList <- lapplyAlias(GRList,
function(x) aggregateSignal(
signal = signal,
signalCoord = signalCoord,
signalCoordType = signalCoordType,
regionSet = x,
signalCol = signalCol,
scoringMetric = scoringMetric,
verbose = verbose,
absVal = absVal, pOlap=NULL,
returnCovInfo=returnCovInfo))
# resultsList <- mapply(FUN = function(x) aggregateSignal(
# signal = signal,
# signalCoord = signalCoord,
# signalCoordType = signalCoordType,
# regionSet = x,
# signalCol = signalCol,
# scoringMetric = scoringMetric,
# verbose = verbose,
# absVal = absVal, pOlap=y,
# returnCovInfo=returnCovInfo), x=GRList, y=pOlapList)
#wilcox.conf.int = wilcox.conf.int,
} else {
# # apply over the list of region sets
# resultsList <- lapplyAlias(olList,
# function(x) aggregateSignal(
# signal = signal,
# signalCoord = signalCoord,
# signalCoordType = signalCoordType,
# regionSet = NULL,
# signalCol = signalCol,
# scoringMetric = scoringMetric,
# verbose = verbose,
# absVal = absVal,
# rsOL=x, pOlap=pOlapList, returnCovInfo=returnCovInfo))
# #wilcox.conf.int = wilcox.conf.int,
if (is.null(pOlapList)) {
# mapply does not work if pOlapList is null and other argument is not null
resultsList <- lapplyAlias(olList, FUN = function(x) aggregateSignal(
signal = signal,
signalCoord = signalCoord,
signalCoordType = signalCoordType,
regionSet = NULL,
signalCol = signalCol,
scoringMetric = scoringMetric,
verbose = verbose,
absVal = absVal,
rsOL=x,
returnCovInfo=returnCovInfo))
} else {
resultsList <- mapply(FUN = function(x, y) aggregateSignal(
signal = signal,
signalCoord = signalCoord,
signalCoordType = signalCoordType,
regionSet = NULL,
signalCol = signalCol,
scoringMetric = scoringMetric,
verbose = verbose,
absVal = absVal,
rsOL=x, pOlap=y,
returnCovInfo=returnCovInfo),
x=olList, y=pOlapList)
}
}
resultsDT <- do.call(rbind, resultsList)
# # add names if they are present
# if (!is.null(names(GRList))) {
# # resultsDT[, rsName := names(GRList)]
# row.names(resultsDT) <- names(GRList)
# }
# Convert the first column to numeric for proper sorting
resultsDT[, 1] <- as.numeric(resultsDT[, 1])
# rsName <- row.names(resultsDT)
resultsDF <- as.data.frame(resultsDT) #, row.names = rsName)
# Sort the results based on the first column
resultsDF <- resultsDF[order(-resultsDF[, 1]), ]
return(resultsDF)
}
# @param dataMat columns of dataMat should be samples/patients, rows should be genomic signal
# (each row corresponds to one genomic coordinate/range)
# @param featureMat Rows should be samples, columns should be "features"
# (whatever you want to get correlation with: eg PC scores),
# all columns in featureMat will be used (subset when passing to function
# in order to not use all columns)
# @param centerDataMat logical object. Should rows in dataMat be centered based on
# their means? (subtracting row mean from each row)
# @param centerFeatureMat logical. Should columns in featureMat be centered based
# on their means? (subtract column mean from each column)
# @param testType character object. Can be "cor" (Pearson correlation),
# "spearmanCor (Spearman correlation)
# "pcor" (partial correlation), "cov" (covariance (Pearson)),
# @param covariate
#
# If a row in dataMat has 0 stand. deviation, correlation will be set to 0
# instead of NA as would be done by cor()
#
# @return returns a matrix where rows are the genomic signal (eg a CpG or region) and
# columns are the columns of featureMat
# @examples dataMat = matrix(rnorm(50), 5, 10)
# featureMat = matrix(rnorm(20), 10, 2)
createCorFeatureMat <- function(dataMat, featureMat,
centerDataMat=TRUE, centerFeatureMat = TRUE,
testType="cor", covariate=NULL) {
featureMat <- as.matrix(featureMat)
featureNames <- colnames(featureMat)
nFeatures <- ncol(featureMat)
nDataDims <- nrow(dataMat)
if (centerDataMat) {
cpgMeans <- rowMeans(dataMat, na.rm = TRUE)
# centering before calculating correlation
dataMat <- apply(X = dataMat, MARGIN = 2, function(x) x - cpgMeans)
}
if (centerFeatureMat) {
featureMeans <- colMeans(featureMat, na.rm = TRUE)
# centering before calculating correlation(also, t() converts to matrix)
featureMat <- t(apply(X = t(featureMat), MARGIN = 2, function(x) x - featureMeans))
if (dim(featureMat)[1] == 1) {
featureMat <- t(featureMat)
}
}
# avoid this copy and/or delay transpose until after calculating correlation?
dataMat <- as.data.frame(t(dataMat))
if (testType == "cor") {
# create feature correlation matrix with PCs (rows: features/CpGs, columns:PCs)
# how much do features correlate with each PC?
# put epigenetic data first in cor()
featurePCCor <- cor(dataMat, featureMat, use="pairwise.complete.obs", method="pearson")
} else if (testType == "spearmanCor") {
# xtfrm(x) ranking
featurePCCor <- cor(dataMat, featureMat, use="pairwise.complete.obs", method="spearman")
# } else if (testType == "pcor") {
# # partial correlation (account for covariates), ppcor package
#
# featurePCCor <- apply(X = featureMat, MARGIN = 2, function(y) apply(X = dataMat, 2,
# FUN = function(x) pcor.test(x = x, y=y,
# z=covariate,
# method="pearson")$estimate))
#
} else if (testType == "cov") {
featurePCCor <- cov(dataMat, featureMat, use="pairwise.complete.obs")
} else {
stop("invalid testType")
}
# if standard deviation of the data was zero, NA will be produced
# set to 0 because no standard deviation means no correlation with attribute of interest
featurePCCor[is.na(featurePCCor)] <- 0
colnames(featurePCCor) <- featureNames
return(featurePCCor)
# corLoadRatio <- signal[, signalCol] / featurePCCor
# hist(corLoadRatio[, "PC10"])
}
#' Create a "meta-region" profile
#'
#' This profile can show enrichment
#' of genomic signals with high feature contribution scores
#' in the region set but not in the
#' surrounding genome, suggesting that variation is linked specifically
#' to that region set.
#'
#' All regions in a given region set
#' are combined into a single aggregate profile. Regions in `regionSet`
#' should be expanded on each side to include a wider area of the genome around
#' the regions of interest (see example and vignettes).
#' To make the profile, first we optionally take
#' the absolute value of `signal` (`absVal` parameter).
#' Then each expanded regionSet region is
#' split into `binNum` bins. The corresponding
#' bins from each region (e.g. all bin1's, all bin2's, etc.) are grouped.
#' All overlapping values from `signal` are aggregated in each bin group
#' according to the `aggrMethod` parameter to get a meta-region profile.
#' Since DNA strand information is not considered,
#' the profile is averaged symmetrically around the center.
#' A peak in the middle of this profile suggests that variability is specific
#' to the region set of interest and is not a product of the surrounding genome.
#' A region set can still be significant even if it does not have a peak.
#' For example, some histone modification region sets may be in large genomic blocks
#' and not show a peak, despite having variation across samples.
#'
#' @template signal
#' @template signalCoord
#' @template regionSet
#' @template signalCol
#' @template signalCoordType
#' @param binNum Number of bins to split each region into when
#' making the aggregate profile. More bins will
#' give a higher resolution but perhaps more noisy profile.
#' @template verbose
#' @templateVar usesAggrMethod TRUE
#' @template scoringMetric
#' @template absVal
#' @return A data.frame with the binned meta-region profile,
#' one row per bin. columns: binID and one column for each target variable
#' in signalCol. The function will return NULL if there
#' is no overlap between signalCoord and any of the bin groups that come
#' from regionSet (e.g. none of the bin1's overlapped signalCoord,
#' NULL returned).
#'
#' @examples
#' data("brcaATACCoord1")
#' data("brcaATACData1")
#' data("esr1_chr1")
#' featureContributionScores <- prcomp(t(brcaATACData1))$rotation
#' esr1_chr1_expanded <- resize(esr1_chr1, 12000, fix="center")
#' mrProfile <- getMetaRegionProfile(signal=featureContributionScores,
#' signalCoord=brcaATACCoord1,
#' regionSet=esr1_chr1_expanded,
#' signalCol=c("PC1", "PC2"),
#' binNum=21)
#' @export
getMetaRegionProfile <- function(signal, signalCoord, regionSet,
signalCol = c("PC1", "PC2"),
signalCoordType = "default",
binNum = 21,
verbose=TRUE,
aggrMethod = "default", absVal=TRUE) {
################### checking inputs #################################
########## check that inputs are the correct class, convert
checkConvertInputClasses(signal=signal,
signalCoord=signalCoord,
regionSet=regionSet,
signalCol = signalCol)
########## check that dimensions of inputs are consistent
# length of signal coord = nrow of signal
if (length(signalCoord) != nrow(signal)) {
stop(cleanws("The number of coordinates in
signalCoord (length(signalCoord)) does not equal the number of
rows in `signal`"))
}
######### check that appropriate columns are present
# signalCol are column names of `signal`
if (!all(signalCol %in% colnames(signal))) {
missingCols = signalCol[!(signalCol %in% colnames(signal))]
stop(cleanws(paste0("Some signalCol are not
columns of signal: ", missingCols)))
}
######## check that aggregation method is appropriate
if (!(aggrMethod %in% getScoringMethods("metaRegionProfile"))) {
stop(cleanws("scoringMetric was not recognized.
Check spelling and available options."))
}
###### check that signalCoordType is appropriate
if (!(signalCoordType %in% c("default", "singleBase", "multiBase"))) {
stop(cleanws("signalCoordType not recognized.
Check spelling/capitalization."))
}
#######
# what happens if there are NAs or Inf in `signal`?
#################################################################
# detect signalCoordType
if (signalCoordType == "default") {
# when signalCoord is a GRanges object
if (any(start(signalCoord) != end(signalCoord))) {
signalCoordType <- "multiBase"
} else {
signalCoordType <- "singleBase"
}
}
# if "default" aggregation method is given, choose based on signalCoordType
if (aggrMethod == "default") {
if (signalCoordType == "singleBase") {
aggrMethod <- "regionMean"
} else if (signalCoordType == "multiBase") {
aggrMethod <- "proportionWeightedMean"
} else {
stop(cleanws("signalCoordType not recognized.
Check spelling/capitalization."))
}
}
##################################################################
# take absolute value or not
if (absVal) {
loadingDT <- as.data.table(abs(signal))
} else {
loadingDT <- as.data.table(signal)
}
GRDT <- grToDt(regionSet)
loadProf <- BSBinAggregate(BSDT = loadingDT,
rangeDT = GRDT,
binCount = binNum,
BSCoord = signalCoord,
byRegionGroup = TRUE,
splitFactor = NULL,
signalCol = signalCol,
aggrMethod = aggrMethod)
# if loadProf is NULL, return NULL from function, otherwise make symmetrical
# it will be NULL when there was no overlap between data and any of the bins
if (!is.null(loadProf)) {
loadProf <- makeSymmetric(loadProf)
loadProf[, regionGroupID := seq_len(binNum)][]
setnames(loadProf, old = "regionGroupID", new="binID")
loadProf <- as.data.frame(loadProf)
} else {
warning("Insufficient overlap between regionSet and signalCoord")
}
return(loadProf)
}
makeSymmetric <- function(prof) {
symProf <- apply(prof, 2, .makeSymmetric)
symProf <- as.data.table(symProf)
return(symProf)
}
.makeSymmetric <- function(vec) {
symVec <- (vec + rev(vec)) / 2
return(symVec)
}
# Produced originally for binning Ewing RRBS data across various region sets
#
# @param BSDT A data.table. For COCOA, a data.table of loading values
# with the PCs to be annotated. One column for the signal of each
# target variable
# and also has columns with the coordinates for the epigenetic features:
# chr (chromosome) and start column (also possibly end)
# @param rangeDT A data.table with the sets of regions to be binned,
# with columns named start, end
# @param binCount Number of bins across the region
# @param BSCoord GRanges. Coordinates for BSDT. If NULL, then "chr" and "start"
# columns must be in BSDT.
# @param byRegionGroup Pass along to binCount (see ?binCount)
# @template signalCol
# @param verbose A "logical" object. Whether progress
# of the function should be shown, one
# bar indicates the region set is completed.
# useful when using BSBinAggregate with 'apply' to do many
# region sets at a time.
# @param aggrMethod see ?getMetaRegionProfile()
BSBinAggregate <- function(BSDT, rangeDT, binCount,
BSCoord=NULL,
byRegionGroup = TRUE,
splitFactor = NULL,
signalCol,
verbose = FALSE,
aggrMethod) {
if (!is(rangeDT, "data.table")) {
stop("rangeDT must be a data.table")
}
if (is.null(BSCoord)) {
BSCoord <- BSdtToGRanges(list(BSDT))
}
seqnamesColName <- "seqnames" # default column name
if (! "seqnames" %in% colnames(rangeDT)) {
if ("chr" %in% colnames(rangeDT)) {
# message("seqnames column name set to: chr")
seqnamesColName <- "chr"
} else {
# Got neither.
stop("rangeDT must have a seqnames column")
}
}
# message("Binning...")
binnedDT <- rangeDT[, MIRA::binRegion(start, end,
binCount, get(seqnamesColName))]
# output is a list of GRanges objects, does not play well with vapply
# one GRanges object for each bin, containing a segment of each original rangeDT region
binnedGR <- sapply(split(binnedDT, binnedDT$binID), dtToGr)
# message("Aggregating...")
if (aggrMethod == "proportionWeightedMean") {
binMeansList <- lapply(X = binnedGR,
FUN = function(x) regionOLWeightedMean(signalMat = BSDT,
# signalGR = dtToGr(BSDT[, .(chr, start, end)]),
signalGR = BSCoord,
regionSet = x,
calcCols = signalCol))
# any bins that had no overlap with data will be NULL
# if any bins had no data, return NULL
if (any(vapply(X = binMeansList, FUN = is.null, FUN.VALUE = TRUE))) {
return(NULL)
}
# rbindlist of all NULL's will still return an (empty) data.table
# otherwise any NULL items will just be skipped and other rows will
# be concatenated
binnedBSDT <- rbindlist(binMeansList)
regionGroupID = 1:length(binMeansList)
# regionGroupID = regionGroupID[!vapply(X = binMeansList, FUN = is.null, FUN.VALUE = TRUE)]
binnedBSDT[, regionGroupID := regionGroupID]
} else if (aggrMethod == "simpleMean") {
binMeansList <- lapply(X = binnedGR,
FUN = function(x) regionOLMean(signalDT = BSDT,
# signalGR = dtToGr(BSDT[, .(chr, start, end)]),
signalGR = BSCoord,
regionSet = x,
calcCols = signalCol))
# any bins that had no overlap with data will be NULL
# if any bins had no data, return NULL
if (any(vapply(X = binMeansList, FUN = is.null, FUN.VALUE = TRUE))) {
return(NULL)
}
binnedBSDT <- rbindlist(binMeansList)
regionGroupID = 1:length(binMeansList)
# regionGroupID = regionGroupID[!vapply(X = binMeansList, FUN = is.null, FUN.VALUE = TRUE)]
binnedBSDT[, regionGroupID := regionGroupID]
} else if (aggrMethod == "regionMean") { # aggrMethod == "regionMean"
# what is output if a region set has no overlap?
binnedBSDT <- BSAggregate(BSDT,
regionsGRL=GRangesList(binnedGR),
BSCoord = BSCoord,
jExpr=buildJ(signalCol,
rep("mean", length(signalCol))),
byRegionGroup = byRegionGroup,
splitFactor = splitFactor)
# if any bins had no data
if (nrow(binnedBSDT) < binCount) {
return(NULL)
}
} else if (aggrMethod == "regionMedian") {
# what is output if a region set has no overlap?
binnedBSDT <- BSAggregate(BSDT,
regionsGRL=GRangesList(binnedGR),
BSCoord = BSCoord,
jExpr=buildJ(signalCol,
rep("median", length(signalCol))),
byRegionGroup = byRegionGroup,
splitFactor = splitFactor)
# if any bins had no data
if (nrow(binnedBSDT) < binCount) {
return(NULL)
}
}
# RGenomeUtils::BSAggregate
# # If we aren't aggregating by bin, then don't restrict to min reads!
# if (byRegionGroup) {
# binnedBSDT <- binnedBSDT[readCount > minReads,]
# }
if (verbose) {
message(".", appendLF=FALSE)
}
return(binnedBSDT)
}
# modification of BSAggregate to just return mean per region
#
# @template signal
# @template signalCoord
# @template regionSet
# Must be from the same reference genome
# as the coordinates for the actual data/samples (signalCoord).
# @template signalCol
# @param returnQuantile "logical" object. If FALSE, return region averages. If TRUE,
# for each region, return the quantile of that region's average value
# based on the distribution of individual genomic signal/feature values
# @template absVal
# @return a data.table with region coordinates and average loading
# values for each region. Has columns chr, start, end, and a column for each
# target variable in signalCol.
# Regions are not in order along the rows of the data.table.
#
# @example averagePerRegion(BSDT = BSDT, regionsGRL,
# jCommand = MIRA:::buildJ(cols = "methylProp", "mean"))
# Devel note: I could add a column for how many cytosines are in each region
averagePerRegion <- function(signal,
signalCoord,
regionSet,
signalCol = c("PC1", "PC2"),
returnQuantile = FALSE,
absVal=TRUE) {
################### checking inputs #################################
########## check that inputs are the correct class and converts
checkConvertInputClasses(signal=signal,
signalCoord=signalCoord,
regionSet=regionSet,
signalCol = signalCol)
########## check that dimensions of inputs are consistent
# length of signal coord = nrow of signal
if (length(signalCoord) != nrow(signal)) {
stop(cleanws("The number of coordinates in
signalCoord (length(signalCoord)) does not equal the number of
rows in `signal`"))
}
######### check that appropriate columns are present
# signalCol are column names of signal
if (!all(signalCol %in% colnames(signal))) {
missingCols = signalCol[!(signalCol %in% colnames(signal))]
stop(cleanws(paste0("Some signalCol are not
columns of signal: ", missingCols)))
}
#######
# what happens if there are NAs or Inf in `signal`?
# any NAs that overlap the regionSet will cause the score to be NA
if (is(signal, "data.table")) {
naRows = apply(X = signal[, signalCol, with=FALSE, drop=FALSE],
MARGIN = 1, FUN = function(x) any(is.na(x)))
} else {
naRows = apply(X = signal[, signalCol, drop=FALSE],
MARGIN = 1, FUN = function(x) any(is.na(x)))
}
if (any(naRows)) {
signal <- signal[!naRows, ]
signalCoord <- signalCoord[!naRows]
warning("Removing rows with NA from `signal`")
}
#################################################################
#################################################################
# determine whether coordinates are single base or a range
# if (!("end" %in% colnames(coordinateDT))) {
# dataCoordType <- "singleBase"
# } else {
if (any(start(signalCoord) != end(signalCoord))) {
dataCoordType <- "multiBase"
} else {
dataCoordType <- "singleBase"
}
# }
# take absolute value or not
if (absVal) {
signalDT <- as.data.table(abs(signal))
} else {
signalDT <- as.data.table(signal)
}
# use different function for single base data and for region data
if (dataCoordType == "singleBase") {
# linking coordinates to loading values, has columns chr start, signalCol
BSDT <- signalDT[, .SD, .SDcols = signalCol]
jExpr <- buildJ(signalCol, rep("mean", length(signalCol)))
avPerRegion <- BSAggregate(BSDT = BSDT,
regionsGRL = regionSet,
BSCoord = signalCoord,
excludeGR = NULL,
regionsGRL.length = NULL,
splitFactor = NULL,
keepCols = NULL,
sumCols = NULL,
jExpr = jExpr,
byRegionGroup = FALSE,
keep.na = FALSE,
returnSD = FALSE,
returnOLInfo = FALSE,
meanPerRegion = TRUE,
returnQuantile = returnQuantile)
} else if (dataCoordType == "multiBase") {
avPerRegion <- weightedAvePerRegion(signalDT = signalDT,
signalCoord=signalCoord,
regionSet=regionSet,
calcCols = signalCol,
returnQuantile = returnQuantile)
} else {
stop("dataCoordType not recognized.")
}
return(avPerRegion)
}
# signalDT
# average by region for region-based data (eg ATAC-seq)
# @return data.table. chr, start, end columns. One column for each calcCols
# that has the average value in each region for that col
weightedAvePerRegion <- function(signalDT,
signalCoord,
regionSet,
calcCols = c("PC1", "PC2"),
returnQuantile = FALSE) {
if (is(signalCoord, "data.frame")) {
signalCoord <- dtToGr(signalCoord)
}
hits <- findOverlaps(query = signalCoord, subject = regionSet)
# if no overlap, return NULL
if (length(hits) == 0) {
return(NULL)
}
olap <- pintersect(signalCoord[queryHits(hits)],
regionSet[subjectHits(hits)])
polap <- width(olap) / width(regionSet[subjectHits(hits)])
# get total proportion overlap per region
# aggregate polap by region
pOlapDT <- data.table(signalDT[queryHits(hits), calcCols, with=FALSE],
rsRegionID = subjectHits(hits),
pOlap = polap)
pOlapByRegionDT <- pOlapDT[, .(regionPOlap = sum(pOlap)), by=rsRegionID]
# specify aggregation operation
# will be done separately for each PC specified
aggrCommand <- paste("list(", paste(paste0(calcCols, "=", "sum", "(",
calcCols, " * pOlap)"), collapse = ", "), ")")
weightedSumByRegionDT <- pOlapDT[, eval(parse(text=aggrCommand)), by=rsRegionID]
# weightedSumByRegionDT and pOlapByRegionDT should be in the same order
regionInd <- weightedSumByRegionDT$rsRegionID
olCoord <- grToDt(regionSet)[regionInd, .(chr, start, end)]
# divide by total proportion overlap to get mean value
jCommand <- paste("list(", paste(paste0(calcCols, "=", calcCols, " / regionPOlap"), collapse = ", "), ")")
meanPerRegion <- cbind(pOlapByRegionDT, weightedSumByRegionDT)
meanPerRegion <- meanPerRegion[, eval(parse(text=jCommand))]
meanPerRegion <- cbind(olCoord, meanPerRegion)
if (returnQuantile) {
for (i in seq_along(calcCols)) {
# perhaps this could be more efficient with mapply
# ecdf example: ecdf(distributionData)(getPercentileOfThis)
meanPerRegion[, c(calcCols[i]) := ecdf(signalDT[[calcCols[i]]])(meanPerRegion[[calcCols[i]]])]
# I tried set() to improve performance but it took about the same time
}
}
# meanPerRegion <- pOlapDT[, .(regionMean = sum(score * (pOlap/sum(pOlap))), by=rsRegionID]
return(meanPerRegion)
}
#' Get regions that are most associated with target variable
#'
#' Get a GRanges with top regions from the region set based on
#' average feature contribution scores
#' for the regions or the quantile of the region's average
#' feature contribution score based on the
#' distribution of all feature contribution scores for the target variable.
#' Returns average feature contribution score or quantile as GRanges metadata.
#'
#' @template signal
#' @template signalCoord
#' @template regionSet
#' @template signalCol
#' @param cutoff Numeric. Only regions with at least this value will be
#' returned (either above this average `signal` value or above this quantile
#' if returnQuantile=TRUE).
#' @param returnQuantile Logical. If FALSE, return region averages. If TRUE,
#' for each region, return the quantile of that region's average value
#' based on the distribution of individual feature values in `signal` for
#' that `signalCol`.
#' @return A GRanges object with region coordinates for regions with
#' scores/quantiles above "cutoff" for any target variable in signalCol.
#' The scores/quantiles
#' for signalCol are given as metadata in the GRanges.
# Are regions in order along the rows of the data.table?
#
#' @examples
#' data("brcaATACCoord1")
#' data("brcaATACData1")
#' data("esr1_chr1")
#' featureContributionScores <- prcomp(t(brcaATACData1))$rotation
#' topRegions <- getTopRegions(signal=featureContributionScores,
#' signalCoord=brcaATACCoord1,
#' regionSet=esr1_chr1,
#' returnQuantile = TRUE)
#' @export
getTopRegions <- function(signal,
signalCoord,
regionSet,
signalCol = c("PC1", "PC2"),
cutoff = 0.8,
returnQuantile=TRUE) {
regionLoadDT <- averagePerRegion(signal=signal,
signalCoord=signalCoord, regionSet=regionSet,
signalCol = signalCol,
returnQuantile = returnQuantile)[]
keepInd <- regionLoadDT[, signalCol, with=FALSE] >= cutoff
# keep region if it is above cutoff in any of the PCs in signalCol
keepInd <- apply(X = keepInd, MARGIN = 1, FUN = any)
highGR <- dtToGr(regionLoadDT[keepInd, ])
values(highGR) <- as.data.frame(regionLoadDT[keepInd, signalCol, with=FALSE])
return(highGR)
}
# different scoring metrics
# remniscent of LOLA:
# support (number of regions, for us regions that have at least 1 cytosine and can be scored)
# mean loading value (or could do ratio of peak to surrounding regions)
# signal to noise ratio, how big is peak compared to noise of surrounding area
# with SNR, even a small peak could have a high SNR
# @param loadingProf
# loadingProfileSNR = function() {
# # magnitude of the peak divided by standard deviation of
# # noise (signal in surrounding areas)
# }
# support is just number of regions from each input region set that overlap at all with
# the cytosines that have loading values
# findOverlaps(regionSet, dtToGr(coordinateDT))
###########################################################################
# BSaggregate -- Aggregate a BSDT across regions or region groups,
# for multiple samples at a time.
# This function is as BScombineByRegion, but can handle not only multiple
# samples in BSDT, but also simultaneously multiple region sets by passing
# a regionsGRL (GRangesList object).
# you can use jExpr to do other functions.
# Given a bisulfite data table as input, with an identifier column for
# different samples; plus a GRanges objects with regions to aggregate.
#
# @param BSDT The bisulfite data.table (output from one of the parsing
# functions for methylation calls) that you wish to aggregate. It can
# be a combined table, with individual samples identified by column passed
# to splitFactor.
# @param regionsGRL Regions across which you want to aggregate. Should be
# from a single region set. eg GRangesList(regionSet)
# @param BSCoord GRanges. Coordinates for BSDT. If NULL, then "chr" and "start"
# columns must be in BSDT.
# @param excludeGR A GRanges object with regions you want to
# exclude from the aggregation function. These regions will be eliminated
# from the input table and not counted.
# @param jExpr You can pass a custom command in the j slot to data.table
# specifying which columns to aggregate, and which functions to use. You
# can use buildJ() to build a jExpr argument easily.
# @param byRegionGroup You can aggregate by regionID or by regionGroupID;
# this reflects the regionsGRL that you pass; by default, BSAggregate will
# aggregate each region individually -- scores will then be contiguous, and
# the output is 1 row per region.
# Turn on this flag to aggregate across all region groups, making the result
# uncontiguous, and resulting in 1 row per *region group*.
# @param returnSD Whether the standard deviation of the columns of interest
# should be returned. Standard deviation for rows that overlap with region set
# and also for rows that do not overlap with region set. Not currently functional.
# @param returnOLInfo if true, include number of overlapping cpgs and number
# of overlapping regions in the result
# @param meanPerRegion Will return the mean value in each region in a data.table
# instead of averaging all regions into a single value, returnOLInfo does
# nothing if meanPerRegion=TRUE (function exits before that code, which
# expects a different data structure)
# @param returnQuantile Only used if meanPerRegion=TRUE, instead of mean
# return the quantile/percentile of the mean of each region
# in relation to the distribution of original values in BSDT
# @template rsOL
#
BSAggregate <- function(BSDT, regionsGRL, BSCoord=NULL, excludeGR=NULL,
regionsGRL.length = NULL, splitFactor=NULL,
keepCols=NULL, sumCols=NULL,
jExpr=NULL, byRegionGroup=FALSE,
keep.na=FALSE, returnSD=FALSE,
returnOLInfo=FALSE, meanPerRegion=FALSE,
returnQuantile=FALSE, rsOL=NULL) {
# Assert that regionsGRL is a GRL.
# If regionsGRL is given as a GRanges, we convert to GRL
if(is(regionsGRL, "GRanges")) {
regionsGRL <- GRangesList(regionsGRL)
} else if (!is(regionsGRL, "GRangesList") && is.null(rsOL)) {
stop("regionsGRL is not a GRanges or GRangesList object")
}
# will cause error if BSDT is only a data.frame
if (is(BSDT, "data.frame") & !is(BSDT, "data.table")) {
BSDT <- as.data.table(BSDT)
}
if (!is(BSDT, "data.table")) {
stop("BSDT must be a data.table")
}
if(! is.null(excludeGR)) {
BSDT <- BSFilter(BSDT, minReads=0, excludeGR)
}
if (returnQuantile) {
# keep all data so quantiles can be calculated
# later code will change BSDT
origBSDT <- data.table::copy(BSDT)
}
# TODO: BSdtToGRanges needs to include end coordinate!!!
if (is.null(BSCoord)) {
bsgr <- BSdtToGRanges(list(BSDT))
} else {
bsgr <- BSCoord
}
additionalColNames <- setdiff(colnames(BSDT),
c("chr","start", "end",
"hitCount","readCount", splitFactor))
# specify that "character" outcome is expected from mode by
# supplying "a" as last vapply argument (any character object would work)
colModes <- vapply(BSDT, mode, "a")
if (is.null(sumCols)) {
sumCols <- setdiff(colnames(BSDT),c("chr", "start", "end",
"strand", splitFactor, keepCols))
# Restrict to numeric columns.
sumCols <- intersect(sumCols,
names(colModes[which(colModes == "numeric")]))
}
if (is.null(rsOL)) {
# It's required to do a findoverlaps on each region individually,
# Not on a GRL, because of the way overlaps with GRLs work. So,
# we must convert the GRL to a GR, but we must keep track of which
# regions came from which group.
regionsGR <- unlist(regionsGRL)
if(is.null(regionsGRL.length)) {
# if (length(regionsGRL) > 100) {
# message("BSAggregate: Calculating sizes. You can speed this up by supplying a regionsGRL.length vector...", appendLF=FALSE)
# }
# specify that output should be numeric with vapply
# (any number would work instead of 1)
regionsGRL.length <- vapply(regionsGRL, length, 1)
# message("Done counting regionsGRL lengths.")
}
# Build a table to keep track of which regions belong to which group
# BIOC note: sapply returns a list where each item is of different length
# therefore, I'm not using vapply
# one regionID per region of "regionsGR", withinGroupID: for when regionsGRL
# has multiple GRanges, Gives the index of each region within its corresponding
# GRanges object. regionGroupID: for when regionsGRL has multiple GRanges,
# indicates which regionsGRL GRanges object the region has come from (
# index of the GRanges object in regionsGRL)
region2group <- data.table(
regionID=seq_along(regionsGR),
chr=as.vector(seqnames(regionsGR)),
start=as.vector(start(regionsGR)),
end=as.vector(end(regionsGR)),
# withinGroupID= as.vector(unlist(sapply(regionsGRL.length, seq))),
regionGroupID=rep(seq_along(regionsGRL), regionsGRL.length)) # repeats each number regionGRL.length times
setkey(region2group, regionID)
#TODO: if overlap method is single, but "end" is present, find center
# and set that to start!
# if (all(start(bsgr[[1]]) != end(bsgr[[length(unique(queryHits(hits)))1]]))) {
if (all(start(bsgr) != end(bsgr))) {
stop("BSDT start and end coordinates are not the same. Choose a different aggrMethod.")
} else {
# fo <- findOverlaps(query = bsgr[[1]], subject = regionsGR)
fo <- findOverlaps(query = bsgr, subject = regionsGR)
}
if (length(subjectHits(fo)) < 1) {
warning("Insufficient overlap between signalCoord and the region set.")
return(NULL)
}
} else {
# Build a table to keep track of which regions belong to which group
# BIOC note: sapply returns a list where each item is of different length
# therefore, I'm not using vapply
# region2group <- data.table(
# regionID=seq_len(regionsGR),
# withinGroupID= as.vector(unlist(sapply(regionsGRL.length, seq))),
# regionGroupID=rep(seq_along(regionsGRL), regionsGRL.length))
if (length(subjectHits(rsOL)) < 1) {
warning("Insufficient overlap between signalCoord and the region set.")
return(NULL)
}
# only works for when one region set is given in regionsGRL (ie does
# not work for metaregion profiles)
region2group <- data.table(
regionID=seq_len(max(subjectHits(rsOL))),
regionGroupID=rep(1, max(subjectHits(rsOL)))) # assumes only 1 region set in regionsGRL
setkey(region2group, regionID)
fo <- rsOL
}
### use info from findOverlaps to see how many individual
# cytosines (or input regions) and how many regions (in region sets)
# overlap with one another
if (returnOLInfo) {
signalCoverage <- length(unique(queryHits(fo)))
regionSetCoverage <- length(unique(subjectHits(fo)))
}
# if("end" %in% colnames(BSDT)){
# setkey(BSDT, chr, start, end)
# } else{
# setkey(BSDT, chr, start)
# }
# Gut check:
# stopifnot(all(elementMetadata(bsgr[[1]])$readCount == BSDT$readCount))
# message("Setting regionIDs...")
BSDT <- BSDT[queryHits(fo),] #restrict the table to CpGs (or input region) in any region set.
BSDT[,regionID:=subjectHits(fo)] #record which region they overlapped.
#BSDT[queryHits(fo),regionID:=subjectHits(fo)]
#if (!keep.na) {
# BSDT = BSDT[queryHits(fo),]
#}
if (is.null(jExpr)) {
cols=c(sumCols, keepCols)
funcs <- c(rep("sum", length(sumCols)), rep("unique", length(keepCols)))
jExpr <- buildJ(cols, funcs)
}
# message("jExpr: ", jExpr)
# Define aggregation column. aggregate by region or by region group?
if (byRegionGroup) {
agCol <- "regionGroupID"
} else {
agCol <- "regionID" # Default
}
# Build the by string
if (is.null(splitFactor)) {
byString <- paste0("list(regionID)")
} else {
byString <- paste0("list(", paste("regionID",
paste0(splitFactor, ""),
collapse=", ", sep=", "), ")")
}
# Now actually do the aggregate: (aggregate within each region)
# message("Combining...")
bsCombined <- BSDT[,eval(parse(text=jExpr)), by=eval(parse(text=byString))]
setkey(bsCombined, regionID)
if (meanPerRegion) {
setkey(region2group, regionID)
avPerRegion <- merge(bsCombined, region2group)
avPerRegion[, c("regionID",
# "withinGroupID",
"regionGroupID") := NULL]
if (returnQuantile) {
for (i in seq_along(sumCols)) {
# perhaps this could be more efficient with mapply
# ecdf example: ecdf(distributionData)(getPercentileOfThis)
avPerRegion[, c(sumCols[i]) := ecdf(origBSDT[[sumCols[i]]])(avPerRegion[[sumCols[i]]])]
# I tried set() to improve performance but it took about the same time
}
}
return(avPerRegion)
}
# Now aggregate across groups.
# I do this in 2 steps to avoid assigning regions to groups,
# which takes awhile. I think this preserves memory and is faster.
# Define aggregation column. aggregate by region or by region group?
# byRegionGroup=TRUE means to aggregate within each individiual GRanges from
# regions.GRL
if (byRegionGroup) {
# must set allow=TRUE here in case there are multiple IDs (splitCol)
bsCombined[region2group, regionGroupID:=regionGroupID, allow=TRUE]
if (! is.null(splitFactor) ) {
byStringGroup <- paste0("list(", paste("regionGroupID",
paste0(splitFactor,
collapse=", "),
sep=", "), ")")
} else {
byStringGroup <- "list(regionGroupID)"
}
# aggregate by regionGroupID (region set/GRanges in regionsGRL)
bsCombined=bsCombined[,eval(parse(text=jExpr)),
by=eval(parse(text=byStringGroup))]
if (returnOLInfo) {
bsCombined[, signalCoverage := signalCoverage]
bsCombined[, regionSetCoverage := regionSetCoverage]
}
return(bsCombined)
} else {
e <- region2group[bsCombined,]
setkey(e, regionID)
return(e)
}
# WARNING: There are now 2^2 ways to aggregate, sum vs mean
# at each level: across regions, then across region sets. This
# doesn't give you a choice at this point.
}
# Given a BSDT (bisulfite data.table), remove any entries that overlap
# regions given in the excludeGR argument and/or filter out sites
# that have lower than a minimum number of reads.
# @param BSDT Bisulfite data.table to filter
# @param minReads Require at least this level of coverage at a cpg.
# @param excludeGR GRanges object with regions to filter.
#
# @return The BSDT with appropriate regions removed.
BSFilter <- function(BSDT, minReads = 10, excludeGR = NULL) {
# First, filter for minimum reads.
if (minReads > 0) {
BSDT <- BSDT[coverage >= minReads, ]
}
if (NROW(BSDT) == 0) { return(data.table(NULL)) }
# Now, filter entries overlapping a region in excludeGR.
if (!is.null(excludeGR)) {
gr <- dtToGr(BSDT)
fo <- findOverlaps(gr, excludeGR)
qh <- unique(queryHits(fo))
length(qh)
nrow(BSDT)
BSDT <- BSDT[-qh, ]
}
return(BSDT)
}
#' Get indices for top scored region sets
#'
#' For each target variable, get index of original region sets
#' but ordered by rsScores
#' ranking for each target variable.
#' The original index refers to that region set's position
#' in the `GRList` param given to `aggregateSignalGRList` which is
#' also that region set's
#' row index in the COCOA output. The first number in a given column
#' of this function's output will be the
#' original index of the region set ranked first for that target variable.
#' The second row for a
#' column will be the original index of the region set that ranked second
#' for that target variable, etc. You can use this function to make it easier
#' when you want to select the top region sets for further analysis or
#' just for sorting the results. Region set scores are sorted in decreasing
#' or increasing order according to the `decreasing` parameter.
#'
#' @template rsScores
#' @templateVar isRSRankingIndex TRUE
#' @template signalCol
#' @param decreasing Logical. Whether to sort rsScores in decreasing
#' or increasing order.
#' @param newColName Character. The names of the columns of the output data.frame.
#' The order should correspond to the order of the
#' input columns given by signalCol.
#' @return A data.frame with one column for each `signalCol`.
#' Column names are given by `signalCol` or `newColName` (if used).
#' Each column has been
#' sorted by score for region sets for that target variable
#' (order given by `decreasing`
#' param).
#' Original indices for region sets that were used to create rsScores
#' are given. Region sets with a score of NA are counted as having the
#' lowest scores and indices for these region sets will be at the bottom of the
#' returned data.frame (na.last=TRUE in sorting)
#' @examples data("rsScores")
#' rsRankInd = rsRankingIndex(rsScores=rsScores,
#' signalCol=c("PC1", "PC2"))
#' # region sets sorted by score for PC1
#' rsScores[rsRankInd$PC1, ]
#' # region sets sorted by score for PC2
#' rsScores[rsRankInd$PC2, ]
#'
#' @export
#'
rsRankingIndex <- function(rsScores, signalCol,
decreasing=TRUE, newColName = signalCol) {
if (!(is(rsScores, "data.frame") || is(rsScores, "matrix"))) {
stop("rsScores should be a data.frame. Check object class.")
}
rsScores <- as.data.table(rsScores)
if (!(is(signalCol, "character") | is(signalCol, "list"))) {
stop("signalCol should be a character vector or list of character vectors.")
}
if (is(signalCol, "list")) {
if (!all(vapply(X = signalCol, FUN = class, FUN.VALUE = "a") %in% "character")) {
stop("Items of signalCol should be character vectors.")
}
if (!length(unique(vapply(signalCol, FUN = length, FUN.VALUE = 2)))) {
stop("Items of signalCol should be the same length as each other.")
}
if (!all(unlist(signalCol) %in% colnames(rsScores))) {
stop("Some column names in signalCol are not present in rsScores.")
}
# the first item in the list is taken for newColName
if (is(newColName, "list")) {
newColName <- signalCol[[1]]
}
} else { # signalCol is character
if (!all(signalCol %in% colnames(rsScores))) {
stop("Some column names in signalCol are not present in rsScores.")
}
}
dtOrder <- rep(-99L, length(decreasing))
# how to sort scores
# -1 for decreasing order of scores
dtOrder[decreasing] <- -1L
# +1 for increasing order of scores
dtOrder[!decreasing] <- 1L
dtOrder <- rep(-99L, length(decreasing))
# how to sort scores
# -1 for decreasing order of scores
dtOrder[decreasing] <- -1L
# +1 for increasing order of scores
dtOrder[!decreasing] <- 1L
# so by references changes will not be a problem
rsScores <- copy(rsScores)
rsScores[, rsIndex := seq_len(nrow(rsScores))]
if (is(signalCol, "list")) {
if (length(newColName) != length(signalCol[[1]])) {
stop("newColName is not the same length as columns given in signalCol.")
}
rsEnSortedInd <- subset(rsScores, select= signalCol[[1]])
setnames(rsEnSortedInd, newColName)
colNameMat <- do.call(rbind, signalCol)
# then scores by each PC and make a column with the original index for sorted region sets
# this object will be used to pull out region sets that were top hits for each PC
for (i in seq_along(signalCol[[1]])) {
theseOrderCols <- colNameMat[, i]
setorderv(rsScores, cols = theseOrderCols, order=dtOrder, na.last=TRUE)
rsEnSortedInd[, newColName[i] := rsScores[, rsIndex]]
}
} else if (is(signalCol, "character")) {
if (length(newColName) != length(signalCol)) {
stop("newColName is not the same length as columns given in signalCol.")
}
signalCol <- signalCol[signalCol %in% colnames(rsScores)]
rsEnSortedInd <- subset(rsScores, select= signalCol)
setnames(rsEnSortedInd, newColName)
# then scores by each PC and make a column with the original index for sorted region sets
# this object will be used to pull out region sets that were top hits for each PC
for (i in seq_along(signalCol)) {
setorderv(rsScores, cols = signalCol[i], order=dtOrder, na.last=TRUE)
rsEnSortedInd[, newColName[i] := rsScores[, rsIndex]]
}
} else {
stop("signalCol should be a character vector or list of character vectors.")
}
# reset order
# setorderv(rsScores, cols = "rsIndex", order=1L)
return(as.data.frame(rsEnSortedInd))
}
#################### Metric functions ########################################
# scores, metrics, or statistical tests
# Instead of averaging within regions first as BSAggregate does,
# this function does a simple average and standard deviation
# for all CpGs that overlap
# with regions of a region set, also does average and
# standard deviation for non overlapping CpGs. Created to
# get metrics of loading values for each PC.
#
# Faster if given total average for each column of interest
#
# @param dataDT a data.table with
# columns to get metrics of eg (PC1, PC2). All columns
# will be considered
# columns to get the metrics from so no unnecessary columns should be
# included.
# @param dataGR GRanges. Coordinates for dataDT.
# @template regionSet
# Metrics will be calculated on
# only coordinates within this region set (and optionally separately
# on those outside this region set with alsoNonOLMet parameter)
# @param signalCol the columns to calculate the metrics on. The
# metrics will be calculated on each one of these columns separately.
# @param metrics character vector with the name of a function or functions
# to calculate on selected cytosines. Function should only require one
# input which should be the values of the cytosines.
# @param alsoNonOLMet also include same metrics
# for non overlapping CpGs (still also returns metrics for overlapping CpGs)
# param columnMeans Not functional/is deprecated. The idea was to use
# the mean of the column to speed up calculations (eg when calculating
# mean of overlapping CpGs, use that info and column mean to get
# mean for non overlapping CpGs without manually calculating it)
#
signalOLMetrics <- function(dataDT,
dataGR,
regionSet,
signalCol = colnames(dataDT)[!(colnames(dataDT) %in% c("chr", "start", "end"))],
metrics=c("mean", "sd"),
alsoNonOLMet=TRUE, rsOL=NULL) {
# convert DT to GR for finding overlaps
# dataGR <- BSdtToGRanges(list(dataDT))[[1]]
if (is.null(rsOL)) {
OL <- findOverlaps(query = dataGR, subject = regionSet)
# region set info
totalRegionNumber <- length(regionSet)
meanRegionSize <- round(mean(width(regionSet)), 1)
} else {
OL <- rsOL
}
# if no overlap, exit
if (length(OL) == 0) {
return(NULL)
}
# get indices for overlapping and non overlapping CpGs
olCpG <- queryHits(OL)
# get info on degree of overlap
# number of CpGs that overlap
signalCoverage <- length(unique(olCpG))
# number of regions that overlap
regionSetCoverage <- length(unique(subjectHits(OL)))
# gets metrics for all columns except chr, start, end
jExpr <- buildJ(cols=rep(signalCol, each=length(metrics)),
funcs=rep(metrics, length(signalCol)),
newColNames = paste0(rep(signalCol,
each=length(metrics)),
"_", metrics))
# getting the metrics
olMetrics <- as.data.frame(dataDT[olCpG, eval(parse(text=jExpr))])
# calculate average of nonOLCpGs based on columnMean if given
# if (!is.null())
#
# formatting so there is one row per PC/testCol
# output is a matrix with ncol = length(metrics)
# for vapply, FUN.VALUE should have length equal to a single output of FUN
olResults <- vapply(X = metrics,
FUN = function(x) as.numeric(olMetrics[, grepl(pattern = x, colnames(olMetrics))]),
as.numeric(seq_along(signalCol)))
olResults <- as.data.table(olResults)
setnames(olResults, old = colnames(olResults), new = paste0(colnames(olResults), "_OL"))
if (alsoNonOLMet) {
nonOLCpG <- (seq_len(nrow(dataDT)))[-olCpG]
# if no OL for region set, don't calculate for non region set
# TODO make conditional on region set having any overlap
nonOLMetrics <- as.data.frame(dataDT[nonOLCpG, eval(parse(text=jExpr))])
nonOLResults <- vapply(X = metrics,
FUN = function(x) as.numeric(nonOLMetrics[, grepl(pattern = x, colnames(nonOLMetrics))]),
as.numeric(seq_along(signalCol)))
nonOLResults <- as.data.table(nonOLResults)
setnames(nonOLResults, old = colnames(nonOLResults), new = paste0(colnames(nonOLResults), "_nonOL"))
if (is.null(rsOL)) {
metricDT <- cbind(data.table(testCol=signalCol),
olResults,
nonOLResults,
data.table(signalCoverage,
regionSetCoverage,
totalRegionNumber,
meanRegionSize))
} else {
metricDT <- cbind(data.table(testCol=signalCol),
olResults,
nonOLResults,
data.table(signalCoverage,
regionSetCoverage))
}
} else {
if (is.null(rsOL)) {
metricDT <- cbind(data.table(testCol=signalCol),
olResults,
data.table(signalCoverage,
regionSetCoverage,
totalRegionNumber,
meanRegionSize))
} else {
metricDT <- cbind(data.table(testCol=signalCol),
olResults,
data.table(signalCoverage,
regionSetCoverage))
}
}
return(metricDT)
}
# Wilcoxon rank sum test for a region set
# @param dataDT a data.table with chr, start, end columns as well
# as columns to get metrics of eg (PC1, PC2). All columns
# except chr, start, and end will be considered
# columns to get the metrics from so no unnecessary columns should be
# included.
# @template regionSet
# @param signalCol the columns of interest. You will do ranksum test separately
# on each of these columns (given test only uses info in one column)
# @param ... Additional parameters of wilcox.test function. See ?wilcox.test.
# For instance specify alternative hypothesis: alternative = "greater".
# @return A vector with a p value for each column other than chr, start or end.
# @examples data("brcaLoadings1")
# data("brcaMCoord1")
# data("nrf1_chr1")
# dataDT = as.data.table(cbind(brcaMCoord1, brcaLoadings1))
# rsWilcox(dataDT = dataDT, regionSet = nrf1_chr1, conf.int=TRUE)
rsWilcox <- function(dataDT,
regionSet,
signalCol = colnames(dataDT)[!(colnames(dataDT) %in% c("chr", "start", "end"))],
conf.int=FALSE,
...) {
# region set info
totalRegionNumber <- length(regionSet)
meanRegionSize <- round(mean(width(regionSet)), 1)
# convert DT to GR for finding overlaps
dataGR <- BSdtToGRanges(list(dataDT))[[1]]
OL <- findOverlaps(query = regionSet, subject = dataGR)
# if no overlap, exit
if (length(OL) == 0) {
return(NULL)
}
# get indices for overlapping and non overlapping CpGs
olCpG <- unique(subjectHits(OL))
nonOLCpG <- (seq_len(nrow(dataDT)))[-olCpG]
# get info on degree of overlap
# number of CpGs that overlap
signalCoverage <- length(unique(olCpG))
# number of regions that overlap
regionSetCoverage <- length(unique(queryHits(OL)))
# each confidence interval has length of 2: [low, high]
if (conf.int) {
# calculate Wilcoxon rank sum test for each column
# additional parameters given with ...
# confIntervals will be [low1, high1, low2, high2, etc.]
confIntervals <- as.numeric(vapply(X = signalCol, FUN = function(x) wilcox.test(x = as.numeric(as.matrix(dataDT[olCpG, x, with=FALSE])),
y = as.numeric(as.matrix(dataDT[nonOLCpG, x, with=FALSE])),
conf.int = conf.int, ...)$conf.int, c(1, 1)))
names(confIntervals) <- paste0(rep(signalCol, each=2), c("_low", "_high"))
wRes <- data.frame(t(confIntervals),
signalCoverage,
regionSetCoverage,
totalRegionNumber,
meanRegionSize)
} else {
# calculate Wilcoxon rank sum test for each column
# additional parameters given with ...
pVals <- vapply(X = signalCol, FUN = function(x) wilcox.test(x = as.numeric(as.matrix(dataDT[olCpG, x, with=FALSE])),
y = as.numeric(as.matrix(dataDT[nonOLCpG, x, with=FALSE])), ...)$p.value, 1)
wRes <- data.frame(t(pVals),
signalCoverage,
regionSetCoverage,
totalRegionNumber,
meanRegionSize)
}
return(wRes)
}
# I optimized upstream so that a matrix would be given to this function
# if this function is rewritten and no longer requires a matrix input,
# then in order to prevent unnecessary object copying,
# rewrite upstream code that converts signalDT to matrix class
# @param signalMat Data to be aggregated (e.g. raw data: ATAC-seq,
# region based DNA methylation or loading values)
# @param signalGR GRanges object with coordinates for signalMat
# @template regionSet
# The region set to score.
# @param calcCols character object. Column names. A weighted sum will be done
# for each of these columns (columns should be numeric).
# @template rsOL
# @param pOlap see "?aggregateSignal"
# @template returnCovInfo
# @value Returns data.frame with columns 'calcCols', signalCoverage col has
# number of signalGR regions that overlapped with any regionSet regions,
# regionSetCoverage has the sum of all proportion overlaps of regions from
# signalGR with regionSet (regionSet region is denominator)
# containing weighted mean for each col.
# Returns NULL if there is no overlap between signalGR and regionSet
regionOLWeightedMean <- function(signalMat, signalGR,
regionSet, calcCols, rsOL=NULL,
pOlap=NULL, returnCovInfo=TRUE) {
if (!is(signalMat, "matrix")) {
signalMat <- as.matrix(signalMat)
}
if (is.null(rsOL)) {
hits <- findOverlaps(query = signalGR, subject = regionSet)
} else {
hits <- rsOL
}
# if no overlap, return NULL
if (length(hits) == 0) {
return(NULL)
}
if (is.null(pOlap)) {
olap <- pintersect(signalGR[queryHits(hits)],
regionSet[subjectHits(hits)])
pOlap <- width(olap) / width(regionSet[subjectHits(hits)])
}
# some rows may be duplicated if a signalMat region overlapped multiple
# regions from signalGR but that is ok
# signalMat <- signalMat[queryHits(hits), ] # done in next step to prevent extra copying
# weight the signalMat values by the proportion overlap (weighted average)
weightedSum <- t(pOlap) %*% signalMat[queryHits(hits), calcCols]
# weighted average
denom <- sum(pOlap)
weightedAve <- as.data.frame(weightedSum / denom)
colnames(weightedAve) <- calcCols
# add columns for coverage info
if (returnCovInfo) {
weightedAve$signalCoverage = length(unique(queryHits(hits)))
weightedAve$regionSetCoverage = length(unique(subjectHits(hits)))
weightedAve$sumProportionOverlap = denom
}
return(weightedAve)
}
# @param signalDT Data to be aggregated (e.g. raw data: ATAC-seq,
# region based DNA methylation or loading values)
# @param signalGR GRanges object with coordinates for signalDT
# @template regionSet
# The region set to score.
# @param calcCols character object. Column names. A mean will be calculated for
# each of these columns (columns should be numeric).
# @param metric character. "mean" or "median"
# @template rsOL
# @template returnCovInfo
# @value Returns data.frame with columns 'calcCols', signalCoverage col has
# number of signalGR regions that overlapped with any regionSet regions,
# regionSetCoverage has the number of regions from
# signalGR that overlapped with regionSet
# Returns NULL if there is no overlap between signalGR and regionSet
regionOLMean <- function(signalDT, signalGR, regionSet,
calcCols, metric="mean", rsOL=NULL, returnCovInfo=TRUE) {
if (is.null(rsOL)) {
hits <- findOverlaps(query = signalGR, subject = regionSet)
} else {
hits <- rsOL
}
# if no overlap, return NULL
if (length(hits) == 0) {
return(NULL)
}
# some rows may be duplicated if a signalDT region overlapped multiple
# regions from signalGR but that is ok
if (metric == "mean") {
# mean of the overlapping signalDT values
signalAve <- as.data.frame(t(colMeans(signalDT[queryHits(hits),..calcCols])))
} else if (metric == "median") {
# median of the overlapping signalDT values
signalAve <- as.data.frame(t(apply(X = signalDT[queryHits(hits),..calcCols], 2, median)))
} else {
stop("Error in regionOLMean function. Invalid metric specified.")
}
if (returnCovInfo) {
# add columns for coverage info
signalAve$signalCoverage <- length(unique(queryHits(hits)))
signalAve$regionSetCoverage <- length(unique(subjectHits(hits)))
}
return(signalAve)
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.