R/workflow.R
In ZarnackGroup/BindingSiteFinder: Binding site defintion based on iCLIP data
Defines functions
assignToGenes
pureClipGeneWiseFilter
pureClipGlobalFilter
Documented in
assignToGenes
pureClipGeneWiseFilter
pureClipGlobalFilter
################################################################################
# Workflow functions
# - these functions are used in the main bindingSiteFinde workflow
# - each function updates the bindingSiteFinder object
# - functions are either designed to work on single nt wide crosslink sites or
# on merged binding sites
# - for each function a related plotting funciton exists to visualize the results
# - the wrapper function BSFind() executes the functions below in the desired order
################################################################################
#' Filter PureCLIP sites by their score distribution
#'
#' Function that applies a filter on the global crosslink site score distribution.
#' The \code{\link{GenomicRanges}} contained in the \code{\link{BSFDataSet}} need to
#' have a meta-column that holds a numeric score value, which is used for filtering.
#' The name of the column can be set with \code{match.score}.
#'
#' The function is part of the standard workflow performed by \code{\link{BSFind}}.
#'
#' @param object a \code{\link{BSFDataSet}} object with stored crosslink ranges
#' of width=1
#' @param cutoff numeric; defines the cutoff for which sites to keep, the
#' smallest step is 1\% (0.01)
#' @param match.score character; meta column name of the crosslink site
#' \code{\link{GenomicRanges}} object that holds the score which is used for
#' sub-setting
#' @param quiet logical; whether to print messages
#'
#' @return an object of class \code{\link{BSFDataSet}} with its ranges filtered
#' by those that passed the threshold set with \code{cutoff}
#'
#' @seealso \code{\link{BSFind}}, \code{\link{pureClipGlobalFilterPlot}}
#'
#' @examples
#' # load data
#' files <- system.file("extdata", package="BindingSiteFinder")
#' load(list.files(files, pattern = ".rda$", full.names = TRUE))
#' # apply 5% filter
#' pureClipGlobalFilter(object = bds, cutoff = 0.05)
#'
#' @export
pureClipGlobalFilter <- function(object, # bindingSiteFinder
cutoff = 0.01, # defines the cutoff for which PureCLIP sites to keep; smallest are 1% (0.01); a cutoff of 0.05 means to remove the lowest 5% PureCLIP sites base on the score
match.score = "score",
quiet = FALSE
) {
# INPUT CHECKS
# --------------------------------------------------------------------------
stopifnot(is(object, "BSFDataSet"))
stopifnot(is.logical(quiet))
# check for correct cutoff value
possibleCutoffs = seq(from = 0, to = 1, by = 0.01)
if (! cutoff %in% possibleCutoffs) {
msg = paste0("Specified cutoff (", cutoff,
"), does not refer to a percentage value. Try cutoffs of the form 0.01, 0.05, ... instead.")
stop(msg)
}
# Check for correct score column
metaNames = colnames(mcols(getRanges(object)))
if (! match.score %in% metaNames) {
msg = paste0("Specified match.score (", match.score, "), is present in the ranges. Add respective column. ")
stop(msg)
}
# Create matching vectors for score column from input source
# --------------------------------------------------------------------------
rngInitial = getRanges(object)
score = mcols(rngInitial)[match(match.score, colnames(mcols(rngInitial)))][[1]]
# ---
# Store function parameters in list
optstr = list(cutoff = cutoff, match.score = match.score)
object@params$pureClipGlobalFilter = optstr
# MAIN COMPUTE
# --------------------------------------------------------------------------
# compute quantile based on PureCLIP score
qCut = quantile(score, probs = seq(0,1, by = 0.01))
# apply user set cutoff
names(qCut) = as.numeric(sub("%", "", names(qCut))) / 100
idx = match(cutoff, names(qCut))
rngFiltered = rngInitial[score >= qCut[idx]]
# ---
# Store for plotting
dfPlot = score
cutoffPoint = qCut[idx]
object@plotData$pureClipGlobalFilter$data = dfPlot
object@plotData$pureClipGlobalFilter$cutoffPoint = cutoffPoint
# ---
# Store for results
resultLine = data.frame(
funName = "pureClipGlobalFilter()", class = "crosslink sites",
nIn = length(rngInitial), nOut = length(rngFiltered),
per = paste0(round(length(rngFiltered)/ length(rngInitial), digits = 2)*100,"%"),
options = paste0("Cutoff=", optstr$cutoff, ", scoreMatchCol=", optstr$match.score)
)
object@results = rbind(object@results, resultLine)
# return BSF object with only ranges above cutoff
objectNew = setRanges(object, rngFiltered, quiet = quiet)
return(objectNew)
}
#' Filter PureCLIP sites by their score distribution per gene
#'
#' Function that applies a filter on the crosslink site score distribution at
#' gene level. This allows to filter for those sites with the strongest signal
#' on each gene. Since scores are tied to the expression level of the hosting
#' transcript this function allows a fair filter for all genes partially
#' independent of the expression level.
#'
#' The \code{\link{GenomicRanges}} contained in the \code{\link{BSFDataSet}} need to
#' have a meta-column that holds a numeric score value, which is used for filtering.
#' The name of the column can be set with \code{scoreCol}.
#'
#' In the case of overlapping gene annotation, a single crosslink site will be
#' attributed to multiple genes. The \code{\link{overlaps}} parameter allows
#' to control these cases. Option `keepSingle` will only keep a single instance
#' of the site; `removeAll` will remove both sites; `keepAll` will keep both
#' sites.
#'
#' The function is part of the standard workflow performed by \code{\link{BSFind}}.
#'
#' @param object a \code{\link{BSFDataSet}} object with stored crosslink ranges
#' of width=1
#' @param cutoff numeric; defines the cutoff for which sites to remove, the
#' smallest step is 1\% (0.01). A cutoff of 5\% will remove the lowest 5\% sites,
#' given their score, on each gene, thus keeping the strongest 95\%.
#' @param overlaps character; how overlapping gene loci should be handled.
#' @param anno.annoDB an object of class \code{OrganismDbi} that contains
#' the gene annotation (!!! Experimental !!!).
#' @param anno.genes an object of class \code{\link{GenomicRanges}} that represents
#' the gene ranges directly
#' @param match.score character; meta column name of the crosslink site
#' \code{\link{GenomicRanges}} object that holds the score which is used for
#' sub-setting
#' @param match.geneID character; meta column name of the genes
#' \code{\link{GenomicRanges}} object that holds a unique geneID
#' @param quiet logical; whether to print messages
#'
#' @return an object of class \code{\link{BSFDataSet}} with its ranges filtered
#' by those that passed the gene-wise threshold set with \code{cutoff}
#'
#' @seealso \code{\link{BSFind}}, \code{\link{estimateBsWidthPlot}}
#'
#' @importFrom dplyr filter reframe group_by mutate mutate_all arrange row_number
#' @importFrom GenomicFeatures genes
#' @importFrom stats density
#' @importFrom grDevices colorRampPalette
#'
#' @examples
#' # load clip data
#' files <- system.file("extdata", package="BindingSiteFinder")
#' load(list.files(files, pattern = ".rda$", full.names = TRUE))
#' # Load GRanges with genes
#' load(list.files(files, pattern = ".rds$", full.names = TRUE)[1])
#' # apply 5% gene-wise filter
#' pureClipGeneWiseFilter(object = bds, anno.genes = gns, cutoff = 0.5, overlaps = "keepSingle")
#'
#' @export
pureClipGeneWiseFilter <- function(object, # bindingSiteFinder
cutoff = 0.05, # defines the cutoff for which PureCLIP sites to keep; smallest steps are in 1% (0.01); cutoff of 0.7 means that the top 70% of PureCLIP sites per gene are retained.
overlaps = c("keepSingle", "removeAll", "keepAll"), # Options to deal with PureCLIP sites on overlapping loci genes; use keepSingle as default; keepAll will inflate the number of binding sites by the number of overlaps, removeAll will kill all duplicated instances
anno.annoDB = NULL, # annotation data base of class OrganismDbi; can be NULL -> but then anno.genes must be provided
anno.genes = NULL, # ranges of the genes of class GenomicRanges; can be NULL -> but then anno.annoDB must be provided
match.score = "score",
match.geneID = "gene_id",
quiet = FALSE) {
# Local variables
datasource <- cutoffs <- quantiles <- name <- NULL
# INPUT CHECKS
# --------------------------------------------------------------------------
# type checks
stopifnot(is(object, "BSFDataSet"))
stopifnot(is.logical(quiet))
# Check if none is specified
if (is.null(anno.annoDB) & is.null(anno.genes)) {
msg = paste0("None of the required annotation sources anno.annoDB or anno.genes was specified. ")
stop(msg)
}
# Check if both are specified
if (!is.null(anno.annoDB) & !is.null(anno.genes)) {
msg = paste0("Both of the required annotation sources anno.annoDB or anno.genes are specified. Please provide only one of the two. ")
stop(msg)
}
# Checks if anno.annoDB should be used
if (!is.null(anno.annoDB) & is.null(anno.genes)) {
stopifnot(is(anno.annoDB, "OrganismDb"))
datasource = "anno.annoDB"
# extract relevant annotation
anno.genes = genes(anno.annoDB, columns = c("GENEID"))
# Create matching vectors for columns from input annotation
selectID = as.character(anno.genes$GENEID)
}
# Checks if anno.genes should be used
if (is.null(anno.annoDB) & !is.null(anno.genes)) {
stopifnot(is(anno.genes, "GenomicRanges"))
datasource = "anno.genes"
# extract relevant annotation
anno.genes = anno.genes
# check correct annotation columns
annoColNames = colnames(mcols(anno.genes))
if (!all(match.geneID %in% annoColNames)) {
msg = paste0("The provided matching column for the geneID (",
match.geneID, ") is not present in the provided annotation. \n")
stop(msg)
}
# Create matching vectors for columns from input annotation
selectID = mcols(anno.genes)[match(match.geneID, colnames(mcols(anno.genes)))][[1]]
}
# check if cutoff was estimated before or if default should be used
if (!is.null(object@params$bsSize) &
!is.null(object@params$geneFilter)) {
cutoff = object@params$geneFilter
}
# handle options (remove, keep)
overlaps = match.arg(overlaps, choices = c("keepSingle", "removeAll", "keepAll"))
# Check for correct score column
# get all ranges
rngInitial = getRanges(object)
metaNames = colnames(mcols(rngInitial))
if (! match.score %in% metaNames) {
# if (! any(metaNames == match.score)) {
msg = paste0("Specified match.score (", match.score, "), is present in the ranges. Add respective column. ")
stop(msg)
}
# ---
# Store function parameters in list
optstr = list(cutoff = cutoff, overlaps = overlaps,
source = datasource, match.geneID = match.geneID)
object@params$geneWiseFilter = optstr
# MAIN COMPUTE
# --------------------------------------------------------------------------
# Create matching vectors for score column from input source
# --------------------------------------------------------------------------
score = mcols(rngInitial)[match(match.score, colnames(mcols(rngInitial)))][[1]]
rngInitial$score = score
# get all intergenic cases
rngIntergeneic = rngInitial[countOverlaps(rngInitial, anno.genes) == 0]
mcols(rngIntergeneic)$geneID = "Intergenic"
# get all cases to solve
# rngToSolve = subset(rngInitial, !rngInitial %in% rngIntergeneic)
# rngToSolve = rngInitial[!rngInitial %in% rngIntergeneic]
rngToSolve = rngInitial[countOverlaps(rngInitial, rngIntergeneic) == 0]
# Gene specific cutoffs
# ---------------------
# set cutpoints
potentialCutoffs = seq(0, 1, by = 0.01)
# group scores by gene range
ols = findOverlaps(anno.genes, rngToSolve) %>%
as.data.frame() %>%
mutate(score = rngToSolve$score[subjectHits]) %>%
group_by(queryHits)
# calculate gene range specific quantile cutoff
quants = ols %>%
dplyr::reframe(quantiles = quantile(score, probs = potentialCutoffs)) %>%
group_by(queryHits) %>%
mutate(cutoffs = potentialCutoffs) %>%
dplyr::filter(as.character(cutoffs) == as.character(cutoff))
# apply cutoff to every PureCLIP peak per gene
peakIdToKeep = ols
peakIdToKeep$quantiles = quants$quantiles[match(peakIdToKeep$queryHits, quants$queryHits)]
peakIdToKeep = peakIdToKeep %>%
dplyr::filter(score >= quantiles)
# add gene information
filteredPerRegion = rngToSolve[peakIdToKeep$subjectHits]
mcols(filteredPerRegion)$geneID = selectID[peakIdToKeep$queryHits]
# ---
# Store for plotting
dfPlot = countOverlaps(filteredPerRegion) %>% table() %>% as.data.frame() %>% rename("#N overlaps" = ".")
object@plotData$geneWiseFilter$data = dfPlot
# Deal with multiple loci
# -----------------------
# check for overlapping loci peaks
potentialOverlaps = findOverlaps(filteredPerRegion)
overlappingLociHits = potentialOverlaps[queryHits(potentialOverlaps) != subjectHits(potentialOverlaps)]
# reduce to a single pair
singlePair = overlappingLociHits %>%
as.data.frame() %>%
dplyr::filter(queryHits > subjectHits)
overlappingLociHits = overlappingLociHits[queryHits(overlappingLociHits) %in% singlePair$queryHits]
# get numbers
overlappingLociGenePairs = paste0(unique(names(filteredPerRegion[queryHits(overlappingLociHits)])),
",",
unique(names(filteredPerRegion[subjectHits(overlappingLociHits)])), " | ")
fractionOfRangesRemoved = paste0(round(length(overlappingLociHits) / length(rngInitial) * 100, digits = 2), "%")
# Handle cases
if (length(overlappingLociHits) > 0) {
msg0 = paste0(fractionOfRangesRemoved, " (", length(overlappingLociHits),"/", length(rngInitial), ")",
" peaks overlap with multiple anno.genes in the given gene annotation. \n")
if (overlaps == "keepSingle") {
msg1 = "A single instance of each peak is kept. This is recommended. \n "
filteredPerRegion = filteredPerRegion[!duplicated(filteredPerRegion)]
# report message
if (!quiet) message(c(msg0, msg1))
}
if (overlaps == "removeAll") {
msg1 = "All peaks on any duplicated range are removed. This is not recommended, since many peaks are lost. \n"
filteredPerRegion = filteredPerRegion[countOverlaps(filteredPerRegion) == 1]
warning(c(msg0, msg1))
}
if (overlaps == "keepAll") {
msg1 = "Duplicated peaks due to overlapping loci are not removed. This is not recommended, since peak numbers are inflated. \n"
warning(msg0, msg1)
}
}
# add intergenic cases again
filteredPerRegion = c(filteredPerRegion, rngIntergeneic)
# ---
# Store for results
resultLine = data.frame(
funName = "pureClipGeneWiseFilter()", class = "crosslink sites",
nIn = length(rngInitial), nOut = length(filteredPerRegion),
per = paste0(round(length(filteredPerRegion)/ length(rngInitial), digits = 2)*100,"%"),
options = paste0("Cutoff=", optstr$cutoff, ", overlapp option=", optstr$overlaps)
)
object@results = rbind(object@results, resultLine)
objectOut = setRanges(object, filteredPerRegion, quiet = quiet)
return(objectOut)
}
#' Assign binding sites to their hosting genes
#'
#' Function that assigns each binding site in the \code{\link{BSFDataSet}} to its
#' hosting gene given a gene annotation (\code{anno.annoDB}, \code{anno.genes}).
#'
#' Regardless of the annotation source that is being used, the respective meta
#' information about the genes have to be present. They can be set by the
#' \code{match.geneID}, \code{match.geneName} and \code{match.geneType} arguments.
#'
#' In the case of overlapping gene annotation, a single binding site will be
#' associated with multiple genes. The \code{\link{overlaps}} parameter allows
#' to decide in these cases. Option `frequency` will take the most frequently
#' observed gene type, option `hierarchy` works in conjunction with a user defined
#' rule (\code{overlaps.rule}). Options `remove` and `keep` will remove or
#' keep all overlapping cases, respectively.
#'
#' Note that if an overlaps exists, but gene types are identical options
#' `frequency` and `hierarchy` will cause the gene that was seen first to be
#' selected as representative.
#'
#' The function is part of the standard workflow performed by \code{\link{BSFind}}.
#'
#' @param object a \code{\link{BSFDataSet}} object with stored binding sites. This
#' means that ranges should be > 1
#' @param overlaps character; how overlapping gene loci should be handled.
#' @param overlaps.rule character vector; a vector of gene type that should
#' be used to handle overlapping cases in a hierarchical manor. The order of the
#' vector is the order of the hierarchy.
#' @param anno.annoDB an object of class \code{OrganismDbi} that contains
#' the gene annotation (!!! Experimental !!!).
#' @param anno.genes an object of class \code{\link{GenomicRanges}} that represents
#' the gene ranges directly
#' @param match.geneID character; meta column name of the gene ID
#' @param match.geneName character; meta column name of the gene name
#' @param match.geneType character; meta column name of the gene type
#' @param quiet logical; whether to print messages
#'
#' @return an object of class \code{\link{BSFDataSet}} with binding sites having
#' hosting gene information added to their meta columns.
#'
#' @seealso \code{\link{BSFind}}, \code{\link{geneOverlapsPlot}},
#' \code{\link{targetGeneSpectrumPlot}}
#'
#' @importFrom dplyr select slice_head group_by count ungroup arrange desc left_join
#' @importFrom GenomicFeatures genes
#' @importFrom dplyr desc
#'
#' @examples
#' # load clip data
#' files <- system.file("extdata", package="BindingSiteFinder")
#' load(list.files(files, pattern = ".rda$", full.names = TRUE))
#' # Load GRanges with genes
#' load(list.files(files, pattern = ".rds$", full.names = TRUE)[1])
#' bds = makeBindingSites(object = bds, bsSize = 9)
#' bds = assignToGenes(bds, anno.genes = gns)
#'
#' @export
assignToGenes <- function(object,
overlaps = c("frequency", "hierarchy", "remove", "keep"),
overlaps.rule = NULL,
anno.annoDB = NULL,
anno.genes = NULL,
match.geneID = "gene_id",
match.geneName = "gene_name",
match.geneType = "gene_type",
quiet = FALSE
) {
# local variables
datasource <- geneIndex <- geneID <- geneName <- value <- bsIndex <- Freq <- geneType <- match.missing <- NULL
selectType <- selectID <- selectName <- NULL
# INPUT CHECKS
# --------------------------------------------------------------------------
stopifnot(is(object, "BSFDataSet"))
stopifnot(is.logical(quiet))
# Check if none is specified
if (is.null(anno.annoDB) & is.null(anno.genes)) {
msg = paste0("None of the required annotation sources anno.annoDB or anno.genes was specified. ")
stop(msg)
}
# Check if both are specified
if (!is.null(anno.annoDB) & !is.null(anno.genes)) {
msg = paste0("Both of the required annotation sources anno.annoDB or anno.genes are specified. Please provide only one of the two. ")
stop(msg)
}
# Checks if anno.annoDB should be used
if (!is.null(anno.annoDB) & is.null(anno.genes)) {
stopifnot(is(anno.annoDB, "OrganismDb"))
if (!is.null(anno.genes)) {
msg = paste0("Parameter anno.annoDB and anno.genes are specified at the same time. Use only one of them.")
stop(msg)
} else {
datasource = "anno.annoDB"
# extract relevant annotation
anno.genes = genes(anno.annoDB, columns=c("ENSEMBL", "GENETYPE", "SYMBOL", "GENEID"))
# Create matching vectors for columns from input annotation
# --------------------------------------------------------------------------
selectID = as.character(anno.genes$GENEID)
selectName = as.character(anno.genes$SYMBOL)
selectType = as.character(anno.genes$GENETYPE)
}
}
# Checks if anno.genes should be used
if (is.null(anno.annoDB) & !is.null(anno.genes)) {
stopifnot(is(anno.genes, "GenomicRanges"))
if (!is.null(anno.annoDB)) {
msg = paste0("Parameter anno.annoDB and anno.genes are specified at the same time. Use only one of them.")
stop(msg)
} else {
datasource = "anno.genes"
# extract relevant annotation
# check for duplicated annotation columns
if (sum(colnames(mcols(anno.genes)) == match.geneID) > 1) {
msg = paste0("The names of multiple columns of the annotation match the match.geneID parameter. Please use a unique column name for match.geneID in anno.genes.")
stop(msg)
}
if (sum(colnames(mcols(anno.genes)) == match.geneName) > 1) {
msg = paste0("The names of multiple columns of the annotation match the match.geneName parameter. Please use a unique column name for match.geneName in anno.genes.")
stop(msg)
}
if (sum(colnames(mcols(anno.genes)) == match.geneType) > 1) {
msg = paste0("The names of multiple columns of the annotation match the match.geneType parameter. Please use a unique column name for match.geneType in anno.genes.")
stop(msg)
}
# check correct annotation columns
inNames = c(match.geneID, match.geneName, match.geneType)
annoColNames = colnames(mcols(anno.genes))
presentNames = inNames[inNames %in% annoColNames]
# extract columns from annotaiton based on available names
present.cols = lapply(presentNames, function(x){
match.col = mcols(anno.genes)[match(x, colnames(mcols(anno.genes)))][[1]]
return(match.col)
})
names(present.cols) = presentNames
# extract values from annotation for all present columns
if (match.geneID %in% names(present.cols)) {
selectID = present.cols[[match(match.geneID, names(present.cols))]]
} else {
msg = paste0("No meta column for ", match.geneID, " present. Creating custom ID.\n")
if (!quiet) message(msg)
selectID = paste0("CUSTOM", seq_along(anno.genes))
}
if (match.geneName %in% names(present.cols)) {
selectName = present.cols[[match(match.geneName, names(present.cols))]]
} else {
msg = paste0("No meta column for ", match.geneName, " present.\n")
if (!quiet) message(msg)
}
if (match.geneType %in% names(present.cols)) {
selectType = present.cols[[match(match.geneType, names(present.cols))]]
} else {
msg = paste0("No meta column for ", match.geneType, " present.\n")
if (!quiet) message(msg)
}
}
}
# handle options (hierarchy, frequency, remove, keep)
overlaps = match.arg(overlaps, choices = c("frequency", "hierarchy", "remove", "keep"))
if (is.null(selectType)) {
# gene type is missing
# -> check for correct overlap option
if (overlaps == "frequency" | overlaps == "hierarchy") {
msg0 = paste0(match.geneType, " is missing in the annotation. Options 'frequency' and 'hierarchy' are not available.\n")
msg1 = paste0("Select one of: 'remove' or 'keep'.\n")
stop(c(msg0, msg1))
}
}
# Check multiple loci options
if (overlaps == "hierarchy" & is.null(overlaps.rule)) {
msg1 = paste0("Binding sites on anno.genes with overlapping annotaitons is set to be handled by 'hierarchy', but no rule is provided. \n")
msg2 = paste0("Change the 'overlaps' option or provide a valid rule. \n")
stop(c(msg1, msg2))
}
# --------------------------------------------------------------------------
# MAIN
# --------------------------------------------------------------------------
# ---
# Store function parameters in list
optstr = list(source = datasource, match.geneID = match.geneID,
match.geneName = match.geneName, match.geneType = match.geneType,
overlaps = overlaps, overlaps.rule = overlaps.rule)
object@params$assignToGenes = optstr
# Match Binding sites and anno.genes
# --------------------------------------------------------------------------
# split ranges in duplicated event and non duplicated events
rngInitial = getRanges(object)
rngInitial$currIdx = seq_along(rngInitial)
# remove possibility of partial overlaps by setting `minoverlap` parameter
# to the size of the binding sites
# -> use the minimum of the size for uneven binding sites sizes
# -> this is not optimal since partial overlaps can be re-introduced
# -> does not affect ranges that are true binding sites
bs.set.range = min(unique(width(rngInitial)))
ols = findOverlaps(anno.genes, rngInitial, minoverlap = bs.set.range)
# ols = findOverlaps(anno.genes, rngInitial)
# Deal with multiple loci
# --------------------------------------------------------------------------
# Count out how many cases there are
countOlsSize = as.data.frame(table(duplicated(subjectHits(ols))))
if (nrow(countOlsSize) == 1){
# no duplications are present
# -> no overlaps have to be resolved
totalBS = countOlsSize$Freq[1]
duplicatedBS = 0
} else {
totalBS = countOlsSize$Freq[1]
duplicatedBS = countOlsSize$Freq[2]
}
duplicatedFraction = paste0(round(duplicatedBS / totalBS * 100, digits = 2), "%")
# ----
# Store results for plotting
# -> select data
if (all(!is.null(selectID) & !is.null(selectType) & !is.null(selectName))) {
# nothing is missing
dfPlot = data.frame(geneIndex = queryHits(ols), bsIndex = subjectHits(ols),
geneType = selectType[queryHits(ols)]) %>% #TODO type
dplyr::mutate(value = 1) %>%
dplyr::select(-geneIndex) %>%
tidyr::pivot_wider(names_from = geneType, values_from = value, values_fn = length, values_fill = 0) %>%
dplyr::select(-bsIndex) %>%
dplyr::mutate_all(~ ifelse(. > 1, 1, .))
} else {
# check what is missing
if (!is.null(selectType)) {
# anything other than gene type is missing
dfPlot = data.frame(geneIndex = queryHits(ols), bsIndex = subjectHits(ols),
geneType = selectType[queryHits(ols)]) %>% #TODO type
dplyr::mutate(value = 1) %>%
dplyr::select(-geneIndex) %>%
tidyr::pivot_wider(names_from = geneType, values_from = value, values_fn = length, values_fill = 0) %>%
dplyr::select(-bsIndex) %>%
dplyr::mutate_all(~ ifelse(. > 1, 1, .))
} else {
# gene type is missing
dfPlot = NULL
}
}
object@plotData$assignToGenes$dataOverlaps = dfPlot
# Case where there are no overlaps
if (duplicatedBS == 0) {
msg = paste0("No binding sites (",
format(totalBS, big.mark = ',', decimal.mark = "."),
") on overlapping gene loci. \n")
msg2 = paste0("Parameter 'overlaps' set to '",
overlaps, "' not in use. \n")
if(!quiet) message(c(msg, msg2))
}
# Cases where there are multiple overlaps
msg0 = paste0(duplicatedFraction, " (", duplicatedBS,"/", totalBS, ")",
" of binding sites overlap with multiple anno.genes in the given gene annotation. \n")
if (overlaps == "hierarchy") {
# apply hierarchy solution
msg1 = paste0("Apply 'hierarchy' solution")
if(!quiet) message(c(msg0, msg1))
ruleMod = data.frame(gene_type = overlaps.rule, idx = seq_along(overlaps.rule))
rngResolved = .resolveGeneOverlapsWithRule(rng = rngInitial, ols = ols, rule = ruleMod,
selectID = selectID, selectName = selectName,
selectType = selectType)
if (length(rngInitial) != length(rngResolved)) {
if(!quiet) message(paste0(format(length(rngInitial)-length(rngResolved), big.mark = ",", decimal.mark = "."),
" binding sites could not be assigned to a gene range.\n"))
}
}
if (overlaps == "frequency") {
# apply frequency solution
msg1 = paste0("Apply 'frequency' solution")
if (!quiet) message(c(msg0, msg1))
ruleFreq = selectType %>% table() %>% as.data.frame() %>%
rename(gene_type = 1) %>% arrange(desc(Freq)) %>%
mutate(idx = row_number()) %>% rename(n = 2)
rngResolved = .resolveGeneOverlapsWithRule(rng = rngInitial, ols = ols, rule = ruleFreq,
selectID = selectID, selectName = selectName, selectType = selectType)
if (length(rngInitial) != length(rngResolved)) {
if (!quiet) message(paste0(format(length(rngInitial)-length(rngResolved), big.mark = ",", decimal.mark = "."),
" binding sites could not be assigned to a gene range.\n"))
}
}
if (overlaps == "remove") {
# apply remove option
msg1 = "Binding sites from overlapping loci are removed. This is not recommended. Please see options 'heirarchy' and 'frequency'. \n"
warning(c(msg0, msg1))
rngResolved = rngInitial[subjectHits(ols)]
# manage meta columns
if (!is.null(selectID)) {
mcols(rngResolved) = cbind(mcols(rngResolved), geneID = selectID[queryHits(ols)])
}
if (!is.null(selectType)) {
mcols(rngResolved) = cbind(mcols(rngResolved), geneType = selectType[queryHits(ols)])
}
if (!is.null(selectName)) {
mcols(rngResolved) = cbind(mcols(rngResolved), geneName = selectName[queryHits(ols)])
}
# remove duplicates
rngResolved = rngResolved[countOverlaps(rngResolved) == 1]
}
if (overlaps == "keep") {
msg1 = "Binding sites from overlapping loci are kept. This is not recommended. Please see options 'heirarchy' and 'frequency'. \n"
warning(c(msg0, msg1))
rngResolved = rngInitial[subjectHits(ols)]
# manage meta columns
if (!is.null(selectID)) {
mcols(rngResolved) = cbind(mcols(rngResolved), geneID = selectID[queryHits(ols)])
}
if (!is.null(selectType)) {
mcols(rngResolved) = cbind(mcols(rngResolved), geneType = selectType[queryHits(ols)])
}
if (!is.null(selectName)) {
mcols(rngResolved) = cbind(mcols(rngResolved), geneName = selectName[queryHits(ols)])
}
}
# ---
# Store results for plotting
if (all(!is.null(selectID) & !is.null(selectType) & !is.null(selectName))) {
# nothing is missing
nCountGenes = rngResolved %>% mcols() %>% as.data.frame() %>%
group_by(geneID) %>% dplyr::count(geneType) %>% ungroup() %>%
dplyr::count(geneType) %>% arrange(desc(n))
nCountBs = rngResolved %>% mcols() %>% as.data.frame() %>%
dplyr::count(geneType) %>% arrange(desc(n))
nCountAll = left_join(nCountBs, nCountGenes, by = "geneType") %>%
rename("nBs" = "n.x", "nGenes" = "n.y")
} else {
# check what is missing
if (!is.null(selectType)) {
# anything other than gene type is missing
nCountGenes = rngResolved %>% mcols() %>% as.data.frame() %>%
group_by(geneID) %>% dplyr::count(geneType) %>% ungroup() %>%
dplyr::count(geneType) %>% arrange(desc(n))
nCountBs = rngResolved %>% mcols() %>% as.data.frame() %>%
dplyr::count(geneType) %>% arrange(desc(n))
nCountAll = left_join(nCountBs, nCountGenes, by = "geneType") %>%
rename("nBs" = "n.x", "nGenes" = "n.y")
}
else {
# gene type is missing
nCountAll = NULL
}
}
object@plotData$assignToGenes$dataSpectrum = nCountAll
# ---
# Store for results
resultLine = data.frame(
funName = "assignToGenes()", class = "binding sites",
nIn = length(rngInitial), nOut = length(rngResolved),
per = paste0(round(length(rngResolved)/ length(rngInitial), digits = 2)*100,"%"),
options = paste0("Overlaps=", optstr$overlaps, ", Source=", optstr$source,
ifelse(!rlang::is_empty(optstr$overlaps.rule),
paste0(", overlaps.rule=",
paste(optstr$overlaps.rule, collapse = ">")), ""))
)
object@results = rbind(object@results, resultLine)
objectOut = setRanges(object, rngResolved, quiet = quiet)
return(objectOut)
}
#' Assign binding sites to their hosting transcript regions
#'
#' Function that assigns each binding site in the \code{\link{BSFDataSet}} to its
#' hosting transcript region given an annotation database (\code{anno.annoDB}), or
#' a GRanges list / \code{\link{CompressedGRangesList}} (\code{anno.transcriptRegionList})
#' that holds the ranges for the transcript regions of interest.
#'
#' Since the assignment is based on the overlaps of annotated transcript ranges
#' and binding sites, no matching meta data is needed.
#'
#' In the case of transcript regions overlaps are very frequent. To resolve these
#' cases the \code{\link{overlaps}} argument can be used. Option `frequency`
#' will take the most frequently observed transcript region, option `hierarchy`
#' works in conjunction with a user defined rule (\code{overlaps.rule}).
#' Options `flag` and `remove` will label binding sites with an ambiguous tag or
#' remove all overlapping cases, respectively.
#'
#' The function is part of the standard workflow performed by \code{\link{BSFind}}.
#'
#' @param object a \code{\link{BSFDataSet}} object with stored binding sites. This
#' means that ranges should be > 1
#' @param overlaps character; how overlapping transcript regions should be handled.
#' @param overlaps.rule character vector; a vector of transcript region
#' names that should be used to handle overlapping cases in a hierarchical manor.
#' The order of the vector is the order of the hierarchy.
#' @param anno.annoDB an object of class \code{OrganismDbi} that contains
#' the transcript region annotation (!!! Experimental !!!).
#' @param anno.transcriptRegionList an object of class \code{\link{CompressedGRangesList}}
#' that holds an ranges for each transcript region
#'
#' @param normalize.exclude.upper numeric; percentage value that indicates the
#' upper boundary for transcript region length to be considered when calculating
#' normalization factors for regions.
#' @param normalize.exclude.lower numeric; percentage value that indicates the
#' lower boundary for transcript region length to be considered when calculating
#' normalization factors for regions.
#'
#' @param quiet logical; whether to print messages
#'
#' @return an object of class \code{\link{BSFDataSet}} with binding sites having
#' hosting transcript region information added to their meta columns.
#'
#' @seealso \code{\link{BSFind}},
#' \code{\link{transcriptRegionOverlapsPlot}},
#' \code{\link{transcriptRegionSpectrumPlot}}
#'
#' @importFrom dplyr rename_with desc
#' @importFrom GenomicFeatures cds intronsByTranscript threeUTRsByTranscript fiveUTRsByTranscript
#'
#' @examples
#' # load clip data
#' files <- system.file("extdata", package="BindingSiteFinder")
#' load(list.files(files, pattern = ".rda$", full.names = TRUE))
#' load(list.files(files, pattern = ".rds$", full.names = TRUE)[1])
#' load(list.files(files, pattern = ".rds$", full.names = TRUE)[2])
#'
#' bds = makeBindingSites(object = bds, bsSize = 9)
#' bds = assignToGenes(bds, anno.genes = gns)
#' bds = assignToTranscriptRegions(object = bds, anno.transcriptRegionList = regions)
#'
#' @export
assignToTranscriptRegions <- function(object, # bindingSiteFinder
overlaps = c("frequency", "hierarchy", "flag", "remove"), # overlapping loci solution options
overlaps.rule = NULL, # rule to apply when option hierarchy is selected
anno.annoDB = NULL, # annotation data base of class OrganismDbi; can be NULL -> but then anno.transcriptRegionList must be provided
anno.transcriptRegionList = NULL, # GRangesList object with the transcript features to be used; can be NULL -> but then anno.annoDB must be provided
normalize.exclude.upper = 0.02,
normalize.exclude.lower = 0.02,
quiet = FALSE
) {
# local variables
datasource <- `.` <- Freq <- transcript_type <- NULL
# INPUT CHECKS
# --------------------------------------------------------------------------
stopifnot(is(object, "BSFDataSet"))
stopifnot(is.logical(quiet))
# Check if none is specified
if (is.null(anno.annoDB) & is.null(anno.transcriptRegionList)) {
msg = paste0("None of the required annotation sources anno.annoDB or anno.transcriptRegionList was specified. ")
stop(msg)
}
# Check if both are specified
if (!is.null(anno.annoDB) & !is.null(anno.transcriptRegionList)) {
msg = paste0("Both of the required annotation sources anno.annoDB or anno.transcriptRegionList are specified. Please provide only one of the two. ")
stop(msg)
}
# anno.annoDB should be used
if (!is.null(anno.annoDB) & is.null(anno.transcriptRegionList)) {
stopifnot(is(anno.annoDB, "OrganismDb"))
if (!is.null(anno.transcriptRegionList)) {
msg = paste0("Parameter anno.annoDB and anno.transcriptRegionList are specified at the same time. Use only one of them.")
stop(msg)
} else {
datasource = "anno.annoDB"
# extract relevant annotation - transcripts
cdseq = cds(anno.annoDB)
intrns = unlist(intronsByTranscript(anno.annoDB))
utrs3 = unlist(threeUTRsByTranscript(anno.annoDB))
utrs5 = unlist(fiveUTRsByTranscript(anno.annoDB))
anno.transcriptRegionList = GRangesList(CDS = cdseq, Intron = intrns, UTR3 = utrs3, UTR5 = utrs5)
names(anno.transcriptRegionList) = toupper(names(anno.transcriptRegionList))
}
}
# anno.transcriptRegionList should be used
if (is.null(anno.annoDB) & !is.null(anno.transcriptRegionList) ) {
stopifnot(is(anno.transcriptRegionList, "GenomicRangesList"))
if (!is.null(anno.annoDB)) {
msg = paste0("Parameter anno.annoDB and anno.transcriptRegionList are specified at the same time. Use only one of them.")
stop(msg)
} else {
datasource = "anno.transcriptRegionList"
# extract relevant annotation - transcripts
anno.transcriptRegionList = anno.transcriptRegionList
names(anno.transcriptRegionList) = toupper(names(anno.transcriptRegionList))
}
}
# handle options (hierarchy, frequency, remove, keep)
overlaps = match.arg(overlaps, choices = c("frequency", "hierarchy", "flag", "remove"))
# Check multiple loci options
if (overlaps == "hierarchy" & is.null(overlaps.rule)) {
msg1 = paste0("Multiple region problem is set to be handled by 'hierarchy', but no overlaps.rule is provided. \n")
msg2 = paste0("Change the 'overlaps' option or provide a valid overlaps.rule. \n")
stop(c(msg1, msg2))
}
overlaps.rule = toupper(overlaps.rule)
if (! all(overlaps.rule %in% names(anno.transcriptRegionList))){
msg = paste0("Regions defined in overlaps.rule (", paste(overlaps.rule, collapse = " > "),
") does not match the input region names (", paste(names(anno.transcriptRegionList), collapse = ", "), ")")
stop(msg)
}
if (overlaps == "hierarchy" & is.null(overlaps.rule)) {
msg1 = paste0("Overlap are set to be solved by 'hierarchy', but no overlaps.rule is provided. \n")
msg2 = paste0("Change the 'overlaps' option or provide a valid overlaps.rule. \n")
stop(c(msg1, msg2))
}
# ---
# Store function parameters in list
optstr = list(source = datasource, overlaps = overlaps, overlaps.rule = overlaps.rule)
object@params$assignToTranscriptRegions = optstr
# MAIN COMPUTE
# --------------------------------------------------------------------------
# Match transcript annotation with binding sites
# --------------------------------------------------------------------------
rngInitial = getRanges(object)
rng = getRanges(object)
cRange = resize(rng, width = 1, fix = "center")
countDf = lapply(anno.transcriptRegionList, function(x) {
countOverlaps(cRange, x)
}) %>% as.data.frame() %>% rename_with(toupper)
# Add intergenic cases
countDf = countDf %>% mutate(OTHER = ifelse(rowSums(.) == 0, 1, 0))
# ---
# Store results for plotting
plotDf = countDf
plotDf[plotDf >= 1] = 1
object@plotData$assignToTranscriptRegions$dataOverlaps = plotDf
# Deal with multiple loci
# --------------------------------------------------------------------------
# Count out how many cases there are
totalBs = nrow(countDf)
numberOfLociHit = rowSums(countDf != 0)
ambigousBs = nrow(countDf[numberOfLociHit != 1,])
ambigousBsFraction = paste0(round(ambigousBs / totalBs * 100, digits = 2), "%")
# Case where there are no overlaps
if (ambigousBs == 0) {
msg0 = paste0("No binding sites (", format(totalBs, big.mark = ',', decimal.mark = "."), ") on overlapping transcript loci. \n")
msg1 = paste0("Parameter 'overlaps' set to '", overlaps, "' not in use. \n")
if(!quiet) message(c(msg0, msg1))
# assign regions
countDfSorted = countDf
countDfSorted[countDfSorted > 0] = 1
col_names = colnames(countDfSorted)
transcriptRegion = ifelse(rowSums(countDfSorted) == 1, col_names[max.col(countDfSorted, ties.method = "first")], "Ambiguous")
mcols(rng) = cbind.data.frame(mcols(rng), transcriptRegion = transcriptRegion)
}
# Cases where we have to resolve overlapping loci
if (ambigousBs > 0) {
msg0 = paste0(ambigousBsFraction, " (", ambigousBs,"/", totalBs, ")",
" of binding sites overlap with multiple different transcript regions in the given annotation. \n")
if (overlaps == "hierarchy") {
# apply hierarchy solution
msg1 = paste0("Apply 'hierarchy' solution. \n ")
if(!quiet) message(c(msg0, msg1))
countDfSorted = countDf[, overlaps.rule]
countDfSorted[countDfSorted > 0] = 1
transcriptRegion = names(countDfSorted)[max.col(countDfSorted, ties.method = "first")]
mcols(rng) = cbind.data.frame(mcols(rng), transcriptRegion = transcriptRegion)
}
if (overlaps == "frequency") {
msg1 = paste0("Apply 'frequency' solution. \n ")
if(!quiet) message(c(msg0, msg1))
countDfSorted = countDf
# countDfSorted[countDfSorted > 0] = 1
# get overlaps.rule by frequency
overlaps.rule = names(countDfSorted)[max.col(countDfSorted, ties.method = "first")] %>%
table() %>% as.data.frame() %>% rename(transcript_type = 1) %>%
arrange(desc(Freq)) %>% pull(transcript_type) %>% as.character()
# apply frequency solution
countDfSorted = countDfSorted[, overlaps.rule]
transcriptRegion = names(countDfSorted)[max.col(countDfSorted, ties.method = "first")]
mcols(rng) = cbind.data.frame(mcols(rng), transcriptRegion = transcriptRegion)
}
if (overlaps == "flag") {
msg1 = paste0("Binding sites with multiple different overlapping transcript regions are being flaged with the tag 'Ambiguous'. \n")
if(!quiet) message(c(msg0, msg1))
countDfSorted = countDf
countDfSorted[countDfSorted > 0] = 1
# flag multiple loci binding sites
col_names = colnames(countDfSorted)
transcriptRegion = ifelse(rowSums(countDfSorted) == 1, col_names[max.col(countDfSorted, ties.method = "first")], "Ambiguous")
mcols(rng) = cbind.data.frame(mcols(rng), transcriptRegion = transcriptRegion)
# remove mulitple loci binding sitges
mcols(rng) = cbind.data.frame(mcols(rng), transcriptRegion = transcriptRegion)
}
if (overlaps == "remove") {
msg1 = paste0("Binding sites with multiple different overlapping transcript regions are being removed. This is not recommended! \n")
if(!quiet) message(c(msg0, msg1))
countDfSorted = countDf
countDfSorted[countDfSorted > 0] = 1
# flag multiple loci binding sites
col_names = colnames(countDfSorted)
transcriptRegion = ifelse(rowSums(countDfSorted) == 1, col_names[max.col(countDfSorted, ties.method = "first")], "Ambiguous")
mcols(rng) = cbind.data.frame(mcols(rng), transcriptRegion = transcriptRegion)
# remove mulitple loci binding sitges
mcols(rng) = cbind.data.frame(mcols(rng), transcriptRegion = transcriptRegion)
rng = rng[rng$transcriptRegion != "Ambiguous"]
}
}
# ---
# Calculate width for length based normalization
## normalize by all expressed regions
debugList = list(object = object, rng = rng)
# saveRDS(debugList, file = "./debugList.rds")
norm.factors = .calcNormalizeFactors(object, rng, anno.transcriptRegionList,
normalize.exclude.lower = normalize.exclude.lower,
normalize.exclude.upper = normalize.exclude.upper)
# ---
# Store results for plotting
dfPlot = rng$transcriptRegion %>%
table() %>%
as.data.frame() %>%
rename("TranscriptRegion" = ".") %>%
arrange(desc(Freq)) %>%
left_join(., y = norm.factors, by = "TranscriptRegion")
object@plotData$assignToTranscriptRegions$dataSpectrum = dfPlot
# ---
# Store for results
resultLine = data.frame(
funName = "assignToTranscriptRegions()", class = "binding sites",
nIn = length(rngInitial), nOut = length(rng),
per = paste0(round(length(rng)/ length(rngInitial), digits = 2)*100,"%"),
options = paste0("Overlaps=", optstr$overlaps, ", Source=", datasource,
ifelse(!rlang::is_empty(optstr$overlaps.rule), paste0(", overlaps.rule=", paste(optstr$overlaps.rule, collapse = ">")), ""))
)
object@results = rbind(object@results, resultLine)
outObject = setRanges(object, rng, quiet = quiet)
return(outObject)
}
#' Function to estimate the appropriate binding site width together with the
#' optimal gene-wise filter level.
#'
#' This function tests different width of binding sites for different gene-wise
#' filtering steps. For each test the signal-to-score ratio is calculated. The
#' mean over all gene-wise filterings at each binding site width is used to
#' extract the optimal width, which serves as anchor to select the optimal
#' gene-wise filter.
#'
#' Parameter estimation is done on a subset of all crosslink sites
#' (\code{est.subsetChromosome}).
#'
#' Gene-level filter can be tested with varying
#' levels of accuracy ranging from `finest` to `coarse`, representing 1% and
#' 20% steps, respectively.
#'
#' Binding site computation at each step can be done on three different accuracy
#' level (\code{bsResolution}). Option `fine` is equal to a normal run
#' of the \code{\link{makeBindingSites}} function. `medium` will perform
#' a shorter version of the binding site computation, skipping some of the
#' refinement steps. Option `coarse` will approximate binding sites by merged
#' crosslinks regions, aligning the center at the site with the highest score.
#'
#' For each binding site in each set given the defined resolutions a signal-to-
#' flank score ratio is calculated and the mean of this score per set is returned.
#' Next a mean of means is created which results in a single score for each
#' binding site width that was tested. The width that yielded the highest score
#' is selected as optimal. In addtion the \code{minimumStepGain} option
#' allows control over the minimum additional gain in the score that a tested
#' width has to have to be selected as the best option.
#'
#' To enhance the sensitivity of the binding site estimation, the sensitivity
#' mode exists. In this mode crosslink sites undergo a pre-filtering and merging
#' step, to exclude potential artifical peaks (experimental-, mapping-biases).
#' If sensitivity mode is activated the \code{est.minWidth} option should be set
#' to 1.
#'
#' The optimal geneFilter is selected as the first one that passes the merged
#' mean of the selected optimal binding site width.
#'
#' The function is part of the standard workflow performed by \code{\link{BSFind}}.
#'
#' @param object a \code{\link{BSFDataSet}} object with stored crosslink sites.
#' This means that ranges should have a width = 1.
#' @param bsResolution character; level of resolution at which different binding
#' site width should be tested
#' @param geneResolution character; level of resolution at which gene-wise filtering
#' steps should be tested
#' @param est.maxBsWidth numeric; the largest binding site width which should
#' considered in the testing
#' @param est.minimumStepGain numeric; the minimum additional gain in the score
#' in percent the next binding site width has to have, to be selected as best option
#' @param est.maxSites numeric; maximum number of PureCLIP sites that are used;
#' @param est.subsetChromosome character; define on which chromosome the
#' estimation should be done in function \code{\link{estimateBsWidth}}
#' @param est.minWidth the minimum size of regions that are subjected to the
#' iterative merging routine, after the initial region concatenation.
#' @param est.offset constant added to the flanking count in the signal-to-flank
#' ratio calculation to avoid division by Zero
#'
#' @param sensitive logical; whether to enable sensitive pre-filtering before
#' binding site merging or not
#' @param sensitive.size numeric; the size (in nucleotides) of the merged
#' sensitive region
#' @param sensitive.minWidth numeric; the minimum size (in nucleoties) of the
#' merged sensitive region
#'
#' @param anno.annoDB an object of class \code{OrganismDbi} that contains
#' the gene annotation (!!! Experimental !!!).
#' @param anno.genes an object of class \code{\link{GenomicRanges}} that represents
#' the gene ranges directly
#' @param bsResolution.steps numeric vector; option to use a user defined threshold
#' for binding site width directly. Overwrites \code{bsResolution}
#' @param geneResolution.steps numeric vector; option to use a user defined threshold
#' vector for gene-wise filtering resolution. Overwrites \code{geneResolution}
#' @param quiet logical; whether to print messages
#' @param veryQuiet logical; whether to suppress all messages
#' @param reportScoresPerBindingSite report the ratio score for each binding site
#' separately. Warning! This is for debugging and testing only. Downstream
#' functions can be impaired.
#' @param ... additional arguments passed to \code{\link{pureClipGeneWiseFilter}}
#'
#' @return an object of class \code{\link{BSFDataSet}} with binding sites with
#' the `params` slots `bsSize` and `geneFilter` being filled
#'
#' @seealso \code{\link{BSFind}},
#' \code{\link{estimateBsWidthPlot}}
#'
#' @importFrom dplyr summarise slice_head slice_tail
#' @importFrom utils head
#' @importFrom GenomicFeatures genes
#' @importFrom stats sd
#'
#' @examples
#' # load clip data
#' files <- system.file("extdata", package="BindingSiteFinder")
#' load(list.files(files, pattern = ".rda$", full.names = TRUE))
#' load(list.files(files, pattern = ".rds$", full.names = TRUE)[1])
#' load(list.files(files, pattern = ".rds$", full.names = TRUE)[2])
#' estimateBsWidth(bds, anno.genes = gns, est.maxBsWidth = 19,
#' geneResolution = "coarse", bsResolution = "coarse", est.subsetChromosome = "chr22")
#'
#' @export
estimateBsWidth <- function(object, # BindingSiteFinder object
bsResolution = c("medium", "fine", "coarse"),
geneResolution = c("medium", "coarse", "fine", "finest"),
est.maxBsWidth = 13,
est.minimumStepGain = 0.02,
est.maxSites = Inf,
est.subsetChromosome = "chr1", #
est.minWidth = 2,
est.offset = 1,
sensitive = FALSE,
sensitive.size = 5,
sensitive.minWidth = 2,
anno.annoDB = NULL,
anno.genes = NULL,
bsResolution.steps = NULL,
geneResolution.steps = NULL,
quiet = TRUE,
veryQuiet = FALSE,
reportScoresPerBindingSite = FALSE,
...
) {
# initialize local variables
bsSize <- ms <- growth_per <- sel <- geneWiseFilter <- est.option <- signalToFlankRatio <- NULL
# INPUT CHECKS
# --------------------------------------------------------------------------
stopifnot(is(object, "BSFDataSet"))
stopifnot(is.logical(quiet))
if (c(est.maxBsWidth %% 2) == 0) {
stop("est.maxBsWidth is even. An odd number is required to have a distinct binding site center.")
}
# Check if none is specified
if (is.null(anno.annoDB) & is.null(anno.genes)) {
msg = paste0("None of the required annotation sources anno.annoDB or anno.genes was specified. ")
stop(msg)
}
# Check if both are specified
if (!is.null(anno.annoDB) & !is.null(anno.genes)) {
msg = paste0("Both of the required annotation sources anno.annoDB or anno.genes are specified. Please provide only one of the two. ")
stop(msg)
}
# Checks if anno.annoDB should be used
if (!is.null(anno.annoDB) & is.null(anno.genes)) {
stopifnot(is(anno.annoDB, "OrganismDb"))
if (!is.null(anno.genes)) {
msg = paste0("Parameter anno.annoDB and anno.genes are specified at the same time. Use only one of them.")
stop(msg)
} else {
# extract relevant annotation
anno.genes = genes(anno.annoDB)
}
}
# Checks if anno.genes should be used
if (is.null(anno.annoDB) & !is.null(anno.genes)) {
stopifnot(is(anno.genes, "GenomicRanges"))
if (!is.null(anno.annoDB)) {
msg = paste0("Parameter anno.annoDB and anno.genes are specified at the same time. Use only one of them.")
stop(msg)
} else {
# extract relevant annotation
anno.genes = anno.genes
}
}
# handle gene filter resolution
geneResolution = match.arg(geneResolution, choices = c("medium", "coarse", "fine", "finest"))
if (!is.null(geneResolution.steps)) {
geneResolution.steps = geneResolution.steps
} else {
if (geneResolution == "coarse") {
geneResolution.steps = seq(from = 0, to = 1, by = 0.2)
geneResolution.steps = geneResolution.steps[1:length(geneResolution.steps)-1]
}
if (geneResolution == "medium") {
geneResolution.steps = seq(from = 0, to = 1, by = 0.1)
geneResolution.steps = geneResolution.steps[1:length(geneResolution.steps)-1]
}
if (geneResolution == "fine") {
geneResolution.steps = seq(from = 0, to = 1, by = 0.05)
geneResolution.steps = geneResolution.steps[1:length(geneResolution.steps)-1]
}
if (geneResolution == "finest") {
geneResolution.steps = seq(from = 0, to = 1, by = 0.01)
geneResolution.steps = geneResolution.steps[1:length(geneResolution.steps)-1]
}
}
# handle binding site compute resolution
bsResolution = match.arg(bsResolution, choices = c("medium", "fine", "coarse"))
# handle est.maxBsWidth
if (!is.null(bsResolution.steps)) {
bsResolution.steps = bsResolution.steps
} else {
bsResolution.steps = seq(from = 3, to = est.maxBsWidth, by = 2)
}
# ---
# Store function parameters in list
optstr = list(bsResolution = bsResolution,
geneResolution = geneResolution,
est.maxBsWidth = est.maxBsWidth,
est.maxSites = est.maxSites,
est.minimumStepGain = est.minimumStepGain,
est.subsetChromosome = est.subsetChromosome,
est.minWidth = est.minWidth,
est.offset = est.offset,
sensitive = sensitive,
sensitive.size = sensitive.size,
sensitive.minWidth = sensitive.minWidth)
object@params$estimateBsWidth = optstr
# PREPARE TEST RANGES + SIGNAL
# --------------------------------------------------------------------------
checkRng = getRanges(object)
# check if subset is part of the seqnames from ranges
if (!all(est.subsetChromosome %in% levels(seqnames(checkRng)))){
msg = paste0("Chromosome to estimate on (", est.subsetChromosome, "), is not included in the ranges: ", paste(levels(seqnames(checkRng)), collapse = ", "))
stop(msg)
}
# limit estimation to a specific chromosome
if (!is.null(est.subsetChromosome)) {
redObj = .subsetByChr(object, chr = est.subsetChromosome, quiet = quiet)
} else {
redObj = object
}
# limit estimation to a maximum number of sites
estRng = getRanges(redObj)
if (length(estRng) > est.maxSites) {
estRng = head(estRng, est.maxSites)
redObj = setRanges(redObj, estRng, quiet = quiet)
}
# limit the clip signal to the ranges used for estimation (plus some extra offset)
if (!is.null(est.subsetChromosome)) {
# reduce signal to frame
maxFrame = ceiling((max(bsResolution.steps) *3))
redObj = .reduceSignalToFrame(redObj, frame = maxFrame, quiet = quiet)
} else {
# don't reduce signal
redObj = redObj
}
# collapse signal from replicates
sgnMerge = .collapseSamples(getSignal(redObj))
# MAIN COMPUTE
# --------------------------------------------------------------------------
# a counter that shows to how many percent the computation is done
counterTotalIterations = length(geneResolution.steps) * length(bsResolution.steps)
counterCurrentIterations = 0
# print start message
msg = paste0("Estimation at: ", round((counterCurrentIterations/counterTotalIterations)*100 ), "%" )
if(!veryQuiet) message(msg)
# calculate binding sites for each filter step and width
scoreAllDf = lapply(geneResolution.steps, function(bsFilterStep){
# bsFilterStep = 0.2
# apply current gene-wise filter (has two modes)
if (isTRUE(sensitive)) {
# -> sensitive mode
# define regions in which sites should be kept
cRng = getRanges(redObj)
cRngMerge = reduce(cRng, min.gapwidth = sensitive.size)
cRngMerge = cRngMerge[width(cRngMerge) >= sensitive.minWidth]
# currFilterObj = pureClipGeneWiseFilter(object = redObj, anno.genes = anno.genes, cutoff = bsFilterStep, quiet = quiet, ...)
currFilterObj = pureClipGeneWiseFilter(object = redObj, anno.genes = anno.genes, cutoff = bsFilterStep, quiet = quiet)
cRngFilter = getRanges(currFilterObj)
currRng = IRanges::subsetByOverlaps(cRngFilter, cRngMerge)
} else {
# -> normal mode
# currFilterObj = pureClipGeneWiseFilter(object = redObj, anno.genes = anno.genes, cutoff = bsFilterStep, quiet = quiet, ...)
currFilterObj = pureClipGeneWiseFilter(object = redObj, anno.genes = anno.genes, cutoff = bsFilterStep, quiet = quiet)
currRng = getRanges(currFilterObj)
}
# calculate binding sites for current bsWidth
rngPerWidth = lapply(bsResolution.steps, function(bsWidthStep){
if (bsResolution == "fine") {
# calculate with full binding sites
currBsObj = makeBindingSites(object = currFilterObj,
bsSize = bsWidthStep,
minWidth = est.minWidth,
quiet = quiet)
currBs = getRanges(currBsObj)
}
if (bsResolution == "medium") {
# approximate binding sites by a single merge and extend round
currBs = .approximateBindingSites_medium(rng = currRng,
sgn = sgnMerge,
bsSize = bsWidthStep,
minWidth = est.minWidth)
}
if (bsResolution == "coarse") {
# approximate binding sites by reduced pureclip sites center
currBs = .approximateBindingSites_coarse(rng = currRng,
bsSize = bsWidthStep,
minWidth = est.minWidth)
}
return(currBs)
})
rngPerWidth = unlist(GRangesList(rngPerWidth))
# handle all plus ranges
cRangePlus = subset(rngPerWidth, strand == "+")
if (length(cRangePlus) > 0) {
bsSumPlus = sum(sgnMerge$signalPlus[cRangePlus])
extendedRangePlus = cRangePlus + cRangePlus$bsSize
exSumPlus = sum(sgnMerge$signalPlus[extendedRangePlus])
# extendedRangePlus$signalToFlankRatio = (bsSumPlus) / (((exSumPlus - bsSumPlus) / 2) + est.offset)
extendedRangePlus$signalToFlankRatio = (bsSumPlus / exSumPlus)
} else {
extendedRangePlus = cRangePlus
}
# handle all minus ranges
cRangeMinus = subset(rngPerWidth, strand == "-")
if (length(cRangeMinus) > 0) {
bsSumMinus = sum(sgnMerge$signalMinus[cRangeMinus])
extendedRangeMinus = cRangeMinus + cRangeMinus$bsSize
exSumMinus = sum(sgnMerge$signalMinus[extendedRangeMinus])
# extendedRangeMinus$signalToFlankRatio = (bsSumMinus) / (((exSumMinus - bsSumMinus) / 2) + est.offset)
extendedRangeMinus$signalToFlankRatio = (bsSumMinus / exSumMinus)
} else {
extendedRangeMinus = cRangeMinus
}
if (isTRUE(reportScoresPerBindingSite)) {
# report the scores for each binding site
# -> does not allow the usage of the standard plot
df = rbind(as.data.frame(mcols(extendedRangeMinus)), as.data.frame(mcols(extendedRangePlus))) %>%
group_by(bsSize) %>%
mutate(geneWiseFilter = bsFilterStep)
} else {
# combine both and make results dataframe
df = rbind(as.data.frame(mcols(extendedRangeMinus)), as.data.frame(mcols(extendedRangePlus))) %>%
group_by(bsSize) %>%
summarize(signalToFlankRatio = median(signalToFlankRatio), .groups = "keep") %>%
mutate(geneWiseFilter = bsFilterStep)
}
# update chunk counter
counterCurrentIterations <<- counterCurrentIterations + length(bsResolution.steps)
msg = paste0("Estimation at: ", round((counterCurrentIterations/counterTotalIterations)*100 ), "%" )
if(!veryQuiet) message(msg)
return(df)
})
df = do.call("rbind", scoreAllDf)
# ---
# Store results for plotting
if (isTRUE(reportScoresPerBindingSite)){
object@plotData$estimateBsWidth$dataEnhanced = df
} else {
object@plotData$estimateBsWidth$data = df
}
resEstimate = .findFirstMaximum(df, est.minimumStepGain = est.minimumStepGain)
est.bsSize = resEstimate$est.bsSize
est.geneFilter = resEstimate$est.geneFilter
est.option = resEstimate$est.option
if (est.option == "local") {
msg = paste0("No global maximum found with est.maxBsWidth=", est.maxBsWidth, ". Using local maximum instead.\n")
if(!veryQuiet) warning(msg)
}
if (est.option == "fallback") {
msg = paste0("No global or local maximum found with est.maxBsWidth=", est.maxBsWidth, ". Falling back to simple maximum instead.\n")
if(!veryQuiet) warning(msg)
}
if (est.option == "error") {
msg = paste0("No maximum found at all. Raise 'est.minimumStepGain' argument/, change bsResolution and/or geneResolution arguments/ raise est.maxBsWidth and check visually. \n")
stop(msg)
}
# Store estimate type in options (global/ local)
optstr$option = est.option
object@params$estimateBsWidth = optstr
# Store estimated parameters
object@params$bsSize = est.bsSize
object@params$geneFilter = est.geneFilter
# ---
# Store for results
rng = getRanges(object)
resultLine = data.frame(
funName = "estimateBsWidth()", class = "estimate",
nIn = length(rng), nOut = length(rng),
per = paste0(round(length(rng)/ length(rng), digits = 2)*100,"%"),
options = paste0("bsResolution=", bsResolution, ", geneResolution=", geneResolution,
", est.maxBsWidth=", est.maxBsWidth, ", est.maxSites=", est.maxSites, ", est.subsetChromosome=", paste(est.subsetChromosome, collapse = ","))
)
object@results = rbind(object@results, resultLine)
return(object)
}
ZarnackGroup/BindingSiteFinder documentation built on May 31, 2024, 3:29 a.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 assignToGenes pureClipGeneWiseFilter pureClipGlobalFilter
Documented in assignToGenes pureClipGeneWiseFilter pureClipGlobalFilter
################################################################################
# Workflow functions
# - these functions are used in the main bindingSiteFinde workflow
# - each function updates the bindingSiteFinder object
# - functions are either designed to work on single nt wide crosslink sites or
# on merged binding sites
# - for each function a related plotting funciton exists to visualize the results
# - the wrapper function BSFind() executes the functions below in the desired order
################################################################################
#' Filter PureCLIP sites by their score distribution
#'
#' Function that applies a filter on the global crosslink site score distribution.
#' The \code{\link{GenomicRanges}} contained in the \code{\link{BSFDataSet}} need to
#' have a meta-column that holds a numeric score value, which is used for filtering.
#' The name of the column can be set with \code{match.score}.
#'
#' The function is part of the standard workflow performed by \code{\link{BSFind}}.
#'
#' @param object a \code{\link{BSFDataSet}} object with stored crosslink ranges
#' of width=1
#' @param cutoff numeric; defines the cutoff for which sites to keep, the
#' smallest step is 1\% (0.01)
#' @param match.score character; meta column name of the crosslink site
#' \code{\link{GenomicRanges}} object that holds the score which is used for
#' sub-setting
#' @param quiet logical; whether to print messages
#'
#' @return an object of class \code{\link{BSFDataSet}} with its ranges filtered
#' by those that passed the threshold set with \code{cutoff}
#'
#' @seealso \code{\link{BSFind}}, \code{\link{pureClipGlobalFilterPlot}}
#'
#' @examples
#' # load data
#' files <- system.file("extdata", package="BindingSiteFinder")
#' load(list.files(files, pattern = ".rda$", full.names = TRUE))
#' # apply 5% filter
#' pureClipGlobalFilter(object = bds, cutoff = 0.05)
#'
#' @export
pureClipGlobalFilter <- function(object, # bindingSiteFinder
cutoff = 0.01, # defines the cutoff for which PureCLIP sites to keep; smallest are 1% (0.01); a cutoff of 0.05 means to remove the lowest 5% PureCLIP sites base on the score
match.score = "score",
quiet = FALSE
) {
# INPUT CHECKS
# --------------------------------------------------------------------------
stopifnot(is(object, "BSFDataSet"))
stopifnot(is.logical(quiet))
# check for correct cutoff value
possibleCutoffs = seq(from = 0, to = 1, by = 0.01)
if (! cutoff %in% possibleCutoffs) {
msg = paste0("Specified cutoff (", cutoff,
"), does not refer to a percentage value. Try cutoffs of the form 0.01, 0.05, ... instead.")
stop(msg)
}
# Check for correct score column
metaNames = colnames(mcols(getRanges(object)))
if (! match.score %in% metaNames) {
msg = paste0("Specified match.score (", match.score, "), is present in the ranges. Add respective column. ")
stop(msg)
}
# Create matching vectors for score column from input source
# --------------------------------------------------------------------------
rngInitial = getRanges(object)
score = mcols(rngInitial)[match(match.score, colnames(mcols(rngInitial)))][[1]]
# ---
# Store function parameters in list
optstr = list(cutoff = cutoff, match.score = match.score)
object@params$pureClipGlobalFilter = optstr
# MAIN COMPUTE
# --------------------------------------------------------------------------
# compute quantile based on PureCLIP score
qCut = quantile(score, probs = seq(0,1, by = 0.01))
# apply user set cutoff
names(qCut) = as.numeric(sub("%", "", names(qCut))) / 100
idx = match(cutoff, names(qCut))
rngFiltered = rngInitial[score >= qCut[idx]]
# ---
# Store for plotting
dfPlot = score
cutoffPoint = qCut[idx]
object@plotData$pureClipGlobalFilter$data = dfPlot
object@plotData$pureClipGlobalFilter$cutoffPoint = cutoffPoint
# ---
# Store for results
resultLine = data.frame(
funName = "pureClipGlobalFilter()", class = "crosslink sites",
nIn = length(rngInitial), nOut = length(rngFiltered),
per = paste0(round(length(rngFiltered)/ length(rngInitial), digits = 2)*100,"%"),
options = paste0("Cutoff=", optstr$cutoff, ", scoreMatchCol=", optstr$match.score)
)
object@results = rbind(object@results, resultLine)
# return BSF object with only ranges above cutoff
objectNew = setRanges(object, rngFiltered, quiet = quiet)
return(objectNew)
}
#' Filter PureCLIP sites by their score distribution per gene
#'
#' Function that applies a filter on the crosslink site score distribution at
#' gene level. This allows to filter for those sites with the strongest signal
#' on each gene. Since scores are tied to the expression level of the hosting
#' transcript this function allows a fair filter for all genes partially
#' independent of the expression level.
#'
#' The \code{\link{GenomicRanges}} contained in the \code{\link{BSFDataSet}} need to
#' have a meta-column that holds a numeric score value, which is used for filtering.
#' The name of the column can be set with \code{scoreCol}.
#'
#' In the case of overlapping gene annotation, a single crosslink site will be
#' attributed to multiple genes. The \code{\link{overlaps}} parameter allows
#' to control these cases. Option `keepSingle` will only keep a single instance
#' of the site; `removeAll` will remove both sites; `keepAll` will keep both
#' sites.
#'
#' The function is part of the standard workflow performed by \code{\link{BSFind}}.
#'
#' @param object a \code{\link{BSFDataSet}} object with stored crosslink ranges
#' of width=1
#' @param cutoff numeric; defines the cutoff for which sites to remove, the
#' smallest step is 1\% (0.01). A cutoff of 5\% will remove the lowest 5\% sites,
#' given their score, on each gene, thus keeping the strongest 95\%.
#' @param overlaps character; how overlapping gene loci should be handled.
#' @param anno.annoDB an object of class \code{OrganismDbi} that contains
#' the gene annotation (!!! Experimental !!!).
#' @param anno.genes an object of class \code{\link{GenomicRanges}} that represents
#' the gene ranges directly
#' @param match.score character; meta column name of the crosslink site
#' \code{\link{GenomicRanges}} object that holds the score which is used for
#' sub-setting
#' @param match.geneID character; meta column name of the genes
#' \code{\link{GenomicRanges}} object that holds a unique geneID
#' @param quiet logical; whether to print messages
#'
#' @return an object of class \code{\link{BSFDataSet}} with its ranges filtered
#' by those that passed the gene-wise threshold set with \code{cutoff}
#'
#' @seealso \code{\link{BSFind}}, \code{\link{estimateBsWidthPlot}}
#'
#' @importFrom dplyr filter reframe group_by mutate mutate_all arrange row_number
#' @importFrom GenomicFeatures genes
#' @importFrom stats density
#' @importFrom grDevices colorRampPalette
#'
#' @examples
#' # load clip data
#' files <- system.file("extdata", package="BindingSiteFinder")
#' load(list.files(files, pattern = ".rda$", full.names = TRUE))
#' # Load GRanges with genes
#' load(list.files(files, pattern = ".rds$", full.names = TRUE)[1])
#' # apply 5% gene-wise filter
#' pureClipGeneWiseFilter(object = bds, anno.genes = gns, cutoff = 0.5, overlaps = "keepSingle")
#'
#' @export
pureClipGeneWiseFilter <- function(object, # bindingSiteFinder
cutoff = 0.05, # defines the cutoff for which PureCLIP sites to keep; smallest steps are in 1% (0.01); cutoff of 0.7 means that the top 70% of PureCLIP sites per gene are retained.
overlaps = c("keepSingle", "removeAll", "keepAll"), # Options to deal with PureCLIP sites on overlapping loci genes; use keepSingle as default; keepAll will inflate the number of binding sites by the number of overlaps, removeAll will kill all duplicated instances
anno.annoDB = NULL, # annotation data base of class OrganismDbi; can be NULL -> but then anno.genes must be provided
anno.genes = NULL, # ranges of the genes of class GenomicRanges; can be NULL -> but then anno.annoDB must be provided
match.score = "score",
match.geneID = "gene_id",
quiet = FALSE) {
# Local variables
datasource <- cutoffs <- quantiles <- name <- NULL
# INPUT CHECKS
# --------------------------------------------------------------------------
# type checks
stopifnot(is(object, "BSFDataSet"))
stopifnot(is.logical(quiet))
# Check if none is specified
if (is.null(anno.annoDB) & is.null(anno.genes)) {
msg = paste0("None of the required annotation sources anno.annoDB or anno.genes was specified. ")
stop(msg)
}
# Check if both are specified
if (!is.null(anno.annoDB) & !is.null(anno.genes)) {
msg = paste0("Both of the required annotation sources anno.annoDB or anno.genes are specified. Please provide only one of the two. ")
stop(msg)
}
# Checks if anno.annoDB should be used
if (!is.null(anno.annoDB) & is.null(anno.genes)) {
stopifnot(is(anno.annoDB, "OrganismDb"))
datasource = "anno.annoDB"
# extract relevant annotation
anno.genes = genes(anno.annoDB, columns = c("GENEID"))
# Create matching vectors for columns from input annotation
selectID = as.character(anno.genes$GENEID)
}
# Checks if anno.genes should be used
if (is.null(anno.annoDB) & !is.null(anno.genes)) {
stopifnot(is(anno.genes, "GenomicRanges"))
datasource = "anno.genes"
# extract relevant annotation
anno.genes = anno.genes
# check correct annotation columns
annoColNames = colnames(mcols(anno.genes))
if (!all(match.geneID %in% annoColNames)) {
msg = paste0("The provided matching column for the geneID (",
match.geneID, ") is not present in the provided annotation. \n")
stop(msg)
}
# Create matching vectors for columns from input annotation
selectID = mcols(anno.genes)[match(match.geneID, colnames(mcols(anno.genes)))][[1]]
}
# check if cutoff was estimated before or if default should be used
if (!is.null(object@params$bsSize) &
!is.null(object@params$geneFilter)) {
cutoff = object@params$geneFilter
}
# handle options (remove, keep)
overlaps = match.arg(overlaps, choices = c("keepSingle", "removeAll", "keepAll"))
# Check for correct score column
# get all ranges
rngInitial = getRanges(object)
metaNames = colnames(mcols(rngInitial))
if (! match.score %in% metaNames) {
# if (! any(metaNames == match.score)) {
msg = paste0("Specified match.score (", match.score, "), is present in the ranges. Add respective column. ")
stop(msg)
}
# ---
# Store function parameters in list
optstr = list(cutoff = cutoff, overlaps = overlaps,
source = datasource, match.geneID = match.geneID)
object@params$geneWiseFilter = optstr
# MAIN COMPUTE
# --------------------------------------------------------------------------
# Create matching vectors for score column from input source
# --------------------------------------------------------------------------
score = mcols(rngInitial)[match(match.score, colnames(mcols(rngInitial)))][[1]]
rngInitial$score = score
# get all intergenic cases
rngIntergeneic = rngInitial[countOverlaps(rngInitial, anno.genes) == 0]
mcols(rngIntergeneic)$geneID = "Intergenic"
# get all cases to solve
# rngToSolve = subset(rngInitial, !rngInitial %in% rngIntergeneic)
# rngToSolve = rngInitial[!rngInitial %in% rngIntergeneic]
rngToSolve = rngInitial[countOverlaps(rngInitial, rngIntergeneic) == 0]
# Gene specific cutoffs
# ---------------------
# set cutpoints
potentialCutoffs = seq(0, 1, by = 0.01)
# group scores by gene range
ols = findOverlaps(anno.genes, rngToSolve) %>%
as.data.frame() %>%
mutate(score = rngToSolve$score[subjectHits]) %>%
group_by(queryHits)
# calculate gene range specific quantile cutoff
quants = ols %>%
dplyr::reframe(quantiles = quantile(score, probs = potentialCutoffs)) %>%
group_by(queryHits) %>%
mutate(cutoffs = potentialCutoffs) %>%
dplyr::filter(as.character(cutoffs) == as.character(cutoff))
# apply cutoff to every PureCLIP peak per gene
peakIdToKeep = ols
peakIdToKeep$quantiles = quants$quantiles[match(peakIdToKeep$queryHits, quants$queryHits)]
peakIdToKeep = peakIdToKeep %>%
dplyr::filter(score >= quantiles)
# add gene information
filteredPerRegion = rngToSolve[peakIdToKeep$subjectHits]
mcols(filteredPerRegion)$geneID = selectID[peakIdToKeep$queryHits]
# ---
# Store for plotting
dfPlot = countOverlaps(filteredPerRegion) %>% table() %>% as.data.frame() %>% rename("#N overlaps" = ".")
object@plotData$geneWiseFilter$data = dfPlot
# Deal with multiple loci
# -----------------------
# check for overlapping loci peaks
potentialOverlaps = findOverlaps(filteredPerRegion)
overlappingLociHits = potentialOverlaps[queryHits(potentialOverlaps) != subjectHits(potentialOverlaps)]
# reduce to a single pair
singlePair = overlappingLociHits %>%
as.data.frame() %>%
dplyr::filter(queryHits > subjectHits)
overlappingLociHits = overlappingLociHits[queryHits(overlappingLociHits) %in% singlePair$queryHits]
# get numbers
overlappingLociGenePairs = paste0(unique(names(filteredPerRegion[queryHits(overlappingLociHits)])),
",",
unique(names(filteredPerRegion[subjectHits(overlappingLociHits)])), " | ")
fractionOfRangesRemoved = paste0(round(length(overlappingLociHits) / length(rngInitial) * 100, digits = 2), "%")
# Handle cases
if (length(overlappingLociHits) > 0) {
msg0 = paste0(fractionOfRangesRemoved, " (", length(overlappingLociHits),"/", length(rngInitial), ")",
" peaks overlap with multiple anno.genes in the given gene annotation. \n")
if (overlaps == "keepSingle") {
msg1 = "A single instance of each peak is kept. This is recommended. \n "
filteredPerRegion = filteredPerRegion[!duplicated(filteredPerRegion)]
# report message
if (!quiet) message(c(msg0, msg1))
}
if (overlaps == "removeAll") {
msg1 = "All peaks on any duplicated range are removed. This is not recommended, since many peaks are lost. \n"
filteredPerRegion = filteredPerRegion[countOverlaps(filteredPerRegion) == 1]
warning(c(msg0, msg1))
}
if (overlaps == "keepAll") {
msg1 = "Duplicated peaks due to overlapping loci are not removed. This is not recommended, since peak numbers are inflated. \n"
warning(msg0, msg1)
}
}
# add intergenic cases again
filteredPerRegion = c(filteredPerRegion, rngIntergeneic)
# ---
# Store for results
resultLine = data.frame(
funName = "pureClipGeneWiseFilter()", class = "crosslink sites",
nIn = length(rngInitial), nOut = length(filteredPerRegion),
per = paste0(round(length(filteredPerRegion)/ length(rngInitial), digits = 2)*100,"%"),
options = paste0("Cutoff=", optstr$cutoff, ", overlapp option=", optstr$overlaps)
)
object@results = rbind(object@results, resultLine)
objectOut = setRanges(object, filteredPerRegion, quiet = quiet)
return(objectOut)
}
#' Assign binding sites to their hosting genes
#'
#' Function that assigns each binding site in the \code{\link{BSFDataSet}} to its
#' hosting gene given a gene annotation (\code{anno.annoDB}, \code{anno.genes}).
#'
#' Regardless of the annotation source that is being used, the respective meta
#' information about the genes have to be present. They can be set by the
#' \code{match.geneID}, \code{match.geneName} and \code{match.geneType} arguments.
#'
#' In the case of overlapping gene annotation, a single binding site will be
#' associated with multiple genes. The \code{\link{overlaps}} parameter allows
#' to decide in these cases. Option `frequency` will take the most frequently
#' observed gene type, option `hierarchy` works in conjunction with a user defined
#' rule (\code{overlaps.rule}). Options `remove` and `keep` will remove or
#' keep all overlapping cases, respectively.
#'
#' Note that if an overlaps exists, but gene types are identical options
#' `frequency` and `hierarchy` will cause the gene that was seen first to be
#' selected as representative.
#'
#' The function is part of the standard workflow performed by \code{\link{BSFind}}.
#'
#' @param object a \code{\link{BSFDataSet}} object with stored binding sites. This
#' means that ranges should be > 1
#' @param overlaps character; how overlapping gene loci should be handled.
#' @param overlaps.rule character vector; a vector of gene type that should
#' be used to handle overlapping cases in a hierarchical manor. The order of the
#' vector is the order of the hierarchy.
#' @param anno.annoDB an object of class \code{OrganismDbi} that contains
#' the gene annotation (!!! Experimental !!!).
#' @param anno.genes an object of class \code{\link{GenomicRanges}} that represents
#' the gene ranges directly
#' @param match.geneID character; meta column name of the gene ID
#' @param match.geneName character; meta column name of the gene name
#' @param match.geneType character; meta column name of the gene type
#' @param quiet logical; whether to print messages
#'
#' @return an object of class \code{\link{BSFDataSet}} with binding sites having
#' hosting gene information added to their meta columns.
#'
#' @seealso \code{\link{BSFind}}, \code{\link{geneOverlapsPlot}},
#' \code{\link{targetGeneSpectrumPlot}}
#'
#' @importFrom dplyr select slice_head group_by count ungroup arrange desc left_join
#' @importFrom GenomicFeatures genes
#' @importFrom dplyr desc
#'
#' @examples
#' # load clip data
#' files <- system.file("extdata", package="BindingSiteFinder")
#' load(list.files(files, pattern = ".rda$", full.names = TRUE))
#' # Load GRanges with genes
#' load(list.files(files, pattern = ".rds$", full.names = TRUE)[1])
#' bds = makeBindingSites(object = bds, bsSize = 9)
#' bds = assignToGenes(bds, anno.genes = gns)
#'
#' @export
assignToGenes <- function(object,
overlaps = c("frequency", "hierarchy", "remove", "keep"),
overlaps.rule = NULL,
anno.annoDB = NULL,
anno.genes = NULL,
match.geneID = "gene_id",
match.geneName = "gene_name",
match.geneType = "gene_type",
quiet = FALSE
) {
# local variables
datasource <- geneIndex <- geneID <- geneName <- value <- bsIndex <- Freq <- geneType <- match.missing <- NULL
selectType <- selectID <- selectName <- NULL
# INPUT CHECKS
# --------------------------------------------------------------------------
stopifnot(is(object, "BSFDataSet"))
stopifnot(is.logical(quiet))
# Check if none is specified
if (is.null(anno.annoDB) & is.null(anno.genes)) {
msg = paste0("None of the required annotation sources anno.annoDB or anno.genes was specified. ")
stop(msg)
}
# Check if both are specified
if (!is.null(anno.annoDB) & !is.null(anno.genes)) {
msg = paste0("Both of the required annotation sources anno.annoDB or anno.genes are specified. Please provide only one of the two. ")
stop(msg)
}
# Checks if anno.annoDB should be used
if (!is.null(anno.annoDB) & is.null(anno.genes)) {
stopifnot(is(anno.annoDB, "OrganismDb"))
if (!is.null(anno.genes)) {
msg = paste0("Parameter anno.annoDB and anno.genes are specified at the same time. Use only one of them.")
stop(msg)
} else {
datasource = "anno.annoDB"
# extract relevant annotation
anno.genes = genes(anno.annoDB, columns=c("ENSEMBL", "GENETYPE", "SYMBOL", "GENEID"))
# Create matching vectors for columns from input annotation
# --------------------------------------------------------------------------
selectID = as.character(anno.genes$GENEID)
selectName = as.character(anno.genes$SYMBOL)
selectType = as.character(anno.genes$GENETYPE)
}
}
# Checks if anno.genes should be used
if (is.null(anno.annoDB) & !is.null(anno.genes)) {
stopifnot(is(anno.genes, "GenomicRanges"))
if (!is.null(anno.annoDB)) {
msg = paste0("Parameter anno.annoDB and anno.genes are specified at the same time. Use only one of them.")
stop(msg)
} else {
datasource = "anno.genes"
# extract relevant annotation
# check for duplicated annotation columns
if (sum(colnames(mcols(anno.genes)) == match.geneID) > 1) {
msg = paste0("The names of multiple columns of the annotation match the match.geneID parameter. Please use a unique column name for match.geneID in anno.genes.")
stop(msg)
}
if (sum(colnames(mcols(anno.genes)) == match.geneName) > 1) {
msg = paste0("The names of multiple columns of the annotation match the match.geneName parameter. Please use a unique column name for match.geneName in anno.genes.")
stop(msg)
}
if (sum(colnames(mcols(anno.genes)) == match.geneType) > 1) {
msg = paste0("The names of multiple columns of the annotation match the match.geneType parameter. Please use a unique column name for match.geneType in anno.genes.")
stop(msg)
}
# check correct annotation columns
inNames = c(match.geneID, match.geneName, match.geneType)
annoColNames = colnames(mcols(anno.genes))
presentNames = inNames[inNames %in% annoColNames]
# extract columns from annotaiton based on available names
present.cols = lapply(presentNames, function(x){
match.col = mcols(anno.genes)[match(x, colnames(mcols(anno.genes)))][[1]]
return(match.col)
})
names(present.cols) = presentNames
# extract values from annotation for all present columns
if (match.geneID %in% names(present.cols)) {
selectID = present.cols[[match(match.geneID, names(present.cols))]]
} else {
msg = paste0("No meta column for ", match.geneID, " present. Creating custom ID.\n")
if (!quiet) message(msg)
selectID = paste0("CUSTOM", seq_along(anno.genes))
}
if (match.geneName %in% names(present.cols)) {
selectName = present.cols[[match(match.geneName, names(present.cols))]]
} else {
msg = paste0("No meta column for ", match.geneName, " present.\n")
if (!quiet) message(msg)
}
if (match.geneType %in% names(present.cols)) {
selectType = present.cols[[match(match.geneType, names(present.cols))]]
} else {
msg = paste0("No meta column for ", match.geneType, " present.\n")
if (!quiet) message(msg)
}
}
}
# handle options (hierarchy, frequency, remove, keep)
overlaps = match.arg(overlaps, choices = c("frequency", "hierarchy", "remove", "keep"))
if (is.null(selectType)) {
# gene type is missing
# -> check for correct overlap option
if (overlaps == "frequency" | overlaps == "hierarchy") {
msg0 = paste0(match.geneType, " is missing in the annotation. Options 'frequency' and 'hierarchy' are not available.\n")
msg1 = paste0("Select one of: 'remove' or 'keep'.\n")
stop(c(msg0, msg1))
}
}
# Check multiple loci options
if (overlaps == "hierarchy" & is.null(overlaps.rule)) {
msg1 = paste0("Binding sites on anno.genes with overlapping annotaitons is set to be handled by 'hierarchy', but no rule is provided. \n")
msg2 = paste0("Change the 'overlaps' option or provide a valid rule. \n")
stop(c(msg1, msg2))
}
# --------------------------------------------------------------------------
# MAIN
# --------------------------------------------------------------------------
# ---
# Store function parameters in list
optstr = list(source = datasource, match.geneID = match.geneID,
match.geneName = match.geneName, match.geneType = match.geneType,
overlaps = overlaps, overlaps.rule = overlaps.rule)
object@params$assignToGenes = optstr
# Match Binding sites and anno.genes
# --------------------------------------------------------------------------
# split ranges in duplicated event and non duplicated events
rngInitial = getRanges(object)
rngInitial$currIdx = seq_along(rngInitial)
# remove possibility of partial overlaps by setting `minoverlap` parameter
# to the size of the binding sites
# -> use the minimum of the size for uneven binding sites sizes
# -> this is not optimal since partial overlaps can be re-introduced
# -> does not affect ranges that are true binding sites
bs.set.range = min(unique(width(rngInitial)))
ols = findOverlaps(anno.genes, rngInitial, minoverlap = bs.set.range)
# ols = findOverlaps(anno.genes, rngInitial)
# Deal with multiple loci
# --------------------------------------------------------------------------
# Count out how many cases there are
countOlsSize = as.data.frame(table(duplicated(subjectHits(ols))))
if (nrow(countOlsSize) == 1){
# no duplications are present
# -> no overlaps have to be resolved
totalBS = countOlsSize$Freq[1]
duplicatedBS = 0
} else {
totalBS = countOlsSize$Freq[1]
duplicatedBS = countOlsSize$Freq[2]
}
duplicatedFraction = paste0(round(duplicatedBS / totalBS * 100, digits = 2), "%")
# ----
# Store results for plotting
# -> select data
if (all(!is.null(selectID) & !is.null(selectType) & !is.null(selectName))) {
# nothing is missing
dfPlot = data.frame(geneIndex = queryHits(ols), bsIndex = subjectHits(ols),
geneType = selectType[queryHits(ols)]) %>% #TODO type
dplyr::mutate(value = 1) %>%
dplyr::select(-geneIndex) %>%
tidyr::pivot_wider(names_from = geneType, values_from = value, values_fn = length, values_fill = 0) %>%
dplyr::select(-bsIndex) %>%
dplyr::mutate_all(~ ifelse(. > 1, 1, .))
} else {
# check what is missing
if (!is.null(selectType)) {
# anything other than gene type is missing
dfPlot = data.frame(geneIndex = queryHits(ols), bsIndex = subjectHits(ols),
geneType = selectType[queryHits(ols)]) %>% #TODO type
dplyr::mutate(value = 1) %>%
dplyr::select(-geneIndex) %>%
tidyr::pivot_wider(names_from = geneType, values_from = value, values_fn = length, values_fill = 0) %>%
dplyr::select(-bsIndex) %>%
dplyr::mutate_all(~ ifelse(. > 1, 1, .))
} else {
# gene type is missing
dfPlot = NULL
}
}
object@plotData$assignToGenes$dataOverlaps = dfPlot
# Case where there are no overlaps
if (duplicatedBS == 0) {
msg = paste0("No binding sites (",
format(totalBS, big.mark = ',', decimal.mark = "."),
") on overlapping gene loci. \n")
msg2 = paste0("Parameter 'overlaps' set to '",
overlaps, "' not in use. \n")
if(!quiet) message(c(msg, msg2))
}
# Cases where there are multiple overlaps
msg0 = paste0(duplicatedFraction, " (", duplicatedBS,"/", totalBS, ")",
" of binding sites overlap with multiple anno.genes in the given gene annotation. \n")
if (overlaps == "hierarchy") {
# apply hierarchy solution
msg1 = paste0("Apply 'hierarchy' solution")
if(!quiet) message(c(msg0, msg1))
ruleMod = data.frame(gene_type = overlaps.rule, idx = seq_along(overlaps.rule))
rngResolved = .resolveGeneOverlapsWithRule(rng = rngInitial, ols = ols, rule = ruleMod,
selectID = selectID, selectName = selectName,
selectType = selectType)
if (length(rngInitial) != length(rngResolved)) {
if(!quiet) message(paste0(format(length(rngInitial)-length(rngResolved), big.mark = ",", decimal.mark = "."),
" binding sites could not be assigned to a gene range.\n"))
}
}
if (overlaps == "frequency") {
# apply frequency solution
msg1 = paste0("Apply 'frequency' solution")
if (!quiet) message(c(msg0, msg1))
ruleFreq = selectType %>% table() %>% as.data.frame() %>%
rename(gene_type = 1) %>% arrange(desc(Freq)) %>%
mutate(idx = row_number()) %>% rename(n = 2)
rngResolved = .resolveGeneOverlapsWithRule(rng = rngInitial, ols = ols, rule = ruleFreq,
selectID = selectID, selectName = selectName, selectType = selectType)
if (length(rngInitial) != length(rngResolved)) {
if (!quiet) message(paste0(format(length(rngInitial)-length(rngResolved), big.mark = ",", decimal.mark = "."),
" binding sites could not be assigned to a gene range.\n"))
}
}
if (overlaps == "remove") {
# apply remove option
msg1 = "Binding sites from overlapping loci are removed. This is not recommended. Please see options 'heirarchy' and 'frequency'. \n"
warning(c(msg0, msg1))
rngResolved = rngInitial[subjectHits(ols)]
# manage meta columns
if (!is.null(selectID)) {
mcols(rngResolved) = cbind(mcols(rngResolved), geneID = selectID[queryHits(ols)])
}
if (!is.null(selectType)) {
mcols(rngResolved) = cbind(mcols(rngResolved), geneType = selectType[queryHits(ols)])
}
if (!is.null(selectName)) {
mcols(rngResolved) = cbind(mcols(rngResolved), geneName = selectName[queryHits(ols)])
}
# remove duplicates
rngResolved = rngResolved[countOverlaps(rngResolved) == 1]
}
if (overlaps == "keep") {
msg1 = "Binding sites from overlapping loci are kept. This is not recommended. Please see options 'heirarchy' and 'frequency'. \n"
warning(c(msg0, msg1))
rngResolved = rngInitial[subjectHits(ols)]
# manage meta columns
if (!is.null(selectID)) {
mcols(rngResolved) = cbind(mcols(rngResolved), geneID = selectID[queryHits(ols)])
}
if (!is.null(selectType)) {
mcols(rngResolved) = cbind(mcols(rngResolved), geneType = selectType[queryHits(ols)])
}
if (!is.null(selectName)) {
mcols(rngResolved) = cbind(mcols(rngResolved), geneName = selectName[queryHits(ols)])
}
}
# ---
# Store results for plotting
if (all(!is.null(selectID) & !is.null(selectType) & !is.null(selectName))) {
# nothing is missing
nCountGenes = rngResolved %>% mcols() %>% as.data.frame() %>%
group_by(geneID) %>% dplyr::count(geneType) %>% ungroup() %>%
dplyr::count(geneType) %>% arrange(desc(n))
nCountBs = rngResolved %>% mcols() %>% as.data.frame() %>%
dplyr::count(geneType) %>% arrange(desc(n))
nCountAll = left_join(nCountBs, nCountGenes, by = "geneType") %>%
rename("nBs" = "n.x", "nGenes" = "n.y")
} else {
# check what is missing
if (!is.null(selectType)) {
# anything other than gene type is missing
nCountGenes = rngResolved %>% mcols() %>% as.data.frame() %>%
group_by(geneID) %>% dplyr::count(geneType) %>% ungroup() %>%
dplyr::count(geneType) %>% arrange(desc(n))
nCountBs = rngResolved %>% mcols() %>% as.data.frame() %>%
dplyr::count(geneType) %>% arrange(desc(n))
nCountAll = left_join(nCountBs, nCountGenes, by = "geneType") %>%
rename("nBs" = "n.x", "nGenes" = "n.y")
}
else {
# gene type is missing
nCountAll = NULL
}
}
object@plotData$assignToGenes$dataSpectrum = nCountAll
# ---
# Store for results
resultLine = data.frame(
funName = "assignToGenes()", class = "binding sites",
nIn = length(rngInitial), nOut = length(rngResolved),
per = paste0(round(length(rngResolved)/ length(rngInitial), digits = 2)*100,"%"),
options = paste0("Overlaps=", optstr$overlaps, ", Source=", optstr$source,
ifelse(!rlang::is_empty(optstr$overlaps.rule),
paste0(", overlaps.rule=",
paste(optstr$overlaps.rule, collapse = ">")), ""))
)
object@results = rbind(object@results, resultLine)
objectOut = setRanges(object, rngResolved, quiet = quiet)
return(objectOut)
}
#' Assign binding sites to their hosting transcript regions
#'
#' Function that assigns each binding site in the \code{\link{BSFDataSet}} to its
#' hosting transcript region given an annotation database (\code{anno.annoDB}), or
#' a GRanges list / \code{\link{CompressedGRangesList}} (\code{anno.transcriptRegionList})
#' that holds the ranges for the transcript regions of interest.
#'
#' Since the assignment is based on the overlaps of annotated transcript ranges
#' and binding sites, no matching meta data is needed.
#'
#' In the case of transcript regions overlaps are very frequent. To resolve these
#' cases the \code{\link{overlaps}} argument can be used. Option `frequency`
#' will take the most frequently observed transcript region, option `hierarchy`
#' works in conjunction with a user defined rule (\code{overlaps.rule}).
#' Options `flag` and `remove` will label binding sites with an ambiguous tag or
#' remove all overlapping cases, respectively.
#'
#' The function is part of the standard workflow performed by \code{\link{BSFind}}.
#'
#' @param object a \code{\link{BSFDataSet}} object with stored binding sites. This
#' means that ranges should be > 1
#' @param overlaps character; how overlapping transcript regions should be handled.
#' @param overlaps.rule character vector; a vector of transcript region
#' names that should be used to handle overlapping cases in a hierarchical manor.
#' The order of the vector is the order of the hierarchy.
#' @param anno.annoDB an object of class \code{OrganismDbi} that contains
#' the transcript region annotation (!!! Experimental !!!).
#' @param anno.transcriptRegionList an object of class \code{\link{CompressedGRangesList}}
#' that holds an ranges for each transcript region
#'
#' @param normalize.exclude.upper numeric; percentage value that indicates the
#' upper boundary for transcript region length to be considered when calculating
#' normalization factors for regions.
#' @param normalize.exclude.lower numeric; percentage value that indicates the
#' lower boundary for transcript region length to be considered when calculating
#' normalization factors for regions.
#'
#' @param quiet logical; whether to print messages
#'
#' @return an object of class \code{\link{BSFDataSet}} with binding sites having
#' hosting transcript region information added to their meta columns.
#'
#' @seealso \code{\link{BSFind}},
#' \code{\link{transcriptRegionOverlapsPlot}},
#' \code{\link{transcriptRegionSpectrumPlot}}
#'
#' @importFrom dplyr rename_with desc
#' @importFrom GenomicFeatures cds intronsByTranscript threeUTRsByTranscript fiveUTRsByTranscript
#'
#' @examples
#' # load clip data
#' files <- system.file("extdata", package="BindingSiteFinder")
#' load(list.files(files, pattern = ".rda$", full.names = TRUE))
#' load(list.files(files, pattern = ".rds$", full.names = TRUE)[1])
#' load(list.files(files, pattern = ".rds$", full.names = TRUE)[2])
#'
#' bds = makeBindingSites(object = bds, bsSize = 9)
#' bds = assignToGenes(bds, anno.genes = gns)
#' bds = assignToTranscriptRegions(object = bds, anno.transcriptRegionList = regions)
#'
#' @export
assignToTranscriptRegions <- function(object, # bindingSiteFinder
overlaps = c("frequency", "hierarchy", "flag", "remove"), # overlapping loci solution options
overlaps.rule = NULL, # rule to apply when option hierarchy is selected
anno.annoDB = NULL, # annotation data base of class OrganismDbi; can be NULL -> but then anno.transcriptRegionList must be provided
anno.transcriptRegionList = NULL, # GRangesList object with the transcript features to be used; can be NULL -> but then anno.annoDB must be provided
normalize.exclude.upper = 0.02,
normalize.exclude.lower = 0.02,
quiet = FALSE
) {
# local variables
datasource <- `.` <- Freq <- transcript_type <- NULL
# INPUT CHECKS
# --------------------------------------------------------------------------
stopifnot(is(object, "BSFDataSet"))
stopifnot(is.logical(quiet))
# Check if none is specified
if (is.null(anno.annoDB) & is.null(anno.transcriptRegionList)) {
msg = paste0("None of the required annotation sources anno.annoDB or anno.transcriptRegionList was specified. ")
stop(msg)
}
# Check if both are specified
if (!is.null(anno.annoDB) & !is.null(anno.transcriptRegionList)) {
msg = paste0("Both of the required annotation sources anno.annoDB or anno.transcriptRegionList are specified. Please provide only one of the two. ")
stop(msg)
}
# anno.annoDB should be used
if (!is.null(anno.annoDB) & is.null(anno.transcriptRegionList)) {
stopifnot(is(anno.annoDB, "OrganismDb"))
if (!is.null(anno.transcriptRegionList)) {
msg = paste0("Parameter anno.annoDB and anno.transcriptRegionList are specified at the same time. Use only one of them.")
stop(msg)
} else {
datasource = "anno.annoDB"
# extract relevant annotation - transcripts
cdseq = cds(anno.annoDB)
intrns = unlist(intronsByTranscript(anno.annoDB))
utrs3 = unlist(threeUTRsByTranscript(anno.annoDB))
utrs5 = unlist(fiveUTRsByTranscript(anno.annoDB))
anno.transcriptRegionList = GRangesList(CDS = cdseq, Intron = intrns, UTR3 = utrs3, UTR5 = utrs5)
names(anno.transcriptRegionList) = toupper(names(anno.transcriptRegionList))
}
}
# anno.transcriptRegionList should be used
if (is.null(anno.annoDB) & !is.null(anno.transcriptRegionList) ) {
stopifnot(is(anno.transcriptRegionList, "GenomicRangesList"))
if (!is.null(anno.annoDB)) {
msg = paste0("Parameter anno.annoDB and anno.transcriptRegionList are specified at the same time. Use only one of them.")
stop(msg)
} else {
datasource = "anno.transcriptRegionList"
# extract relevant annotation - transcripts
anno.transcriptRegionList = anno.transcriptRegionList
names(anno.transcriptRegionList) = toupper(names(anno.transcriptRegionList))
}
}
# handle options (hierarchy, frequency, remove, keep)
overlaps = match.arg(overlaps, choices = c("frequency", "hierarchy", "flag", "remove"))
# Check multiple loci options
if (overlaps == "hierarchy" & is.null(overlaps.rule)) {
msg1 = paste0("Multiple region problem is set to be handled by 'hierarchy', but no overlaps.rule is provided. \n")
msg2 = paste0("Change the 'overlaps' option or provide a valid overlaps.rule. \n")
stop(c(msg1, msg2))
}
overlaps.rule = toupper(overlaps.rule)
if (! all(overlaps.rule %in% names(anno.transcriptRegionList))){
msg = paste0("Regions defined in overlaps.rule (", paste(overlaps.rule, collapse = " > "),
") does not match the input region names (", paste(names(anno.transcriptRegionList), collapse = ", "), ")")
stop(msg)
}
if (overlaps == "hierarchy" & is.null(overlaps.rule)) {
msg1 = paste0("Overlap are set to be solved by 'hierarchy', but no overlaps.rule is provided. \n")
msg2 = paste0("Change the 'overlaps' option or provide a valid overlaps.rule. \n")
stop(c(msg1, msg2))
}
# ---
# Store function parameters in list
optstr = list(source = datasource, overlaps = overlaps, overlaps.rule = overlaps.rule)
object@params$assignToTranscriptRegions = optstr
# MAIN COMPUTE
# --------------------------------------------------------------------------
# Match transcript annotation with binding sites
# --------------------------------------------------------------------------
rngInitial = getRanges(object)
rng = getRanges(object)
cRange = resize(rng, width = 1, fix = "center")
countDf = lapply(anno.transcriptRegionList, function(x) {
countOverlaps(cRange, x)
}) %>% as.data.frame() %>% rename_with(toupper)
# Add intergenic cases
countDf = countDf %>% mutate(OTHER = ifelse(rowSums(.) == 0, 1, 0))
# ---
# Store results for plotting
plotDf = countDf
plotDf[plotDf >= 1] = 1
object@plotData$assignToTranscriptRegions$dataOverlaps = plotDf
# Deal with multiple loci
# --------------------------------------------------------------------------
# Count out how many cases there are
totalBs = nrow(countDf)
numberOfLociHit = rowSums(countDf != 0)
ambigousBs = nrow(countDf[numberOfLociHit != 1,])
ambigousBsFraction = paste0(round(ambigousBs / totalBs * 100, digits = 2), "%")
# Case where there are no overlaps
if (ambigousBs == 0) {
msg0 = paste0("No binding sites (", format(totalBs, big.mark = ',', decimal.mark = "."), ") on overlapping transcript loci. \n")
msg1 = paste0("Parameter 'overlaps' set to '", overlaps, "' not in use. \n")
if(!quiet) message(c(msg0, msg1))
# assign regions
countDfSorted = countDf
countDfSorted[countDfSorted > 0] = 1
col_names = colnames(countDfSorted)
transcriptRegion = ifelse(rowSums(countDfSorted) == 1, col_names[max.col(countDfSorted, ties.method = "first")], "Ambiguous")
mcols(rng) = cbind.data.frame(mcols(rng), transcriptRegion = transcriptRegion)
}
# Cases where we have to resolve overlapping loci
if (ambigousBs > 0) {
msg0 = paste0(ambigousBsFraction, " (", ambigousBs,"/", totalBs, ")",
" of binding sites overlap with multiple different transcript regions in the given annotation. \n")
if (overlaps == "hierarchy") {
# apply hierarchy solution
msg1 = paste0("Apply 'hierarchy' solution. \n ")
if(!quiet) message(c(msg0, msg1))
countDfSorted = countDf[, overlaps.rule]
countDfSorted[countDfSorted > 0] = 1
transcriptRegion = names(countDfSorted)[max.col(countDfSorted, ties.method = "first")]
mcols(rng) = cbind.data.frame(mcols(rng), transcriptRegion = transcriptRegion)
}
if (overlaps == "frequency") {
msg1 = paste0("Apply 'frequency' solution. \n ")
if(!quiet) message(c(msg0, msg1))
countDfSorted = countDf
# countDfSorted[countDfSorted > 0] = 1
# get overlaps.rule by frequency
overlaps.rule = names(countDfSorted)[max.col(countDfSorted, ties.method = "first")] %>%
table() %>% as.data.frame() %>% rename(transcript_type = 1) %>%
arrange(desc(Freq)) %>% pull(transcript_type) %>% as.character()
# apply frequency solution
countDfSorted = countDfSorted[, overlaps.rule]
transcriptRegion = names(countDfSorted)[max.col(countDfSorted, ties.method = "first")]
mcols(rng) = cbind.data.frame(mcols(rng), transcriptRegion = transcriptRegion)
}
if (overlaps == "flag") {
msg1 = paste0("Binding sites with multiple different overlapping transcript regions are being flaged with the tag 'Ambiguous'. \n")
if(!quiet) message(c(msg0, msg1))
countDfSorted = countDf
countDfSorted[countDfSorted > 0] = 1
# flag multiple loci binding sites
col_names = colnames(countDfSorted)
transcriptRegion = ifelse(rowSums(countDfSorted) == 1, col_names[max.col(countDfSorted, ties.method = "first")], "Ambiguous")
mcols(rng) = cbind.data.frame(mcols(rng), transcriptRegion = transcriptRegion)
# remove mulitple loci binding sitges
mcols(rng) = cbind.data.frame(mcols(rng), transcriptRegion = transcriptRegion)
}
if (overlaps == "remove") {
msg1 = paste0("Binding sites with multiple different overlapping transcript regions are being removed. This is not recommended! \n")
if(!quiet) message(c(msg0, msg1))
countDfSorted = countDf
countDfSorted[countDfSorted > 0] = 1
# flag multiple loci binding sites
col_names = colnames(countDfSorted)
transcriptRegion = ifelse(rowSums(countDfSorted) == 1, col_names[max.col(countDfSorted, ties.method = "first")], "Ambiguous")
mcols(rng) = cbind.data.frame(mcols(rng), transcriptRegion = transcriptRegion)
# remove mulitple loci binding sitges
mcols(rng) = cbind.data.frame(mcols(rng), transcriptRegion = transcriptRegion)
rng = rng[rng$transcriptRegion != "Ambiguous"]
}
}
# ---
# Calculate width for length based normalization
## normalize by all expressed regions
debugList = list(object = object, rng = rng)
# saveRDS(debugList, file = "./debugList.rds")
norm.factors = .calcNormalizeFactors(object, rng, anno.transcriptRegionList,
normalize.exclude.lower = normalize.exclude.lower,
normalize.exclude.upper = normalize.exclude.upper)
# ---
# Store results for plotting
dfPlot = rng$transcriptRegion %>%
table() %>%
as.data.frame() %>%
rename("TranscriptRegion" = ".") %>%
arrange(desc(Freq)) %>%
left_join(., y = norm.factors, by = "TranscriptRegion")
object@plotData$assignToTranscriptRegions$dataSpectrum = dfPlot
# ---
# Store for results
resultLine = data.frame(
funName = "assignToTranscriptRegions()", class = "binding sites",
nIn = length(rngInitial), nOut = length(rng),
per = paste0(round(length(rng)/ length(rngInitial), digits = 2)*100,"%"),
options = paste0("Overlaps=", optstr$overlaps, ", Source=", datasource,
ifelse(!rlang::is_empty(optstr$overlaps.rule), paste0(", overlaps.rule=", paste(optstr$overlaps.rule, collapse = ">")), ""))
)
object@results = rbind(object@results, resultLine)
outObject = setRanges(object, rng, quiet = quiet)
return(outObject)
}
#' Function to estimate the appropriate binding site width together with the
#' optimal gene-wise filter level.
#'
#' This function tests different width of binding sites for different gene-wise
#' filtering steps. For each test the signal-to-score ratio is calculated. The
#' mean over all gene-wise filterings at each binding site width is used to
#' extract the optimal width, which serves as anchor to select the optimal
#' gene-wise filter.
#'
#' Parameter estimation is done on a subset of all crosslink sites
#' (\code{est.subsetChromosome}).
#'
#' Gene-level filter can be tested with varying
#' levels of accuracy ranging from `finest` to `coarse`, representing 1% and
#' 20% steps, respectively.
#'
#' Binding site computation at each step can be done on three different accuracy
#' level (\code{bsResolution}). Option `fine` is equal to a normal run
#' of the \code{\link{makeBindingSites}} function. `medium` will perform
#' a shorter version of the binding site computation, skipping some of the
#' refinement steps. Option `coarse` will approximate binding sites by merged
#' crosslinks regions, aligning the center at the site with the highest score.
#'
#' For each binding site in each set given the defined resolutions a signal-to-
#' flank score ratio is calculated and the mean of this score per set is returned.
#' Next a mean of means is created which results in a single score for each
#' binding site width that was tested. The width that yielded the highest score
#' is selected as optimal. In addtion the \code{minimumStepGain} option
#' allows control over the minimum additional gain in the score that a tested
#' width has to have to be selected as the best option.
#'
#' To enhance the sensitivity of the binding site estimation, the sensitivity
#' mode exists. In this mode crosslink sites undergo a pre-filtering and merging
#' step, to exclude potential artifical peaks (experimental-, mapping-biases).
#' If sensitivity mode is activated the \code{est.minWidth} option should be set
#' to 1.
#'
#' The optimal geneFilter is selected as the first one that passes the merged
#' mean of the selected optimal binding site width.
#'
#' The function is part of the standard workflow performed by \code{\link{BSFind}}.
#'
#' @param object a \code{\link{BSFDataSet}} object with stored crosslink sites.
#' This means that ranges should have a width = 1.
#' @param bsResolution character; level of resolution at which different binding
#' site width should be tested
#' @param geneResolution character; level of resolution at which gene-wise filtering
#' steps should be tested
#' @param est.maxBsWidth numeric; the largest binding site width which should
#' considered in the testing
#' @param est.minimumStepGain numeric; the minimum additional gain in the score
#' in percent the next binding site width has to have, to be selected as best option
#' @param est.maxSites numeric; maximum number of PureCLIP sites that are used;
#' @param est.subsetChromosome character; define on which chromosome the
#' estimation should be done in function \code{\link{estimateBsWidth}}
#' @param est.minWidth the minimum size of regions that are subjected to the
#' iterative merging routine, after the initial region concatenation.
#' @param est.offset constant added to the flanking count in the signal-to-flank
#' ratio calculation to avoid division by Zero
#'
#' @param sensitive logical; whether to enable sensitive pre-filtering before
#' binding site merging or not
#' @param sensitive.size numeric; the size (in nucleotides) of the merged
#' sensitive region
#' @param sensitive.minWidth numeric; the minimum size (in nucleoties) of the
#' merged sensitive region
#'
#' @param anno.annoDB an object of class \code{OrganismDbi} that contains
#' the gene annotation (!!! Experimental !!!).
#' @param anno.genes an object of class \code{\link{GenomicRanges}} that represents
#' the gene ranges directly
#' @param bsResolution.steps numeric vector; option to use a user defined threshold
#' for binding site width directly. Overwrites \code{bsResolution}
#' @param geneResolution.steps numeric vector; option to use a user defined threshold
#' vector for gene-wise filtering resolution. Overwrites \code{geneResolution}
#' @param quiet logical; whether to print messages
#' @param veryQuiet logical; whether to suppress all messages
#' @param reportScoresPerBindingSite report the ratio score for each binding site
#' separately. Warning! This is for debugging and testing only. Downstream
#' functions can be impaired.
#' @param ... additional arguments passed to \code{\link{pureClipGeneWiseFilter}}
#'
#' @return an object of class \code{\link{BSFDataSet}} with binding sites with
#' the `params` slots `bsSize` and `geneFilter` being filled
#'
#' @seealso \code{\link{BSFind}},
#' \code{\link{estimateBsWidthPlot}}
#'
#' @importFrom dplyr summarise slice_head slice_tail
#' @importFrom utils head
#' @importFrom GenomicFeatures genes
#' @importFrom stats sd
#'
#' @examples
#' # load clip data
#' files <- system.file("extdata", package="BindingSiteFinder")
#' load(list.files(files, pattern = ".rda$", full.names = TRUE))
#' load(list.files(files, pattern = ".rds$", full.names = TRUE)[1])
#' load(list.files(files, pattern = ".rds$", full.names = TRUE)[2])
#' estimateBsWidth(bds, anno.genes = gns, est.maxBsWidth = 19,
#' geneResolution = "coarse", bsResolution = "coarse", est.subsetChromosome = "chr22")
#'
#' @export
estimateBsWidth <- function(object, # BindingSiteFinder object
bsResolution = c("medium", "fine", "coarse"),
geneResolution = c("medium", "coarse", "fine", "finest"),
est.maxBsWidth = 13,
est.minimumStepGain = 0.02,
est.maxSites = Inf,
est.subsetChromosome = "chr1", #
est.minWidth = 2,
est.offset = 1,
sensitive = FALSE,
sensitive.size = 5,
sensitive.minWidth = 2,
anno.annoDB = NULL,
anno.genes = NULL,
bsResolution.steps = NULL,
geneResolution.steps = NULL,
quiet = TRUE,
veryQuiet = FALSE,
reportScoresPerBindingSite = FALSE,
...
) {
# initialize local variables
bsSize <- ms <- growth_per <- sel <- geneWiseFilter <- est.option <- signalToFlankRatio <- NULL
# INPUT CHECKS
# --------------------------------------------------------------------------
stopifnot(is(object, "BSFDataSet"))
stopifnot(is.logical(quiet))
if (c(est.maxBsWidth %% 2) == 0) {
stop("est.maxBsWidth is even. An odd number is required to have a distinct binding site center.")
}
# Check if none is specified
if (is.null(anno.annoDB) & is.null(anno.genes)) {
msg = paste0("None of the required annotation sources anno.annoDB or anno.genes was specified. ")
stop(msg)
}
# Check if both are specified
if (!is.null(anno.annoDB) & !is.null(anno.genes)) {
msg = paste0("Both of the required annotation sources anno.annoDB or anno.genes are specified. Please provide only one of the two. ")
stop(msg)
}
# Checks if anno.annoDB should be used
if (!is.null(anno.annoDB) & is.null(anno.genes)) {
stopifnot(is(anno.annoDB, "OrganismDb"))
if (!is.null(anno.genes)) {
msg = paste0("Parameter anno.annoDB and anno.genes are specified at the same time. Use only one of them.")
stop(msg)
} else {
# extract relevant annotation
anno.genes = genes(anno.annoDB)
}
}
# Checks if anno.genes should be used
if (is.null(anno.annoDB) & !is.null(anno.genes)) {
stopifnot(is(anno.genes, "GenomicRanges"))
if (!is.null(anno.annoDB)) {
msg = paste0("Parameter anno.annoDB and anno.genes are specified at the same time. Use only one of them.")
stop(msg)
} else {
# extract relevant annotation
anno.genes = anno.genes
}
}
# handle gene filter resolution
geneResolution = match.arg(geneResolution, choices = c("medium", "coarse", "fine", "finest"))
if (!is.null(geneResolution.steps)) {
geneResolution.steps = geneResolution.steps
} else {
if (geneResolution == "coarse") {
geneResolution.steps = seq(from = 0, to = 1, by = 0.2)
geneResolution.steps = geneResolution.steps[1:length(geneResolution.steps)-1]
}
if (geneResolution == "medium") {
geneResolution.steps = seq(from = 0, to = 1, by = 0.1)
geneResolution.steps = geneResolution.steps[1:length(geneResolution.steps)-1]
}
if (geneResolution == "fine") {
geneResolution.steps = seq(from = 0, to = 1, by = 0.05)
geneResolution.steps = geneResolution.steps[1:length(geneResolution.steps)-1]
}
if (geneResolution == "finest") {
geneResolution.steps = seq(from = 0, to = 1, by = 0.01)
geneResolution.steps = geneResolution.steps[1:length(geneResolution.steps)-1]
}
}
# handle binding site compute resolution
bsResolution = match.arg(bsResolution, choices = c("medium", "fine", "coarse"))
# handle est.maxBsWidth
if (!is.null(bsResolution.steps)) {
bsResolution.steps = bsResolution.steps
} else {
bsResolution.steps = seq(from = 3, to = est.maxBsWidth, by = 2)
}
# ---
# Store function parameters in list
optstr = list(bsResolution = bsResolution,
geneResolution = geneResolution,
est.maxBsWidth = est.maxBsWidth,
est.maxSites = est.maxSites,
est.minimumStepGain = est.minimumStepGain,
est.subsetChromosome = est.subsetChromosome,
est.minWidth = est.minWidth,
est.offset = est.offset,
sensitive = sensitive,
sensitive.size = sensitive.size,
sensitive.minWidth = sensitive.minWidth)
object@params$estimateBsWidth = optstr
# PREPARE TEST RANGES + SIGNAL
# --------------------------------------------------------------------------
checkRng = getRanges(object)
# check if subset is part of the seqnames from ranges
if (!all(est.subsetChromosome %in% levels(seqnames(checkRng)))){
msg = paste0("Chromosome to estimate on (", est.subsetChromosome, "), is not included in the ranges: ", paste(levels(seqnames(checkRng)), collapse = ", "))
stop(msg)
}
# limit estimation to a specific chromosome
if (!is.null(est.subsetChromosome)) {
redObj = .subsetByChr(object, chr = est.subsetChromosome, quiet = quiet)
} else {
redObj = object
}
# limit estimation to a maximum number of sites
estRng = getRanges(redObj)
if (length(estRng) > est.maxSites) {
estRng = head(estRng, est.maxSites)
redObj = setRanges(redObj, estRng, quiet = quiet)
}
# limit the clip signal to the ranges used for estimation (plus some extra offset)
if (!is.null(est.subsetChromosome)) {
# reduce signal to frame
maxFrame = ceiling((max(bsResolution.steps) *3))
redObj = .reduceSignalToFrame(redObj, frame = maxFrame, quiet = quiet)
} else {
# don't reduce signal
redObj = redObj
}
# collapse signal from replicates
sgnMerge = .collapseSamples(getSignal(redObj))
# MAIN COMPUTE
# --------------------------------------------------------------------------
# a counter that shows to how many percent the computation is done
counterTotalIterations = length(geneResolution.steps) * length(bsResolution.steps)
counterCurrentIterations = 0
# print start message
msg = paste0("Estimation at: ", round((counterCurrentIterations/counterTotalIterations)*100 ), "%" )
if(!veryQuiet) message(msg)
# calculate binding sites for each filter step and width
scoreAllDf = lapply(geneResolution.steps, function(bsFilterStep){
# bsFilterStep = 0.2
# apply current gene-wise filter (has two modes)
if (isTRUE(sensitive)) {
# -> sensitive mode
# define regions in which sites should be kept
cRng = getRanges(redObj)
cRngMerge = reduce(cRng, min.gapwidth = sensitive.size)
cRngMerge = cRngMerge[width(cRngMerge) >= sensitive.minWidth]
# currFilterObj = pureClipGeneWiseFilter(object = redObj, anno.genes = anno.genes, cutoff = bsFilterStep, quiet = quiet, ...)
currFilterObj = pureClipGeneWiseFilter(object = redObj, anno.genes = anno.genes, cutoff = bsFilterStep, quiet = quiet)
cRngFilter = getRanges(currFilterObj)
currRng = IRanges::subsetByOverlaps(cRngFilter, cRngMerge)
} else {
# -> normal mode
# currFilterObj = pureClipGeneWiseFilter(object = redObj, anno.genes = anno.genes, cutoff = bsFilterStep, quiet = quiet, ...)
currFilterObj = pureClipGeneWiseFilter(object = redObj, anno.genes = anno.genes, cutoff = bsFilterStep, quiet = quiet)
currRng = getRanges(currFilterObj)
}
# calculate binding sites for current bsWidth
rngPerWidth = lapply(bsResolution.steps, function(bsWidthStep){
if (bsResolution == "fine") {
# calculate with full binding sites
currBsObj = makeBindingSites(object = currFilterObj,
bsSize = bsWidthStep,
minWidth = est.minWidth,
quiet = quiet)
currBs = getRanges(currBsObj)
}
if (bsResolution == "medium") {
# approximate binding sites by a single merge and extend round
currBs = .approximateBindingSites_medium(rng = currRng,
sgn = sgnMerge,
bsSize = bsWidthStep,
minWidth = est.minWidth)
}
if (bsResolution == "coarse") {
# approximate binding sites by reduced pureclip sites center
currBs = .approximateBindingSites_coarse(rng = currRng,
bsSize = bsWidthStep,
minWidth = est.minWidth)
}
return(currBs)
})
rngPerWidth = unlist(GRangesList(rngPerWidth))
# handle all plus ranges
cRangePlus = subset(rngPerWidth, strand == "+")
if (length(cRangePlus) > 0) {
bsSumPlus = sum(sgnMerge$signalPlus[cRangePlus])
extendedRangePlus = cRangePlus + cRangePlus$bsSize
exSumPlus = sum(sgnMerge$signalPlus[extendedRangePlus])
# extendedRangePlus$signalToFlankRatio = (bsSumPlus) / (((exSumPlus - bsSumPlus) / 2) + est.offset)
extendedRangePlus$signalToFlankRatio = (bsSumPlus / exSumPlus)
} else {
extendedRangePlus = cRangePlus
}
# handle all minus ranges
cRangeMinus = subset(rngPerWidth, strand == "-")
if (length(cRangeMinus) > 0) {
bsSumMinus = sum(sgnMerge$signalMinus[cRangeMinus])
extendedRangeMinus = cRangeMinus + cRangeMinus$bsSize
exSumMinus = sum(sgnMerge$signalMinus[extendedRangeMinus])
# extendedRangeMinus$signalToFlankRatio = (bsSumMinus) / (((exSumMinus - bsSumMinus) / 2) + est.offset)
extendedRangeMinus$signalToFlankRatio = (bsSumMinus / exSumMinus)
} else {
extendedRangeMinus = cRangeMinus
}
if (isTRUE(reportScoresPerBindingSite)) {
# report the scores for each binding site
# -> does not allow the usage of the standard plot
df = rbind(as.data.frame(mcols(extendedRangeMinus)), as.data.frame(mcols(extendedRangePlus))) %>%
group_by(bsSize) %>%
mutate(geneWiseFilter = bsFilterStep)
} else {
# combine both and make results dataframe
df = rbind(as.data.frame(mcols(extendedRangeMinus)), as.data.frame(mcols(extendedRangePlus))) %>%
group_by(bsSize) %>%
summarize(signalToFlankRatio = median(signalToFlankRatio), .groups = "keep") %>%
mutate(geneWiseFilter = bsFilterStep)
}
# update chunk counter
counterCurrentIterations <<- counterCurrentIterations + length(bsResolution.steps)
msg = paste0("Estimation at: ", round((counterCurrentIterations/counterTotalIterations)*100 ), "%" )
if(!veryQuiet) message(msg)
return(df)
})
df = do.call("rbind", scoreAllDf)
# ---
# Store results for plotting
if (isTRUE(reportScoresPerBindingSite)){
object@plotData$estimateBsWidth$dataEnhanced = df
} else {
object@plotData$estimateBsWidth$data = df
}
resEstimate = .findFirstMaximum(df, est.minimumStepGain = est.minimumStepGain)
est.bsSize = resEstimate$est.bsSize
est.geneFilter = resEstimate$est.geneFilter
est.option = resEstimate$est.option
if (est.option == "local") {
msg = paste0("No global maximum found with est.maxBsWidth=", est.maxBsWidth, ". Using local maximum instead.\n")
if(!veryQuiet) warning(msg)
}
if (est.option == "fallback") {
msg = paste0("No global or local maximum found with est.maxBsWidth=", est.maxBsWidth, ". Falling back to simple maximum instead.\n")
if(!veryQuiet) warning(msg)
}
if (est.option == "error") {
msg = paste0("No maximum found at all. Raise 'est.minimumStepGain' argument/, change bsResolution and/or geneResolution arguments/ raise est.maxBsWidth and check visually. \n")
stop(msg)
}
# Store estimate type in options (global/ local)
optstr$option = est.option
object@params$estimateBsWidth = optstr
# Store estimated parameters
object@params$bsSize = est.bsSize
object@params$geneFilter = est.geneFilter
# ---
# Store for results
rng = getRanges(object)
resultLine = data.frame(
funName = "estimateBsWidth()", class = "estimate",
nIn = length(rng), nOut = length(rng),
per = paste0(round(length(rng)/ length(rng), digits = 2)*100,"%"),
options = paste0("bsResolution=", bsResolution, ", geneResolution=", geneResolution,
", est.maxBsWidth=", est.maxBsWidth, ", est.maxSites=", est.maxSites, ", est.subsetChromosome=", paste(est.subsetChromosome, collapse = ","))
)
object@results = rbind(object@results, resultLine)
return(object)
}
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.