-
Notifications
You must be signed in to change notification settings - Fork 1
/
4C.R
57 lines (56 loc) · 1.9 KB
/
4C.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
#' Computing virtual 4C profiles
#'
#' From a (m)cool pre-imported in memory, computes a 4C profile
#' using a user-specified `viewpoint`.
#'
#' @param x a `HiCExperiment` object
#' @param viewpoint viewpoint, defined as a GRanges
#' @param use.scores use.scores
#' @return A tibble with the contact frequency of the viewpoint, per bin
#' along the imported genomic range.
#'
#' @import tibble
#' @importFrom scales unit_format
#' @importFrom S4Vectors queryHits
#' @importFrom IRanges IRanges
#' @importFrom GenomicRanges GRanges
#' @importFrom GenomicRanges start
#' @importFrom GenomicRanges end
#' @importFrom GenomicRanges findOverlaps
#' @export
#' @examples
#' library(HiContacts)
#' contacts_yeast <- contacts_yeast()
#' v4C <- virtual4C(contacts_yeast, GenomicRanges::GRanges('II:490001-510000'))
#' v4C
virtual4C <- function(x, viewpoint, use.scores = 'balanced') {
gis <- InteractionSet::interactions(x)
gis$score <- HiCExperiment::scores(x, use.scores)
cm <- Matrix::as.matrix(gi2cm(gis))
cm <- base::as.matrix(cm)
regions <- InteractionSet::regions(gis)
regions_in_viewpoint <- seq_along(regions) %in% S4Vectors::queryHits(
GenomicRanges::findOverlaps(regions, viewpoint)
)
if (sum(regions_in_viewpoint) > 1) {
score <- rowSums(cm[, regions_in_viewpoint], na.rm = TRUE)
}
else {
score <- cm[, regions_in_viewpoint]
}
gr <- GenomicRanges::GRanges(
seqnames = as.vector(GenomicRanges::seqnames(regions)),
IRanges::IRanges(
GenomicRanges::start(regions),
GenomicRanges::end(regions)
),
score = score,
center = GenomicRanges::start(regions) +
(GenomicRanges::end(regions) - GenomicRanges::start(regions))/2,
in_viewpoint = regions_in_viewpoint
)
if (length(as.character(viewpoint)) == 1) {
gr$viewpoint <- as.character(viewpoint)
}
return(gr)
}