diff --git a/.Rbuildignore b/.Rbuildignore index f702134..4186bcf 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -6,3 +6,4 @@ ^docs$ ^pkgdown$ ^\.github$ +^.*\.BiocCheck$ diff --git a/.github/workflows/pkgdown.yaml b/.github/workflows/pkgdown.yaml index ed7650c..d69f3b2 100644 --- a/.github/workflows/pkgdown.yaml +++ b/.github/workflows/pkgdown.yaml @@ -30,6 +30,12 @@ jobs: with: use-public-rspm: true + - name: Install loopcityData from GitHub + run: | + if (!requireNamespace("remotes", quietly = TRUE)) install.packages("remotes") + remotes::install_github("jpflores-13/loopcityData") + shell: Rscript {0} + - uses: r-lib/actions/setup-r-dependencies@v2 with: extra-packages: any::pkgdown, local::. diff --git a/.gitignore b/.gitignore index a69277f..1e5132c 100644 --- a/.gitignore +++ b/.gitignore @@ -9,3 +9,12 @@ outputs/* inst/doc data/processed/ vignettes/*.png +*.BiocCheck/ +*.tar.gz +*.Rcheck/ +.batch_logs/ + +# Reviewer-response materials (BioC resubmission) — kept locally, not in the package +inst/scripts/tad_null_model_sparsity.R +figures/tad_null_model_sparsity.pdf +inst/scripts/call_tads_arrowhead.sbatch diff --git a/DESCRIPTION b/DESCRIPTION index dee0c2b..e747eb4 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -2,9 +2,13 @@ Package: loopcity Type: Package Title: Call Chromatin Communities From Hi-C Data Version: 0.1.0 -Authors@R: person("Sarah", "Parker", role = c("aut", "cre"), - email = "smp3800@gmail.com", - comment = c(ORCID = "0000-0002-2700-3979")) +Authors@R: c( + person("Sarah", "Parker", role = c("aut", "cre"), + email = "smp3800@gmail.com", + comment = c(ORCID = "0000-0002-2700-3979")), + person("Howard Hughes Medical Institute", role = "fnd"), + person("National Institutes of Health", role = "fnd"), + person("National Science Foundation", role = "fnd")) Description: Tools for assigning and visualizing chromatin community interactions from Hi-C data. License: AGPL (>= 3) Encoding: UTF-8 @@ -34,7 +38,7 @@ Imports: stats, S4Vectors Depends: - R (>= 4.1.0) + R (>= 4.6.0) biocViews: Software, HiC, Sequencing, Epigenetics, Network, Clustering, Visualization BugReports: https://github.com/sarmapar/loopcity/issues @@ -42,4 +46,6 @@ URL: https://sarmapar.com/loopcity/, https://github.com/sarmapar/loopcity VignetteBuilder: knitr +Config/Needs/website: github::jpflores-13/loopcityData Config/roxygen2/version: 8.0.0 +RoxygenNote: 7.3.2 diff --git a/NEWS.md b/NEWS.md new file mode 100644 index 0000000..27e6bdd --- /dev/null +++ b/NEWS.md @@ -0,0 +1,4 @@ +# loopcity 0.1.0 + +* Initial release. Tools for assigning and visualizing chromatin community + interactions from Hi-C data. diff --git a/R/assignCommunities.R b/R/assignCommunities.R index 8e1d95e..42ac1af 100644 --- a/R/assignCommunities.R +++ b/R/assignCommunities.R @@ -8,7 +8,7 @@ #' "walktrap", "leiden", "infomap", "label_prop", or "edge_betweenness" #' @param leidenResolution numeric between 0 and 1, higher resolutions lead to #' more smaller communities, while lower resolutions lead to fewer larger -#' communities. See `resolution_parameter` in `igraph::cluster_leiden` +#' communities. See `resolution` in `igraph::cluster_leiden` #' @param pruneUnder numeric, added loops with scores less than this value will #' be removed #' @@ -21,138 +21,179 @@ #' @export #' #' @examples -#' hicFile <- "inst/extdata/GM12878_chr22.hic" -#' -#' mergedLoops <- mergeAnchors(GM12878_10KbLoops, 1) -#' connections <- connectLoopAnchors(mergedLoops, 1e6) -#' scores <- scoreInteractions(connections, hicFile, mergedLoops) -#' -#' assignCommunities(interactions(scores)) +#' data(GM12878_10KbLoopsChr21) +#' loops <- GM12878_10KbLoopsChr21 +#' loops$score <- loops$value +#' assignCommunities(loops, pruneUnder = 0) assignCommunities <- function(loops, scores, clusterType = "leiden", leidenResolution = 0.1, pruneUnder) { - ## Suppress NSE notes in R CMD check anchorCommunity <- chr <- NULL -## Parameter checking ---------------------------------------------------------- + ## Parameter checking ---------------------------------------- call_args <- match.call() ## Check if loops is a GInteractions object, with helpful message for ## converting data frame to GInteractions - if(!is(loops,"GInteractions")){ - if(is(loops,"data.frame") | is(loops,"data.table") | - is(loops,"DataFrame")){ + if (!is(loops, "GInteractions")) { + if (is(loops, "data.frame") | is(loops, "data.table") | + is(loops, "DataFrame")) { rlang::abort(c("`loops` must be a `GInteractions` object", - x=glue("`class(loops)` is {class(loops)}"), - i=glue("See `mariner::as_ginteractions()` to convert", - " a data-frame-like object to", - " a GInteractions object"))) + x = glue("`class(loops)` is {class(loops)}"), + i = glue( + "See `mariner::as_ginteractions()` to convert", + " a data-frame-like object to", + " a GInteractions object" + ) + )) } else { rlang::abort(c("`loops` must be a `GInteractions` object", - x=glue("`class(loops)` is {class(loops)}"))) + x = glue("`class(loops)` is {class(loops)}") + )) } } ## Set scores to loops$score or abort if missing or wrong type/length - if(missing(scores)){ + if (missing(scores)) { scores <- loops$score - if(is.null(scores)){ - rlang::abort(c(glue("argument `scores` is missing, and `loops` ", - "does not contain a column named `score`"))) + if (is.null(scores)) { + rlang::abort(c(glue( + "argument `scores` is missing, and `loops` ", + "does not contain a column named `score`" + ))) } - } else if(!is(scores, "numeric")){ + } else if (!is(scores, "numeric")) { rlang::abort(c("`scores` must be a numeric vector", - x = glue("`class(scores)` is {class(scores)}"))) - } else if(length(scores) != length(loops)){ + x = glue("`class(scores)` is {class(scores)}") + )) + } else if (length(scores) != length(loops)) { rlang::abort(c("`scores` must be the same length as `loops`", - x = glue("`length(scores)` is {length(scores)} and ", - "`length(loops)` is {length(loops)}"))) + x = glue( + "`length(scores)` is {length(scores)} and ", + "`length(loops)` is {length(loops)}" + ) + )) } ## Check that clusterType is valid - possibleClusterTypes <- c("fast_greedy","walktrap", "leiden", "infomap", - "label_prop", "edge_betweenness", - "connected_components") - if(!clusterType %in% possibleClusterTypes){ - rlang::abort(c(glue("`clusterType`=\"{clusterType}\" is not a valid", - " clustering algorithm option"), - i=glue("`clusterType` must be one of \"", - glue_collapse({possibleClusterTypes}, sep ="\", \""), - "\""))) + possibleClusterTypes <- c( + "fast_greedy", "walktrap", "leiden", "infomap", + "label_prop", "edge_betweenness", + "connected_components" + ) + if (!clusterType %in% possibleClusterTypes) { + rlang::abort(c( + glue( + "`clusterType`=\"{clusterType}\" is not a valid", + " clustering algorithm option" + ), + i = glue( + "`clusterType` must be one of \"", + glue_collapse( + { + possibleClusterTypes + }, + sep = "\", \"" + ), + "\"" + ) + )) } ## Check that leiden resolution is a numeric - if(!is(leidenResolution, "numeric")){ + if (!is(leidenResolution, "numeric")) { rlang::abort(c("`leidenResolution` must be numeric", - x=glue("`class(leidenResolution)` is ", - "{class(leidenResolution)}"))) - } else if(length(leidenResolution) > 1){ - rlang::warn(glue("More than one `leidenResolution` provided, only", - " the first will be used")) + x = glue( + "`class(leidenResolution)` is ", + "{class(leidenResolution)}" + ) + )) + } else if (length(leidenResolution) > 1) { + rlang::warn(glue( + "More than one `leidenResolution` provided, only", + " the first will be used" + )) leidenResolution <- leidenResolution[1] } ## Give a warning if leidenResolution is provided but leiden is not ## used as clustering algorithm - if(clusterType != "leiden"){ - if("leidenResolution" %in% names(call_args)){ - rlang::warn(glue("`leidenResolution` is provided but will not be", - " used because `clusterType` is not \"leiden\"")) + if (clusterType != "leiden") { + if ("leidenResolution" %in% names(call_args)) { + rlang::warn(glue( + "`leidenResolution` is provided but will not be", + " used because `clusterType` is not \"leiden\"" + )) } } - #### Filter Scores ------------------------------------------------------------- + #### Filter Scores ------------------------------------------------- ## calculate pruneUnder value if both original and added loops are present ## in "source" column - if(missing(pruneUnder)){ - if(!is.null(loops$source)){ - tempLoops <- data.frame(source = loops$source, - score = scores) + if (missing(pruneUnder)) { + if (!is.null(loops$source)) { + tempLoops <- data.frame( + source = loops$source, + score = scores + ) subset <- split(tempLoops, tempLoops$source) originalLoops <- subset$original newLoops <- subset$added - if(length(originalLoops) == 0 & length(newLoops == 0)){ + if (length(originalLoops) == 0 & length(newLoops == 0)) { ## error if no value provided and source labels are unexpected - rlang::abort(c(glue("`pruneUnder` could not be calculated, ", - "expected \"original\" and \"added\" ", - "in `source` column"), - i=glue("Use a `pruneUnder` value of 0 ", - "to keep all loops."))) + rlang::abort(c( + glue( + "`pruneUnder` could not be calculated, ", + "expected \"original\" and \"added\" ", + "in `source` column" + ), + i = glue( + "Use a `pruneUnder` value of 0 ", + "to keep all loops." + ) + )) } ## set pruneUnder to the median score of original loops pruneUnder <- median(originalLoops$score) - } else { ## error if no value provided and no source column rlang::abort(c("`pruneUnder` must be provided since there is no `source` column in `loops`", - i=glue("Use a `pruneUnder` value of 0 to keep all loops."))) + i = glue("Use a `pruneUnder` value of 0 to keep all loops.") + )) } } -#### Build network ------------------------------------------------------------- + #### Build network ------------------------------------------------- ## Remove scores from any added loops with a score less than pruneUnder - if(!is.null(loops$source)){ + if (!is.null(loops$source)) { scores[scores < pruneUnder & - loops$source == "added"] <- NA + loops$source == "added"] <- NA } else { scores[scores < pruneUnder] <- NA } message(glue("Pruning added loops with a score less than {pruneUnder}")) + ## Ensure loops$score reflects the (possibly pruned) scores so that + ## .compareNodeToComm can compute medians even when scores were passed + ## separately rather than as a column in loops + loops$score <- scores + ## Set nodes to anchors and edges to scores - anc <- InteractionSet::anchors(loops, type="both", id=TRUE) + anc <- InteractionSet::anchors(loops, type = "both", id = TRUE) - relations <- data.frame(from = anc$first, - to = anc$second, - weights = scores, - chr = mariner::seqnames1(loops)) + relations <- data.frame( + from = anc$first, + to = anc$second, + weights = scores, + chr = mariner::seqnames1(loops) + ) relations[is.na(relations)] <- 0 @@ -169,71 +210,82 @@ assignCommunities <- function(loops, GenomicRanges::seqnames() |> unique() - for(chrom in chroms){ + for (chrom in chroms) { # filter to loops in chromosome with scores above pruneUnder - relations_chrom <- dplyr::filter(relations, - chr == chrom) + relations_chrom <- dplyr::filter( + relations, + chr == chrom + ) g <- igraph::graph_from_data_frame(relations_chrom, - directed=FALSE) + directed = FALSE + ) - #### Assign anchors to communities------------------------------------------ + #### Assign anchors to communities ------------------------- # transform weights to log2 - relations_chrom$weights = log2(relations_chrom$weights) + relations_chrom$weights <- log2(relations_chrom$weights) # remove edges with scores 0 or less g_filter <- igraph::delete_edges(g, which(relations_chrom$weights <= 0)) relations_filter <- dplyr::filter(relations_chrom, weights > 0) ## generate communities using chosen algorithm - if(clusterType == "edge_betweenness"){ + if (clusterType == "edge_betweenness") { maxScore <- max(relations_filter$weights) - flippedWeights <- 1 - ((relations_filter$weights)/maxScore) - flippedWeights[which(flippedWeights==0)] <- 1e-10 + flippedWeights <- 1 - ((relations_filter$weights) / maxScore) + flippedWeights[which(flippedWeights == 0)] <- 1e-10 } between_g <- switch(clusterType, - fast_greedy = - igraph::cluster_fast_greedy(g_filter, - weights = - relations_filter$weights), - walktrap = - igraph::cluster_walktrap(g_filter, - weights = - relations_filter$weights), - leiden = - igraph::cluster_leiden(g_filter, - weights = relations_filter$weights, - objective_function = "CPM", - resolution_parameter = - leidenResolution, - n_iterations = 5), - infomap = - igraph::cluster_infomap(g_filter, - e.weights = - relations_filter$weights), - label_prop = - igraph::cluster_label_prop(g_filter, - weights = - relations_filter$weights), - edge_betweenness = - igraph::cluster_edge_betweenness(g_filter, - weights = flippedWeights, - directed = FALSE), - connected_components = - igraph::components(g_filter) + fast_greedy = + igraph::cluster_fast_greedy(g_filter, + weights = + relations_filter$weights + ), + walktrap = + igraph::cluster_walktrap(g_filter, + weights = + relations_filter$weights + ), + leiden = + igraph::cluster_leiden(g_filter, + weights = relations_filter$weights, + objective_function = "CPM", + resolution = leidenResolution, + n_iterations = 5 + ), + infomap = + igraph::cluster_infomap(g_filter, + e.weights = + relations_filter$weights + ), + label_prop = + igraph::cluster_label_prop(g_filter, + weights = + relations_filter$weights + ), + edge_betweenness = + igraph::cluster_edge_betweenness(g_filter, + weights = flippedWeights, + directed = FALSE + ), + connected_components = + igraph::components(g_filter) ) ## put community numbers in genomic order anchorCommunities <- - match(as.numeric(igraph::membership(between_g)), - unique(as.numeric(igraph::membership(between_g)))) + match( + as.numeric(igraph::membership(between_g)), + unique(as.numeric(igraph::membership(between_g))) + ) ## assign loop anchors to communities # add totalCommunities to assigned communities to ensure unique # community numbers InteractionSet::regions(loops)$anchorCommunity[ - as.numeric(names(igraph::membership(between_g)))] <- + as.numeric(names(igraph::membership(between_g))) + ] <- anchorCommunities + totalCommunities ## add number of new unique communities to total communities @@ -241,7 +293,7 @@ assignCommunities <- function(loops, length(unique(as.numeric(igraph::membership(between_g)))) } -#### adjust border nodes-------------------------------------------------------- + #### adjust border nodes ------------------------------------------- ## check if nodes on borders of communities fit well in neighboring ## communities based on scores anchors <- regions(loops) |> as.data.frame() @@ -255,12 +307,12 @@ assignCommunities <- function(loops, dplyr::slice_tail(n = 1) ## TODO replace for loop with map - communityNums <- unique(loops$loopCommunity) |> + communityNums <- unique(InteractionSet::regions(loops)$anchorCommunity) |> na.omit() |> as.numeric() movement <- GenomicRanges::GRanges() - for(comm in communityNums){ + for (comm in communityNums) { nodeLeft <- bordersLeft |> dplyr::filter(anchorCommunity == comm) |> plyranges::as_granges() @@ -269,38 +321,44 @@ assignCommunities <- function(loops, dplyr::filter(anchorCommunity == comm) |> plyranges::as_granges() - ln <- .compareNodeToComm(loops, nodeLeft, TRUE) - rn <- .compareNodeToComm(loops, nodeRight, FALSE) + ln <- .compareNodeToComm(loops, nodeLeft, TRUE, comm) + rn <- .compareNodeToComm(loops, nodeRight, FALSE, comm) movement <- c(movement, ln, rn) } ## filter to only anchors that moved - if(length(movement) > 0){ - movement <- movement[which(vapply(movement$anchorCommunity, - function(x) length(x) > 1, - logical(1)))] + if (length(movement) > 0) { + movement <- movement[which(vapply( + movement$anchorCommunity, + function(x) length(x) > 1, + logical(1) + ))] } ## reassign anchor communities - overlaps <- InteractionSet::findOverlaps(movement, - InteractionSet::regions(loops)) + overlaps <- InteractionSet::findOverlaps( + movement, + InteractionSet::regions(loops) + ) InteractionSet::regions(loops)[S4Vectors::subjectHits(overlaps)] <- movement[S4Vectors::queryHits(overlaps)] ## set loop community to any communities shared by both anchors anchor1com <- InteractionSet::regions(loops)[ - InteractionSet::anchorIds(loops)$first]$anchorCommunity + InteractionSet::anchorIds(loops)$first + ]$anchorCommunity anchor2com <- InteractionSet::regions(loops)[ - InteractionSet::anchorIds(loops)$second]$anchorCommunity + InteractionSet::anchorIds(loops)$second + ]$anchorCommunity loops$loopCommunity <- - mapply(function(x, y) intersect(x,y), anchor1com, anchor2com) + mapply(function(x, y) intersect(x, y), anchor1com, anchor2com) ## remove added loops under the pruning value - if(!is.null(loops$source)){ + if (!is.null(loops$source)) { loops <- loops[which(loops$source == "original" | - as.logical(scores >= pruneUnder))] + as.logical(scores >= pruneUnder))] } else { loops <- loops[which(as.logical(scores >= pruneUnder))] } diff --git a/R/connectLoopAnchors.R b/R/connectLoopAnchors.R index 9ea8d99..9ec638c 100644 --- a/R/connectLoopAnchors.R +++ b/R/connectLoopAnchors.R @@ -1,17 +1,19 @@ #' connectLoopAnchors #' #' @param loops GInteractions object of loops -#' @param overlapDist an integer representing the farthest distance in base pairs -#' to make a connection +#' @param overlapDist an integer representing the farthest distance in +#' base pairs to make a connection #' #' @return GInteractions object #' @export #' #' @examples -#' mergedLoops <- mergeAnchors(GM12878_10KbLoops, 1) -#' -#' connectLoopAnchors(mergedLoops, 1e6) -connectLoopAnchors <- function(loops, overlapDist){ +#' if (requireNamespace("loopcityData", quietly = TRUE)) { +#' data(GM12878_10KbLoops, package = "loopcityData") +#' mergedLoops <- mergeAnchors(GM12878_10KbLoops, 1) +#' connectLoopAnchors(mergedLoops, 1e6) +#' } +connectLoopAnchors <- function(loops, overlapDist) { ## Suppress NSE notes in R CMD check queryHits <- subjectHits <- NULL @@ -31,7 +33,8 @@ connectLoopAnchors <- function(loops, overlapDist){ ## Remove duplicates newInteractions <- InteractionSet::GInteractions( anchor1 = anchors[overlaps$queryHits], - anchor2 = anchors[overlaps$subjectHits]) |> + anchor2 = anchors[overlaps$subjectHits] + ) |> InteractionSet::swapAnchors(mode = "order") |> unique() @@ -46,5 +49,3 @@ connectLoopAnchors <- function(loops, overlapDist){ return(newInteractions) } - - diff --git a/R/data.R b/R/data.R index e69de29..e570390 100644 --- a/R/data.R +++ b/R/data.R @@ -0,0 +1,45 @@ +#' GM12878 10Kb Hi-C Loop Calls (Chr22) +#' +#' A \code{GInteractions} object of 10kb resolution Hi-C loop calls for +#' chromosome 22 from the GM12878 human B-lymphoblastoid cell line. +#' +#' @format A \code{GInteractions} object with 175 interactions and 9 +#' metadata columns: +#' \describe{ +#' \item{color}{RGB color string for visualization} +#' \item{APScoreAvg}{Average anchor pair score} +#' \item{ProbabilityofEnrichment}{Probability of enrichment at loop site} +#' \item{RegAPScoreAvg}{Regional average anchor pair score} +#' \item{Avg_diffMaxNeihgboor_1}{Average difference from max neighbor, +#' anchor 1} +#' \item{Avg_diffMaxNeihgboor_2}{Average difference from max neighbor, +#' anchor 2} +#' \item{avg}{Average Hi-C contact count} +#' \item{std}{Standard deviation of Hi-C contact count} +#' \item{value}{Loop enrichment score} +#' } +#' @usage data(GM12878_10KbLoops) +"GM12878_10KbLoops" + +#' GM12878 10Kb Hi-C Loop Calls (Chr21) +#' +#' A \code{GInteractions} object of 10kb resolution Hi-C loop calls for +#' chromosome 21 from the GM12878 human B-lymphoblastoid cell line. +#' +#' @format A \code{GInteractions} object with 141 interactions and 9 +#' metadata columns: +#' \describe{ +#' \item{color}{RGB color string for visualization} +#' \item{APScoreAvg}{Average anchor pair score} +#' \item{ProbabilityofEnrichment}{Probability of enrichment at loop site} +#' \item{RegAPScoreAvg}{Regional average anchor pair score} +#' \item{Avg_diffMaxNeihgboor_1}{Average difference from max neighbor, +#' anchor 1} +#' \item{Avg_diffMaxNeihgboor_2}{Average difference from max neighbor, +#' anchor 2} +#' \item{avg}{Average Hi-C contact count} +#' \item{std}{Standard deviation of Hi-C contact count} +#' \item{value}{Loop enrichment score} +#' } +#' @usage data(GM12878_10KbLoopsChr21) +"GM12878_10KbLoopsChr21" diff --git a/R/mergeAnchors.R b/R/mergeAnchors.R index e3db67c..472e56e 100644 --- a/R/mergeAnchors.R +++ b/R/mergeAnchors.R @@ -21,10 +21,11 @@ #' @export #' #' @examples -#' mergeAnchors(loops = GM12878_10KbLoops, pixelOverlap = 1) -#' - -mergeAnchors <- function(loops, pixelOverlap = 1, dropDups = TRUE){ +#' if (requireNamespace("loopcityData", quietly = TRUE)) { +#' data(GM12878_10KbLoops, package = "loopcityData") +#' mergeAnchors(loops = GM12878_10KbLoops, pixelOverlap = 1) +#' } +mergeAnchors <- function(loops, pixelOverlap = 1, dropDups = TRUE) { ## Suppress NSE notes in R CMD check subjectHits <- queryHits <- V1 <- anchorId <- overlapIds <- NULL sameGroup <- group <- count <- NULL @@ -37,27 +38,30 @@ mergeAnchors <- function(loops, pixelOverlap = 1, dropDups = TRUE){ InteractionSet::regions(loops) |> IRanges::ranges() |> IRanges::width() |> - unique()-1 + unique() - 1 - if(length(resolution) > 1){ + if (length(resolution) > 1) { rlang::abort(glue("Interaction ranges must have the same width ", - "to use pixelOverlap.", - "x"=glue("There are {length(resolution)} different widths ", - "in `ranges(regions(loops))`."), - "i"="Consider using bpOverlap instead.")) + "to use pixelOverlap.", + "x" = glue( + "There are {length(resolution)} different widths ", + "in `ranges(regions(loops))`." + ), + "i" = "Consider using bpOverlap instead." + )) } ## Create data table of all overlapping anchor indexes ## (one line per overlap) overlaps <- loops |> InteractionSet::regions() |> - IRanges::findOverlaps(maxgap = pixelOverlap*resolution-2) |> + IRanges::findOverlaps(maxgap = pixelOverlap * resolution - 2) |> data.table::as.data.table() ## Condense overlaps into list of all overlapping anchor indexes ## (one line per unique anchor) ## & give each chain a group number - overlapGroups <- overlaps[, .(list(subjectHits)), by=queryHits] |> + overlapGroups <- overlaps[, .(list(subjectHits)), by = queryHits] |> dplyr::rename(anchorId = queryHits, overlapIds = V1) |> dplyr::rowwise() |> dplyr::mutate(sameGroup = anchorId %in% dplyr::lead(overlapIds)) |> @@ -75,17 +79,23 @@ mergeAnchors <- function(loops, pixelOverlap = 1, dropDups = TRUE){ ## helper function to pick the middle value of a vector ## (or the largest of 2 middle vals) - middle <- function(x) {x[floor(length(x)/2) + 1]} + middle <- function(x) { + x[floor(length(x) / 2) + 1] + } - ## For each group of anchors, collect all anchor IDs with the maximum counts (modes) + ## For each group of anchors, collect all anchor IDs with the + ## maximum counts (modes) ## If there is only one mode, that is the representative anchor ## If there is more than one mode, choose the middle-est one representativeAnchors <- overlapGroups |> dplyr::group_by(group) |> - dplyr::mutate(repAnchorId = {modes <- anchorId[count == max(count)]; - ifelse(length(modes) > 1, - middle(modes), - modes)}) + dplyr::mutate(repAnchorId = { + modes <- anchorId[count == max(count)] + ifelse(length(modes) > 1, + middle(modes), + modes + ) + }) ## Replace original loop anchor IDs with new anchor IDs and create ## new GInteractions object @@ -94,20 +104,22 @@ mergeAnchors <- function(loops, pixelOverlap = 1, dropDups = TRUE){ InteractionSet::regions(loops)[representativeAnchors$repAnchorId] mergedLoops <- InteractionSet::GInteractions( - anchor1 = InteractionSet::anchors(mergedLoops, 'first'), - anchor2 = InteractionSet::anchors(mergedLoops, 'second')) + anchor1 = InteractionSet::anchors(mergedLoops, "first"), + anchor2 = InteractionSet::anchors(mergedLoops, "second") + ) ## Remove duplicates and give a message of how many interactions were ## dropped, if `dropDups` == TRUE - if(dropDups){ - if(length(mergedLoops) != length(unique(mergedLoops))){ - message(glue("Duplicates dropped. {length(mergedLoops)} ", - "interactions reduced to {length(unique(mergedLoops))}", - " interactions.")) + if (dropDups) { + if (length(mergedLoops) != length(unique(mergedLoops))) { + message(glue( + "Duplicates dropped. {length(mergedLoops)} ", + "interactions reduced to {length(unique(mergedLoops))}", + " interactions." + )) } return(unique(mergedLoops)) - } else { return(mergedLoops) } diff --git a/R/plotHicCommunities.R b/R/plotHicCommunities.R index 7d1c4b3..c9cf7ae 100644 --- a/R/plotHicCommunities.R +++ b/R/plotHicCommunities.R @@ -19,16 +19,14 @@ #' @export #' #' @examples +#' data(GM12878_10KbLoops) #' hicFile <- system.file("extdata", "GM12878_chr22.hic", package = "loopcity") -#' -#' \donttest{ #' communities <- GM12878_10KbLoops |> #' mergeAnchors() |> #' connectLoopAnchors(overlapDist = 1e6) |> #' scoreInteractions(hicFile = hicFile) |> #' InteractionSet::interactions() |> #' assignCommunities() -#' #' plotHicCommunities( #' pdfName = tempfile(fileext = ".pdf"), #' communities = communities, @@ -37,104 +35,119 @@ #' starts = 25e6, #' ends = 30e6 #' ) -#' } plotHicCommunities <- function(pdfName, communities, hicFile, norm = "SCALE", - chroms, starts, ends, zmax, colorPalette){ -#### parameter checking -------------------------------------------------------- + chroms, starts, ends, zmax, colorPalette) { + #### parameter checking -------------------------------------------- ## if only one chromosome, replicate for each region - if(length(chroms) == 1){ - chroms = rep(chroms, length(starts)) + if (length(chroms) == 1) { + chroms <- rep(chroms, length(starts)) } - if(length(chroms) != length(starts) | length(chroms) != length(ends)){ + if (length(chroms) != length(starts) | length(chroms) != length(ends)) { abort("`chroms`, `starts`, and `ends` must be the same length") } - #todo check if colorPalette is a color palette - if(missing(colorPalette)){ - colorPalette = loopcityColors + # todo check if colorPalette is a color palette + if (missing(colorPalette)) { + colorPalette <- loopcityColors } -#### plan plotgardener page ---------------------------------------------------- - plotW = 9 - hicH = 2 - labelH = 0.5 - archH = 0.75 - archSpace = 0.25 - border = 0.5 + #### plan plotgardener page ---------------------------------------- + plotW <- 9 + hicH <- 2 + labelH <- 0.5 + archH <- 0.75 + archSpace <- 0.25 + border <- 0.5 - pageHeight = border + hicH + + pageHeight <- border + hicH + archH + archSpace + labelH + border - pageWidth = plotW + 2*border + pageWidth <- plotW + 2 * border pdf(pdfName, width = pageWidth, height = pageHeight) - for(i in seq_along(chroms)){ + for (i in seq_along(chroms)) { # set region to current chr, start, end - chr = chroms[i] - start = starts[i] - end = ends[i] + chr <- chroms[i] + start <- starts[i] + end <- ends[i] # params regions <- InteractionSet::regions(communities) - regions <-regions[ + regions <- regions[ GenomicRanges::seqnames(regions) == chr & GenomicRanges::start(regions) >= start & - GenomicRanges::end(regions) <= end] - - if(missing(zmax)){ - z = plotgardener::readHic(file = hicFile, - chrom = chr, - chromstart = start, - chromend = end, - norm = norm, - resolution = 10e3) - - if(nrow(z) > 0){ - zmax = as.numeric(ceiling(quantile(z[z[,3] >= 10,3])[4])) + 100 - } else{ - zmax = 0 - rlang::abort("`zmax` could not be calculated, please provide a value.") + GenomicRanges::end(regions) <= end + ] + + if (missing(zmax)) { + z <- plotgardener::readHic( + file = hicFile, + chrom = chr, + chromstart = start, + chromend = end, + norm = norm, + resolution = 10e3 + ) + + if (nrow(z) > 0) { + zmax <- as.numeric( + ceiling(quantile(z[z[, 3] >= 10, 3])[4]) + ) + 100 + } else { + zmax <- 0 + rlang::abort(paste0( + "`zmax` could not be calculated, ", + "please provide a value." + )) } } # make plot - plotgardener::pageCreate(width = pageWidth, - height = pageHeight, - showGuides = FALSE) + plotgardener::pageCreate( + width = pageWidth, + height = pageHeight, + showGuides = FALSE + ) # hic plot - hicPlot <- plotgardener::plotHicRectangle(data = hicFile, - chrom = chr, - chromstart = start, - chromend = end, - norm = norm, - resolution = 10e3, - zrange = c(0,zmax), - x = border, - y = border, - width = plotW, - height = hicH) + hicPlot <- plotgardener::plotHicRectangle( + data = hicFile, + chrom = chr, + chromstart = start, + chromend = end, + norm = norm, + resolution = 10e3, + zrange = c(0, zmax), + x = border, + y = border, + width = plotW, + height = hicH + ) # filter for # - loops in region # - scores greater than 0 (also removes NAs) # - loops in a community - windowLoops_i <- GenomicRanges::findOverlaps(regions, - communities) |> + windowLoops_i <- GenomicRanges::findOverlaps( + regions, + communities + ) |> S4Vectors::to() windowLoops <- communities[windowLoops_i] - windowLoops <- windowLoops[which(as.logical(windowLoops$score> 0) & - lengths(windowLoops$loopCommunity)>0)] + windowLoops <- windowLoops[which(as.logical(windowLoops$score > 0) & + lengths(windowLoops$loopCommunity) > 0)] # if `source` column exists, annotate pixels # gray box = added loop, black box = original loop noSource <- FALSE - if(is.null(windowLoops$source)){ - rlang::warn(glue::glue("Cannot differentiate added vs original ", - "loops in annotated pixels or arches,", - " no `source` column in `communities`.")) + if (is.null(windowLoops$source)) { + rlang::warn(glue::glue( + "Cannot differentiate added vs original ", + "loops in annotated pixels or arches,", + " no `source` column in `communities`." + )) tempLoops <- windowLoops S4Vectors::mcols(tempLoops) <- NULL @@ -144,65 +157,73 @@ plotHicCommunities <- function(pdfName, communities, hicFile, norm = "SCALE", } else { # annotating original loops originalLoops <- windowLoops[which(as.logical( - windowLoops$source == "original"))] + windowLoops$source == "original" + ))] S4Vectors::mcols(originalLoops) <- NULL plotgardener::annoPixels(hicPlot, originalLoops) # annotating added loops addedLoops <- windowLoops[which(as.logical( - windowLoops$source == "added"))] + windowLoops$source == "added" + ))] S4Vectors::mcols(addedLoops) <- NULL - if(length(addedLoops) > 0){ + if (length(addedLoops) > 0) { plotgardener::annoPixels(hicPlot, addedLoops, col = "gray") } } # plot genome label and legend for hicplot plotgardener::annoGenomeLabel(hicPlot, - x = border, - y = pageHeight - border - labelH, - scale = "Mb") + x = border, + y = pageHeight - border - labelH, + scale = "Mb" + ) - plotgardener::annoHeatmapLegend(hicPlot, x = border-0.25, y = border, - width = 0.2, height = hicH*(3/5)) + plotgardener::annoHeatmapLegend(hicPlot, + x = border - 0.25, y = border, + width = 0.2, height = hicH * (3 / 5) + ) ## Plot arches, color by community and height based on log2(score) # Map unique values to colors currColors <- colorPalette(length(unique(windowLoops$loopCommunity))) - value_to_color <- setNames(currColors, - unique(windowLoops$loopCommunity)) + value_to_color <- setNames( + currColors, + unique(windowLoops$loopCommunity) + ) # Assign colors based on values windowLoops$color <- value_to_color[ - as.character(windowLoops$loopCommunity)] + as.character(windowLoops$loopCommunity) + ] - #calculate log 2 of scores for plotting arches + # calculate log 2 of scores for plotting arches windowLoops$score_lg <- log2(windowLoops$score) # find max score max_score <- max(windowLoops$score_lg) # plot arches - if(max_score > 0){ - arches <- plotgardener::plotPairsArches(windowLoops, - chrom = chr, - chromstart = start, - chromend = end, - flip = TRUE, - clip = TRUE, - archHeight = "score_lg", - range = c(0, max_score), - x = border, - y = border + hicH, - height = archH, - width = plotW, - fill = windowLoops$color) + if (max_score > 0) { + arches <- plotgardener::plotPairsArches(windowLoops, + chrom = chr, + chromstart = start, + chromend = end, + flip = TRUE, + clip = TRUE, + archHeight = "score_lg", + range = c(0, max_score), + x = border, + y = border + hicH, + height = archH, + width = plotW, + fill = windowLoops$color + ) } } dev.off() message(glue::glue("Plots created in "), pdfName) } - diff --git a/R/scoreInteractions.R b/R/scoreInteractions.R index 73f4afb..dab85ec 100644 --- a/R/scoreInteractions.R +++ b/R/scoreInteractions.R @@ -6,31 +6,34 @@ #' @returns a MatrixSelection object corresponding to the keyword #' @noRd #' @keywords internal -.pickSelection <- function(x = c("bowtie","donut"), n, buffer){ +.pickSelection <- function(x = c("bowtie", "donut"), n, buffer) { x <- match.arg(x) - switch (x, - bowtie = mariner::selectTopLeft(n = n, buffer = buffer) + - mariner::selectBottomRight(n = n, buffer = buffer), - donut = mariner::selectOuter(n = n, buffer = buffer) + switch(x, + bowtie = mariner::selectTopLeft(n = n, buffer = buffer) + + mariner::selectBottomRight(n = n, buffer = buffer), + donut = mariner::selectOuter(n = n, buffer = buffer) ) } ## Generic for scoreInteractions -#' @rdname calcLoopEnrichment +#' @rdname scoreInteractions #' @export -setGeneric("scoreInteractions", - function(x, - hicFile, - fgSize = 0, - bg = "bowtie", - bgGap = 0, - bgSize = 5, - pseudo = 5, - pruneUnder = 0, - truncate = FALSE, - resolution, - norm = "NONE") - standardGeneric("scoreInteractions")) +setGeneric( + "scoreInteractions", + function(x, + hicFile, + fgSize = 0, + bg = "bowtie", + bgGap = 0, + bgSize = 5, + pseudo = 5, + pruneUnder = 0, + truncate = FALSE, + resolution, + norm = "NONE") { + standardGeneric("scoreInteractions") + } +) #' Internal for scoreInteractionsFromMats @@ -38,23 +41,27 @@ setGeneric("scoreInteractions", #' @importFrom mariner calcLoopEnrichment #' @noRd #' @keywords internal -.scoreInteractionsFromMats <- function(x,fgSize,bg,bgGap,bgSize, - pseudo,pruneUnder,truncate){ - +.scoreInteractionsFromMats <- function(x, fgSize, bg, bgGap, bgSize, + pseudo, pruneUnder, truncate) { ## Get buffer from counts matrices - buffer <- (dim(BiocGenerics::counts(x))[1] - 1)/2 + buffer <- (dim(BiocGenerics::counts(x))[1] - 1) / 2 ## Check that buffer is large enough for bg and fg - if(buffer < fgSize + bgGap + bgSize){ + if (buffer < fgSize + bgGap + bgSize) { abort(c("`x` dimensions are too small", - "x"=glue('`dim(counts(x))` must be at least {fgSize+bgGap+bgSize}'), - "i"=glue('`dim(counts(x))` is {dim(counts(x))}'), - "i"="Try using a smaller value for fgSize, bgGap, or bgSize, - or re-run with loop and hi-c files.")) + "x" = glue( + "`dim(counts(x))` must be at least {fgSize+bgGap+bgSize}" + ), + "i" = glue("`dim(counts(x))` is {dim(counts(x))}"), + "i" = "Try using a smaller value for fgSize, bgGap, or bgSize, + or re-run with loop and hi-c files." + )) } ## Select center pixel/foreground - cp_sel <- mariner::selectCenterPixel(mhDist = 1:fgSize, buffer = buffer) + cp_sel <- mariner::selectCenterPixel( + mhDist = seq_len(fgSize), buffer = buffer + ) ## Select background based on bg keyword bg_sel <- .pickSelection(bg, n = bgSize, buffer = buffer) @@ -63,26 +70,32 @@ setGeneric("scoreInteractions", x$pseudocounts <- pseudo ## Function to divide median of fg and bg with pseudocounts - scoringFunction <- function(fg,bg,pseudo=pseudo){ - median((fg + 1 + pseudo)/(bg + 1 + pseudo), na.rm = TRUE) + scoringFunction <- function(fg, bg, pseudo = pseudo) { + median((fg + 1 + pseudo) / (bg + 1 + pseudo), na.rm = TRUE) } ## Calculate enrichment scores - scores <- mariner::calcLoopEnrichment(x = x, fg = cp_sel, bg = bg_sel, - FUN = scoringFunction) + scores <- mariner::calcLoopEnrichment( + x = x, fg = cp_sel, bg = bg_sel, + FUN = scoringFunction + ) ## Prune away scores where the center pixels have ## raw counts less than pruneUnder - centerPx <- mariner::calcLoopEnrichment(x = x, - fg = mariner::selectCenterPixel( - 0, buffer = buffer), - FUN = function(fg,bg) fg) + centerPx <- mariner::calcLoopEnrichment( + x = x, + fg = mariner::selectCenterPixel( + 0, + buffer = buffer + ), + FUN = function(fg, bg) fg + ) scores[centerPx < pruneUnder] <- NA ## truncate high scores to 99th percentile - if(truncate){ - newMax <- quantile(scores,.99) - scores[which(scores>newMax),] <- newMax + if (truncate) { + newMax <- quantile(scores, .99) + scores[which(scores > newMax), ] <- newMax } GenomicRanges::mcols(x)$score <- scores @@ -95,23 +108,23 @@ setGeneric("scoreInteractions", #' @importFrom mariner calcLoopEnrichment #' @noRd #' @keywords internal -.scoreInteractionsFromFiles <- function(x,hicFile,fgSize,bg,bgGap, - bgSize,pseudo,pruneUnder,truncate, - resolution,norm){ - +.scoreInteractionsFromFiles <- function(x, hicFile, fgSize, bg, bgGap, + bgSize, pseudo, pruneUnder, truncate, + resolution, norm) { ## get resolution from loops - if(missing(resolution)){ + if (missing(resolution)) { uniqueWidths <- InteractionSet::regions(x) |> IRanges::width() |> unique() - if(length(uniqueWidths) > 1){ + if (length(uniqueWidths) > 1) { abort(c("Regions in `x` must be the same resolution.", - "x"=glue('`width(regions(x))` include - {length(uniqueWidths)} different widths.'), - "i"=glue("Providing a value for `resolution` will put - regions into equal sized bins."))) + "x" = glue("`width(regions(x))` include + {length(uniqueWidths)} different widths."), + "i" = glue("Providing a value for `resolution` will put + regions into equal sized bins.") + )) } - resolution = uniqueWidths-1 + resolution <- uniqueWidths - 1 } ## pull hi-c contacts @@ -119,11 +132,11 @@ setGeneric("scoreInteractions", buffer <- fgSize + bgGap + bgSize ## Extract matrices around each loop - mats <- mariner::pixelsToMatrices(x=x, buffer=buffer) |> + mats <- mariner::pixelsToMatrices(x = x, buffer = buffer) |> mariner::pullHicMatrices( - files=hicFile, - binSize=resolution, - half="upper", + files = hicFile, + binSize = resolution, + half = "upper", norm = norm ) @@ -131,14 +144,16 @@ setGeneric("scoreInteractions", InteractionSet::interactions(mats) <- x ## Pass in parameters, including optional ones, to the next function - args <- list(x = mats, - fgSize = fgSize, - bg = bg, - bgGap = bgGap, - bgSize = bgSize, - pruneUnder = pruneUnder, - truncate = truncate, - pseudo = pseudo) + args <- list( + x = mats, + fgSize = fgSize, + bg = bg, + bgGap = bgGap, + bgSize = bgSize, + pruneUnder = pruneUnder, + truncate = truncate, + pseudo = pseudo + ) do.call(".scoreInteractionsFromMats", args) } @@ -186,17 +201,19 @@ setGeneric("scoreInteractions", #' @export #' #' @examples -#' hicFile <- "inst/extdata/GM12878_chr22.hic" -#' +#' data(GM12878_10KbLoops) +#' hicFile <- system.file("extdata", "GM12878_chr22.hic", package = "loopcity") #' mergedLoops <- mergeAnchors(GM12878_10KbLoops, 1) #' connections <- connectLoopAnchors(mergedLoops, 1e6) -#' -#' scoreInteractions(connections, hicFile, mergedLoops) +#' scoreInteractions(connections, hicFile) #' setMethod("scoreInteractions", - signature(x="GInteractions", - hicFile="character"), - definition=.scoreInteractionsFromFiles) + signature( + x = "GInteractions", + hicFile = "character" + ), + definition = .scoreInteractionsFromFiles +) #' Calculate loop enrichment over background. @@ -204,7 +221,9 @@ setMethod("scoreInteractions", #' @rdname scoreInteractions #' @export setMethod("scoreInteractions", - signature(x="InteractionArray", - hicFile="missing"), - definition=.scoreInteractionsFromMats) - + signature( + x = "InteractionArray", + hicFile = "missing" + ), + definition = .scoreInteractionsFromMats +) diff --git a/R/utils.R b/R/utils.R index 8ab9591..0f8c88e 100644 --- a/R/utils.R +++ b/R/utils.R @@ -7,11 +7,11 @@ #' #' @noRd #' @keywords internal -.optimizePseudocounts <- function(mats,values,loopCalls,fgSize,bg,bgGap, bgSize, - pseudo,pruneUnder,truncate){ - +.optimizePseudocounts <- function(mats, values, loopCalls, fgSize, + bg, bgGap, bgSize, pseudo, + pruneUnder, truncate) { pseudoAUCs <- data.table::data.table(auc = numeric(), psc = numeric()) - for(pseudo in values){ + for (pseudo in values) { loopScores <- scoreInteractions( pseudo = pseudo, x = mats, @@ -20,7 +20,8 @@ bgGap = bgGap, bgSize = bgSize, pruneUnder = pruneUnder, - truncate = truncate) + truncate = truncate + ) ## pull out interactions from InteractionArray inters <- InteractionSet::interactions(loopScores) @@ -30,13 +31,18 @@ inters$called_loop[inters %in% loopCalls] <- TRUE ## calculate area under the roc curve - roc <- pROC::roc(predictor = as.numeric(inters$score), - response = inters$called_loop) + roc <- pROC::roc( + predictor = as.numeric(inters$score), + response = inters$called_loop + ) auc <- pROC::auc(roc) - pseudoAUCs <- dplyr::bind_rows(pseudoAUCs, - data.table::data.table( - auc = as.numeric(auc), - psc = pseudo)) + pseudoAUCs <- dplyr::bind_rows( + pseudoAUCs, + data.table::data.table( + auc = as.numeric(auc), + psc = pseudo + ) + ) } ## return pseudocounts with maximum AUC @@ -44,7 +50,6 @@ } - #' Helper function to determine if the node on the border of a community should #' be assigned to both its current community and the neighboring community #' @param loops GInteractions object containing all loops of interest, must also @@ -57,64 +62,80 @@ #' #' @noRd #' @keywords internal -.compareNodeToComm <- function(loops, borderNode, left){ - ## Suppress NSE notes in R CMD check - queryHits <- loopCommunity <- comm <- score <- NULL +.compareNodeToComm <- function(loops, borderNode, left, comm) { ## Find loops that contain the borderNode (overlaps 1) and anchors in ## the neighboring community (overlaps 2) - if(left){ + if (left) { overlaps1 <- InteractionSet::findOverlaps( - InteractionSet::anchors(loops,"second"), borderNode) + InteractionSet::anchors(loops, "second"), borderNode + ) overlaps2 <- InteractionSet::findOverlaps( - InteractionSet::anchors(loops,"first"), + InteractionSet::anchors(loops, "first"), InteractionSet::regions(loops)[which( - InteractionSet::regions(loops)$anchorCommunity == - borderNode$anchorCommunity-1)]) + InteractionSet::regions(loops)$anchorCommunity == + borderNode$anchorCommunity - 1 + )] + ) } else { overlaps1 <- InteractionSet::findOverlaps( - InteractionSet::anchors(loops,"first"), borderNode) + InteractionSet::anchors(loops, "first"), borderNode + ) overlaps2 <- InteractionSet::findOverlaps( - InteractionSet::anchors(loops,"second"), + InteractionSet::anchors(loops, "second"), InteractionSet::regions(loops)[which( - InteractionSet::regions(loops)$anchorCommunity == - borderNode$anchorCommunity+1)]) + InteractionSet::regions(loops)$anchorCommunity == + borderNode$anchorCommunity + 1 + )] + ) } # Keep only interactions between border node and neighboring community nodes interactions_to_keep <- overlaps1[ - S4Vectors::queryHits(overlaps1) %in% S4Vectors::queryHits(overlaps2)] + S4Vectors::queryHits(overlaps1) %in% S4Vectors::queryHits(overlaps2) + ] ## Compare the median score of all interactions within the neighboring ## community (comm_med) to all interactions between the border node and ## anchors in the neighboring community (node_med) - filtered_scores <- loops$score[queryHits(interactions_to_keep)] - node_med <- filtered_scores |> median() - - comm_med <- loops |> - as.data.frame() |> - dplyr::filter(loopCommunity == comm) |> - dplyr::pull(score) |> - median() + filtered_scores <- loops$score[S4Vectors::queryHits(interactions_to_keep)] + node_med <- median(filtered_scores, na.rm = TRUE) + + ## Derive comm_med from anchorCommunity on regions — loopCommunity is + ## not yet set on loops at the time this function is called + anc1_comm <- InteractionSet::regions(loops)[ + InteractionSet::anchorIds(loops)$first + ]$anchorCommunity + anc2_comm <- InteractionSet::regions(loops)[ + InteractionSet::anchorIds(loops)$second + ]$anchorCommunity + in_comm <- which(anc1_comm == comm & anc2_comm == comm) + comm_med <- median(loops$score[in_comm], na.rm = TRUE) ## If node_med >= comm_med, assign the border node to both communities - if(!is.na(node_med) & node_med >= comm_med){ - if(left){ - borderNode$anchorCommunity <- list(c(comm-1, - borderNode$anchorCommunity)) + if (isTRUE(node_med >= comm_med)) { + if (left) { + borderNode$anchorCommunity <- list(c( + comm - 1, + borderNode$anchorCommunity + )) } else { - borderNode$anchorCommunity <- list(c(borderNode$anchorCommunity, - comm+1)) + borderNode$anchorCommunity <- list(c( + borderNode$anchorCommunity, + comm + 1 + )) } } return(borderNode) } -##todo document -loopcityColors <- function(n){ - baseColors <- c("chartreuse3", "deepskyblue3", "darkorange", - "darkorchid2", "deeppink3") +## todo document +loopcityColors <- function(n) { + baseColors <- c( + "chartreuse3", "deepskyblue3", "darkorange", + "darkorchid2", "deeppink3" + ) return(rep(baseColors, length.out = n)) } diff --git a/R/zzz.R b/R/zzz.R new file mode 100644 index 0000000..2657e21 --- /dev/null +++ b/R/zzz.R @@ -0,0 +1,13 @@ +.onLoad <- function(libname, pkgname) { + # Seqinfo 1.0.0 was split from GenomeInfoDb, creating two separate seqinfo + # generics. InteractionSet registers its GInteractions method against + # Seqinfo::seqinfo, but mariner imports seqinfo from GenomeInfoDb (a different + # generic object). Register a bridge method against GenomeInfoDb::seqinfo so + # that mariner::pullHicMatrices can resolve chromosome information correctly. + gid_seqinfo <- GenomeInfoDb::seqinfo + if (!existsMethod(gid_seqinfo, "GInteractions")) { + setMethod(gid_seqinfo, "GInteractions", function(x) { + GenomeInfoDb::seqinfo(InteractionSet::regions(x)) + }) + } +} diff --git a/README.Rmd b/README.Rmd index 26c6f09..fd6fc0f 100644 --- a/README.Rmd +++ b/README.Rmd @@ -25,8 +25,9 @@ loopcity is currently being submitted to [Bioconductor](https://www.bioconductor You can install the development version of loopcity from [GitHub](https://github.com/) with: ```{r installation, eval = F} -if (!requireNamespace("remotes", quietly = TRUE)) +if (!requireNamespace("remotes", quietly = TRUE)) { install.packages("remotes") +} remotes::install_github("sarmapar/loopcity") ``` diff --git a/README.md b/README.md index 2663e2d..573b770 100644 --- a/README.md +++ b/README.md @@ -51,7 +51,7 @@ remotes::install_github("sarmapar/loopcity") ## Overview -![](man/figures/Fig1.png) **Loopcity detects chromatin communities from +![Overview of the loopcity pipeline showing chromatin community detection from Hi-C data](man/figures/Fig1.png) **Loopcity detects chromatin communities from Hi-C data.** (a) Schematic overview of the main loopcity functions. (From top to bottom) Nearby anchors are merged to create network nodes, edges are added to connect all nearby nodes, edge weights are added diff --git a/docs/reference/figures/Fig1.png b/docs/reference/figures/Fig1.png index 82f96b5..7334225 100644 Binary files a/docs/reference/figures/Fig1.png and b/docs/reference/figures/Fig1.png differ diff --git a/inst/scripts/border_fix_genomewide_diff.R b/inst/scripts/border_fix_genomewide_diff.R new file mode 100644 index 0000000..700ab09 --- /dev/null +++ b/inst/scripts/border_fix_genomewide_diff.R @@ -0,0 +1,203 @@ +# Border-adjustment bug fix: genome-wide before/after diff ---------------- +# Author: JP Flores +# Date: 2026-07-13 +# Project: loopcity +# Description: Quantifies the real-data impact of the assignCommunities() +# border-adjustment fix (commits fdc1731, 527eb78, ce26123), +# addressing PI question "are we sure this doesn't change the +# results?" The chr22 toy-data test (188 loops) showed the fix +# is not a no-op (1 previously-orphaned loop correctly rescued), +# but that dataset is too small to know the real magnitude. +# +# Isolates ONLY the border-adjustment logic: mergeAnchors, +# connectLoopAnchors, and scoreInteractions are run once via the +# current package source (identical for both arms), and the +# RNG seed is pinned immediately before each assignCommunities() +# call so Leiden clustering draws the same random sequence. +# Only .compareNodeToComm / assignCommunities differ between +# the "old" (pre-fix, commit b252df2) and "new" (current) runs. +# +# Also flags whether any changed loop falls inside the +# chr1:110-115Mb region used for Figure 1b/c in +# survey_communities.R, since that's the one figure already +# built from this pipeline's output. +# +# Input: data/processed/hic/GSE214123_MEGA_K562_WT_0_inter.hic +# data/processed/hic/5kbLoops_0.txt (K562 loop calls) +# Output: data/processed/border_fix_diff.csv (all loops whose +# loopCommunity assignment changed) +# prints summary + manuscript sentence to console +# --------------------------------------------------------------------------- + + +# Parameters ---------------------------------------------------------------- + +## Path to the loopcity package SOURCE checkout (a git clone of +## jpflores-13/loopcity with full history) — used only to pull the pre-fix +## versions of assignCommunities.R/utils.R via `git show`. This is likely a +## different path than `project_dir` below if data and package source live +## in separate directories on the cluster. +loopcity_repo_dir <- "/work/users/j/p/jpflores/projects/loopcity" + +project_dir <- "/work/users/j/p/jpflores/projects/loopcity" +data_dir <- file.path(project_dir, "data") + +hic_file <- file.path(data_dir, "processed/hic/GSE214123_MEGA_K562_WT_0_inter.hic") +loops_file <- file.path(data_dir, "processed/hic/5kbLoops_0.txt") + +hic_resolution <- 10e3 +hic_norm <- "NONE" +pixel_overlap <- 1 +overlap_dist <- 1e6 +leiden_res <- 0.5 ## matches survey_communities.R / benchmark_timing.R + +pre_fix_commit <- "b252df2" ## last commit before the border-adjustment fixes + +## Figure 1b/c region (survey_communities.R) — flagged separately below +survey_chrom <- "1" +survey_chromstart <- 110e6 +survey_chromend <- 115e6 + + +# Libraries ------------------------------------------------------------- + +suppressPackageStartupMessages({ + devtools::load_all(loopcity_repo_dir, quiet = TRUE) + library(data.table) + library(GenomeInfoDb) + library(GenomicRanges) + library(glue) + library(InteractionSet) + library(mariner) +}) + + +# Run shared pipeline steps ONCE (identical for both arms) ----------------- + +message("Running mergeAnchors -> connectLoopAnchors -> scoreInteractions ...") + +merged <- data.table::fread(loops_file) |> + mariner::assignToBins(binSize = 10e3) |> + mergeAnchors(pixelOverlap = pixel_overlap) +GenomeInfoDb::seqlevelsStyle(merged) <- "ENSEMBL" + +connected <- connectLoopAnchors(merged, overlapDist = overlap_dist) +scored <- scoreInteractions(connected, + hicFile = hic_file, + resolution = hic_resolution, + norm = hic_norm) |> + InteractionSet::interactions() + + +# Run 1: current (fixed) assignCommunities ----------------------------------- + +fixed_assignCommunities <- assignCommunities + +set.seed(1) +new_result <- fixed_assignCommunities(scored, leidenResolution = leiden_res) +message(glue("NEW (fixed) run: {length(new_result)} loops")) + + +# Run 2: pre-fix assignCommunities, pulled from git and sourced to shadow ---- + +tmp_dir <- tempdir() +old_utils_path <- file.path(tmp_dir, "old_utils.R") +old_assign_path <- file.path(tmp_dir, "old_assignCommunities.R") + +get_old_file <- function(relpath, out_path) { + ok <- system2("git", + args = c("-C", shQuote(loopcity_repo_dir), "show", + glue("{pre_fix_commit}:{relpath}")), + stdout = out_path) + if (!isTRUE(ok == 0) || file.size(out_path) == 0) { + rlang::abort(glue( + "Could not retrieve {relpath} at {pre_fix_commit} from ", + "{loopcity_repo_dir}. Make sure that directory is a git checkout ", + "with this commit (try `git fetch origin` there)." + )) + } +} +get_old_file("R/utils.R", old_utils_path) +get_old_file("R/assignCommunities.R", old_assign_path) + +source(old_utils_path) +source(old_assign_path) +stopifnot(!identical(environment(assignCommunities), + environment(fixed_assignCommunities))) + +set.seed(1) +old_result <- assignCommunities(scored, leidenResolution = leiden_res) +message(glue("OLD (pre-fix) run: {length(old_result)} loops")) + + +# Compare --------------------------------------------------------------- + +mk_key <- function(x) paste( + GenomeInfoDb::seqnames(InteractionSet::anchors(x, "first")), + start(InteractionSet::anchors(x, "first")), + start(InteractionSet::anchors(x, "second")) +) + +old_df <- data.frame( + key = mk_key(old_result), + chr = as.character(GenomeInfoDb::seqnames(InteractionSet::anchors(old_result, "first"))), + start1 = start(InteractionSet::anchors(old_result, "first")), + start2 = start(InteractionSet::anchors(old_result, "second")), + loopCommunity_old = I(lapply(old_result$loopCommunity, sort)) +) +new_df <- data.frame( + key = mk_key(new_result), + loopCommunity_new = I(lapply(new_result$loopCommunity, sort)) +) + +merged_df <- merge(old_df, new_df, by = "key") +changed <- mapply(function(a, b) !identical(a, b), + merged_df$loopCommunity_old, merged_df$loopCommunity_new) + +diff_df <- merged_df[changed, c("chr", "start1", "start2", + "loopCommunity_old", "loopCommunity_new")] +diff_df$loopCommunity_old <- vapply(diff_df$loopCommunity_old, + paste, character(1), collapse = ",") +diff_df$loopCommunity_new <- vapply(diff_df$loopCommunity_new, + paste, character(1), collapse = ",") +diff_df$rescued <- diff_df$loopCommunity_old == "" + +## Flag whether any changed loop falls in the Figure 1b/c survey region +in_survey_region <- diff_df$chr == survey_chrom & + diff_df$start1 >= survey_chromstart & diff_df$start2 <= survey_chromend +diff_df$in_fig1bc_region <- in_survey_region + +data.table::fwrite(diff_df, file.path(data_dir, "processed/border_fix_diff.csv")) + + +# Summary + manuscript sentence ------------------------------------------- + +n_total <- nrow(merged_df) +n_changed <- sum(changed) +n_rescued <- sum(diff_df$rescued) +n_lost <- sum(!diff_df$rescued & diff_df$loopCommunity_new == "") +n_in_fig <- sum(diff_df$in_fig1bc_region) + +cat(glue("\n=== Genome-wide border-adjustment fix impact ===\n")) +cat(glue("Total loops compared: {n_total}\n")) +cat(glue("Loops with changed loopCommunity assignment: {n_changed} ", + "({round(100 * n_changed / n_total, 3)}%)\n")) +cat(glue(" - Rescued (previously unassigned, now assigned): {n_rescued}\n")) +cat(glue(" - Lost (previously assigned, now unassigned): {n_lost}\n")) +cat(glue(" - Changed loops overlapping Fig 1b/c region ", + "(chr1:110-115Mb): {n_in_fig}\n\n")) + +cat(glue( + "During revision we identified and fixed a bug in assignCommunities()'s ", + "border-adjustment step ... [see commits fdc1731/527eb78/ce26123]. ", + "Comparing pre-fix and post-fix outputs on the full genome-wide K562 ", + "loop set (identical inputs, seed pinned), {n_changed} of {n_total} ", + "loops ({round(100 * n_changed / n_total, 3)}%) changed community ", + "assignment, {n_rescued} of which were previously-orphaned border loops ", + "now correctly assigned. {ifelse(n_in_fig == 0, ", + "'None of these loops fall within the region shown in Figure 1b/c.', ", + "glue('{n_in_fig} of these loops fall within the region shown in ", + "Figure 1b/c and that panel should be regenerated.'))}" +), "\n\n") + +sessionInfo() diff --git a/inst/scripts/call_tads_arrowhead.sbatch b/inst/scripts/call_tads_arrowhead.sbatch new file mode 100644 index 0000000..930e2e5 --- /dev/null +++ b/inst/scripts/call_tads_arrowhead.sbatch @@ -0,0 +1,31 @@ +#!/bin/bash +#SBATCH -J arrowhead_tads_0h +#SBATCH -p general +#SBATCH -N 1 +#SBATCH -n 4 +#SBATCH -t 240 +#SBATCH --mem=32g +#SBATCH -o .batch_logs/arrowhead_tads_0h.%j.out +#SBATCH -e .batch_logs/arrowhead_tads_0h.%j.err +#----------------------------------------------------------------------------- +# Calls Arrowhead TADs (SCALE-normalized, 10kb) on the K562 MEGA 0h .hic file, +# matching the resolution used by inst/scripts/tad_null_model_sparsity.R. +# Output is renamed to SCALE_0.bedpe, the path that script expects. +# +# Run from the project root: +# mkdir -p .batch_logs +# sbatch inst/scripts/call_tads_arrowhead.sbatch +#----------------------------------------------------------------------------- +set -euo pipefail + +JUICER_TOOLS=/proj/phanstiel_lab/Software/juicer/scripts/juicer_tools_2.13.07.jar +HIC_FILE=/work/users/j/p/jpflores/projects/loopcity/data/processed/hic/GSE214123_MEGA_K562_WT_0_inter.hic +OUT_DIR=/work/users/j/p/jpflores/projects/loopcity/data/processed/hic/tads + +mkdir -p "$OUT_DIR" + +java -jar "$JUICER_TOOLS" arrowhead -r 10000 -k SCALE "$HIC_FILE" "$OUT_DIR" + +mv "$OUT_DIR/10000_blocks.bedpe" "$OUT_DIR/SCALE_0.bedpe" + +echo "Done. TAD calls written to $OUT_DIR/SCALE_0.bedpe" diff --git a/inst/scripts/tad_null_model_sparsity.R b/inst/scripts/tad_null_model_sparsity.R new file mode 100644 index 0000000..063cc87 --- /dev/null +++ b/inst/scripts/tad_null_model_sparsity.R @@ -0,0 +1,272 @@ +# TAD null-model sparsity comparison --------------------------------------- +# Author: JP Flores +# Date: 2026-07-13 +# Project: loopcity +# Description: Addresses reviewer comment on comparing loopcity's community +# network against a naive null model: assume all merged loop +# anchors within a single TAD form a fully connected graph. +# Quantifies how much sparser the network becomes after the +# loopcity community-detection step, demonstrating specificity +# of called loops versus mere spatial (TAD) proximity. +# +# TADs come from Arrowhead calls on the same K562 MEGA dataset +# used throughout the paper (Bond et al. 2023), 0h condition, +# matching the .hic file used in survey_communities.R and +# benchmark_timing.R. See: +# github.com/mbond0718/MEGA2023/blob/main/scripts/processing/makeTadRanges.R +# for how the source bedpe files were generated (Arrowhead, +# SCALE-normalized, one file per timepoint). +# +# Nested/hierarchical TAD calls (Arrowhead reports sub-TADs +# inside larger domains) are flattened to only the outermost, +# non-nested domain at each locus before computing sparsity, +# so the null model isn't built on an arbitrarily small sub-TAD. +# +# Sparsity is computed with clustAnalytics::internal_density(), +# already a loopcity dependency (used for conductance elsewhere). +# Given a graph and a vertex grouping, internal_density() returns +# observed edges / all possible pairs within each group — i.e. +# exactly the density of the naive-null complete graph's actual +# fill, per TAD. The null's own density is 1 by construction, so +# 1 - internal_density() is directly "how much sparser." +# +# Input: data/processed/hic/GSE214123_MEGA_K562_WT_0_inter.hic +# data/processed/hic/5kbLoops_0.txt (K562 loop calls) +# data/processed/hic/tads/SCALE_0.bedpe (Arrowhead TAD calls, +# 0h condition — copy from MEGA2023 project data if not +# already present in this project's data dir) +# Output: figures/tad_null_model_sparsity.pdf +# data/processed/tad_sparsity.csv +# prints manuscript sentence to console +# --------------------------------------------------------------------------- + + +# Parameters ------------------------------------------------------------ + +project_dir <- "/work/users/j/p/jpflores/projects/loopcity" +data_dir <- file.path(project_dir, "data") +output_dir <- file.path(project_dir, "figures") + +hic_file <- file.path(data_dir, "processed/hic/GSE214123_MEGA_K562_WT_0_inter.hic") +loops_file <- file.path(data_dir, "processed/hic/5kbLoops_0.txt") +tad_file <- file.path(data_dir, "processed/hic/tads/SCALE_0.bedpe") + +## Communities cache from the same pipeline used elsewhere — avoids +## re-running mergeAnchors/connectLoopAnchors/scoreInteractions/assignCommunities +communities_rds <- file.path(data_dir, "processed/loopcity_communities.rds") + +hic_resolution <- 10e3 +hic_norm <- "NONE" +pixel_overlap <- 1 +overlap_dist <- 1e6 +leiden_res <- 0.5 + + +# Libraries --------------------------------------------------------------- + +library(clustAnalytics) +library(data.table) +library(GenomeInfoDb) +library(GenomicRanges) +library(ggplot2) +library(glue) +library(igraph) +library(InteractionSet) +library(loopcity) +library(mariner) + + +# Load / run loopcity pipeline --------------------------------------------- + +if (file.exists(communities_rds)) { + message("Loading cached communities from: ", communities_rds) + communities <- readRDS(communities_rds) +} else { + merged <- data.table::fread(loops_file) |> + mariner::assignToBins(binSize = 10e3) |> + mergeAnchors(pixelOverlap = pixel_overlap) + GenomeInfoDb::seqlevelsStyle(merged) <- "ENSEMBL" + connected <- connectLoopAnchors(merged, overlapDist = overlap_dist) + scored <- scoreInteractions(connected, + hicFile = hic_file, + resolution = hic_resolution, + norm = hic_norm) + communities <- assignCommunities(InteractionSet::interactions(scored), + leidenResolution = leiden_res) + saveRDS(communities, communities_rds) +} + + +# Load TAD calls and flatten nested domains -------------------------------- + +## Arrowhead bedpe: V1:V3 = anchor1 (chr,start,end), V4:V6 = anchor2, +## V12 = score. A TAD's genomic span is anchor1$start to anchor2$end. +## skip = 2: juicer_tools 2.x prefixes output with a "#chr1 x1 x2..." column +## header line and a "# juicer_tools version..." comment line, both of which +## confuse fread's header autodetection if left in. +tad_raw <- data.table::fread(tad_file, skip = 2, header = FALSE) +tad_gr <- GenomicRanges::GRanges( + seqnames = tad_raw$V1, + ranges = IRanges::IRanges(start = tad_raw$V2 + 1, end = tad_raw$V6), + score = tad_raw$V12 +) +GenomeInfoDb::seqlevelsStyle(tad_gr) <- "ENSEMBL" + +## Keep only outermost domains: drop any TAD fully nested within a larger one. +## Iterate largest-to-smallest, keeping a TAD only if it isn't `within` an +## already-kept (necessarily larger) domain. +flattenNestedTads <- function(gr) { + gr <- gr[order(-GenomicRanges::width(gr))] + keep <- logical(length(gr)) + kept_gr <- GenomicRanges::GRanges() + for (i in seq_along(gr)) { + nested <- length(kept_gr) > 0 && + length(IRanges::subsetByOverlaps(gr[i], kept_gr, type = "within")) > 0 + if (!nested) { + keep[i] <- TRUE + kept_gr <- c(kept_gr, gr[i]) + } + } + gr[keep] +} + +top_tads <- flattenNestedTads(tad_gr) +top_tads$tadId <- seq_along(top_tads) + +message(glue( + "{length(tad_gr)} raw TAD calls flattened to ", + "{length(top_tads)} outermost (non-nested) domains" +)) + + +# Build genome-wide anchor graph from final loopcity network --------------- + +## Nodes = every merged anchor (regardless of degree); edges = interactions +## retained by assignCommunities (i.e. the network *after* community +## detection / pruning — the graph the reviewer wants compared to the null). +all_anchors <- InteractionSet::regions(communities) +anc <- InteractionSet::anchorIds(communities) + +edge_df <- data.frame(from = anc$first, to = anc$second) +g <- igraph::graph_from_data_frame( + edge_df, + directed = FALSE, + vertices = data.frame(name = seq_along(all_anchors)) +) + +## Assign each anchor to a flattened top-level TAD (NA if it falls outside +## all TAD calls — those anchors are excluded from the sparsity comparison) +overlaps <- GenomicRanges::findOverlaps(all_anchors, top_tads) +tad_membership <- rep(NA_integer_, length(all_anchors)) +tad_membership[S4Vectors::queryHits(overlaps)] <- + top_tads$tadId[S4Vectors::subjectHits(overlaps)] + +## Restrict to anchors assigned to exactly one TAD +keep_vertices <- which(!is.na(tad_membership)) +g_sub <- igraph::induced_subgraph(g, keep_vertices) +com_sub <- tad_membership[keep_vertices] + +## clustAnalytics's community functions require contiguous 1:k labels (they +## index communities positionally); tadId has gaps because many TADs have no +## overlapping anchors. Remap to contiguous labels for the call, keeping +## tad_ids_present to translate back to genomic TAD identity in the output. +tad_ids_present <- sort(unique(com_sub)) +com_contig <- match(com_sub, tad_ids_present) + +## clustAnalytics::internal_density()/edges_inside() build a weighted copy of +## the graph internally but then (a package bug) compute the edge list from +## the original, unweighted graph if it has no "weight" attribute -- silently +## returning garbage values instead of erroring. Set weight explicitly so +## they never hit that path. +g_sub <- igraph::set_edge_attr(g_sub, "weight", value = 1) + + +# Compute sparsity per TAD -------------------------------------------------- + +## internal_density = observed edges / all possible pairs within the group. +## The naive null (fully connected graph) has density 1 by construction, so +## 1 - density is exactly "how much sparser" loopcity's network is. +tad_density <- clustAnalytics::internal_density(g_sub, com_contig) +tad_edges <- clustAnalytics::edges_inside(g_sub, com_contig) +tad_sizes <- as.integer(table(factor(com_contig, levels = seq_along(tad_ids_present)))) + +tad_sparsity <- data.frame( + tadId = tad_ids_present, + n_anchors = tad_sizes, + possible_edges = choose(tad_sizes, 2), + observed_edges = tad_edges, + density = tad_density +) |> + dplyr::filter(n_anchors >= 2) + +tad_sparsity$sparsity <- 1 - tad_sparsity$density + +data.table::fwrite(tad_sparsity, + file.path(data_dir, "processed/tad_sparsity.csv")) + + +# Summary stats + manuscript sentence --------------------------------------- + +mean_density <- mean(tad_sparsity$density, na.rm = TRUE) +mean_sparsity <- 1 - mean_density + +message(glue( + "Mean density across {nrow(tad_sparsity)} TADs: ", + "{round(mean_density, 4)} ", + "(naive fully-connected null = 1)" +)) + +cat(glue( + "Compared to a naive null model in which all merged loop anchors within ", + "a TAD form a fully connected graph (mean density = 1 by construction), ", + "loopcity's community-detected network had a mean density of only ", + "{round(mean_density, 2)} across {nrow(tad_sparsity)} TADs, a ", + "{round(mean_sparsity * 100, 1)}% reduction in edge density relative to ", + "the naive null." +), "\n\n") + + +# Visualization ------------------------------------------------------------- + +loopcity_color <- "#2a78d6" # blue -- loopcity (observed) + +## Density by TAD-size bin vs. the null's density of 1. A raw +## histogram of density is misleading here: small TADs have few possible +## anchor pairs (choose(k,2)), so density can only take a handful of discrete +## values (e.g. exactly 0 or 1 for a 2-anchor TAD) -- producing spiky, +## multimodal-looking bars that are a combinatorial artifact, not biology. +## Binning by TAD size and boxplotting shows sparsity holds across size +## classes, not just among small, discretization-prone TADs. +tad_sparsity$size_bin <- cut( + tad_sparsity$n_anchors, + breaks = c(1, 3, 6, 10, Inf), + labels = c("2-3", "4-6", "7-10", "11+") +) + +## Direct-label each box with its median so the headline numbers are legible +## without having to visually read quartiles off the plot. +box_medians <- tad_sparsity |> + dplyr::group_by(size_bin) |> + dplyr::summarize(med = median(density), .groups = "drop") + +p2 <- ggplot(tad_sparsity, aes(x = size_bin, y = density)) + + geom_boxplot(fill = loopcity_color, alpha = 0.6, outlier.alpha = 0.3, width = 0.6) + + geom_hline(yintercept = 1, linetype = "dashed", color = "grey40") + + annotate("text", x = Inf, y = 1, label = "naive null density = 1", + hjust = 1.05, vjust = 1.5, size = 3, color = "grey40") + + geom_text(data = box_medians, aes(x = size_bin, y = med, label = sprintf("%.2f", med)), + vjust = -1, size = 3.2, fontface = "bold", color = "grey20", + inherit.aes = FALSE) + + scale_y_continuous(limits = c(0, 1.08), breaks = seq(0, 1, 0.25)) + + labs(x = "Merged anchors per TAD", y = "Internal density (observed / possible edges)", + title = "Network density by TAD size") + + theme_bw(base_size = 11) + +pdf(file.path(output_dir, "tad_null_model_sparsity.pdf"), width = 6, height = 4.5) +print(p2) +dev.off() + + +# Session info -------------------------------------------------------------- + +sessionInfo() diff --git a/man/GM12878_10KbLoops.Rd b/man/GM12878_10KbLoops.Rd new file mode 100644 index 0000000..81e78b2 --- /dev/null +++ b/man/GM12878_10KbLoops.Rd @@ -0,0 +1,31 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/data.R +\docType{data} +\name{GM12878_10KbLoops} +\alias{GM12878_10KbLoops} +\title{GM12878 10Kb Hi-C Loop Calls (Chr22)} +\format{ +A \code{GInteractions} object with 175 interactions and 9 +metadata columns: +\describe{ + \item{color}{RGB color string for visualization} + \item{APScoreAvg}{Average anchor pair score} + \item{ProbabilityofEnrichment}{Probability of enrichment at loop site} + \item{RegAPScoreAvg}{Regional average anchor pair score} + \item{Avg_diffMaxNeihgboor_1}{Average difference from max neighbor, + anchor 1} + \item{Avg_diffMaxNeihgboor_2}{Average difference from max neighbor, + anchor 2} + \item{avg}{Average Hi-C contact count} + \item{std}{Standard deviation of Hi-C contact count} + \item{value}{Loop enrichment score} +} +} +\usage{ +data(GM12878_10KbLoops) +} +\description{ +A \code{GInteractions} object of 10kb resolution Hi-C loop calls for +chromosome 22 from the GM12878 human B-lymphoblastoid cell line. +} +\keyword{datasets} diff --git a/man/GM12878_10KbLoopsChr21.Rd b/man/GM12878_10KbLoopsChr21.Rd new file mode 100644 index 0000000..c307295 --- /dev/null +++ b/man/GM12878_10KbLoopsChr21.Rd @@ -0,0 +1,31 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/data.R +\docType{data} +\name{GM12878_10KbLoopsChr21} +\alias{GM12878_10KbLoopsChr21} +\title{GM12878 10Kb Hi-C Loop Calls (Chr21)} +\format{ +A \code{GInteractions} object with 141 interactions and 9 +metadata columns: +\describe{ + \item{color}{RGB color string for visualization} + \item{APScoreAvg}{Average anchor pair score} + \item{ProbabilityofEnrichment}{Probability of enrichment at loop site} + \item{RegAPScoreAvg}{Regional average anchor pair score} + \item{Avg_diffMaxNeihgboor_1}{Average difference from max neighbor, + anchor 1} + \item{Avg_diffMaxNeihgboor_2}{Average difference from max neighbor, + anchor 2} + \item{avg}{Average Hi-C contact count} + \item{std}{Standard deviation of Hi-C contact count} + \item{value}{Loop enrichment score} +} +} +\usage{ +data(GM12878_10KbLoopsChr21) +} +\description{ +A \code{GInteractions} object of 10kb resolution Hi-C loop calls for +chromosome 21 from the GM12878 human B-lymphoblastoid cell line. +} +\keyword{datasets} diff --git a/man/assignCommunities.Rd b/man/assignCommunities.Rd index ad12ad0..e752907 100644 --- a/man/assignCommunities.Rd +++ b/man/assignCommunities.Rd @@ -36,11 +36,8 @@ GInteractions object with community metadata Assign communities to interactions } \examples{ -hicFile <- "inst/extdata/GM12878_chr22.hic" - -mergedLoops <- mergeAnchors(GM12878_10KbLoops, 1) -connections <- connectLoopAnchors(mergedLoops, 1e6) -scores <- scoreInteractions(connections, hicFile, mergedLoops) - -assignCommunities(interactions(scores)) +data(GM12878_10KbLoopsChr21) +loops <- GM12878_10KbLoopsChr21 +loops$score <- loops$value +assignCommunities(loops, pruneUnder = 0) } diff --git a/man/connectLoopAnchors.Rd b/man/connectLoopAnchors.Rd index 271ea73..3a1949a 100644 --- a/man/connectLoopAnchors.Rd +++ b/man/connectLoopAnchors.Rd @@ -9,8 +9,8 @@ connectLoopAnchors(loops, overlapDist) \arguments{ \item{loops}{GInteractions object of loops} -\item{overlapDist}{an integer representing the farthest distance in base pairs -to make a connection} +\item{overlapDist}{an integer representing the farthest distance in +base pairs to make a connection} } \value{ GInteractions object @@ -19,7 +19,9 @@ GInteractions object connectLoopAnchors } \examples{ -mergedLoops <- mergeAnchors(GM12878_10KbLoops, 1) - -connectLoopAnchors(mergedLoops, 1e6) +if (requireNamespace("loopcityData", quietly = TRUE)) { + data(GM12878_10KbLoops, package = "loopcityData") + mergedLoops <- mergeAnchors(GM12878_10KbLoops, 1) + connectLoopAnchors(mergedLoops, 1e6) +} } diff --git a/man/figures/Fig1.png b/man/figures/Fig1.png index 82f96b5..7334225 100644 Binary files a/man/figures/Fig1.png and b/man/figures/Fig1.png differ diff --git a/man/loopcity-package.Rd b/man/loopcity-package.Rd index 49022f7..020dc71 100644 --- a/man/loopcity-package.Rd +++ b/man/loopcity-package.Rd @@ -13,12 +13,26 @@ Tools for assigning and visualizing chromatin community interactions from Hi-C d \seealso{ Useful links: \itemize{ - \item \url{http://www.sarmapar.com/loopcity/} + \item \url{https://sarmapar.com/loopcity/} + \item \url{https://github.com/sarmapar/loopcity} + \item Report bugs at \url{https://github.com/sarmapar/loopcity/issues} } } \author{ \strong{Maintainer}: Sarah Parker \email{smp3800@gmail.com} (\href{https://orcid.org/0000-0002-2700-3979}{ORCID}) +Authors: +\itemize{ + \item Sarah Parker \email{smp3800@gmail.com} (\href{https://orcid.org/0000-0002-2700-3979}{ORCID}) +} + +Other contributors: +\itemize{ + \item Howard Hughes Medical Institute [funder] + \item National Institutes of Health [funder] + \item National Science Foundation [funder] +} + } \keyword{internal} diff --git a/man/mergeAnchors.Rd b/man/mergeAnchors.Rd index 5503954..7f43a29 100644 --- a/man/mergeAnchors.Rd +++ b/man/mergeAnchors.Rd @@ -5,7 +5,7 @@ \title{Merge neighboring loop anchor positions into one representative anchor.} \usage{ -mergeAnchors(loops, pixelOverlap = 1, dropDups = T) +mergeAnchors(loops, pixelOverlap = 1, dropDups = TRUE) } \arguments{ \item{loops}{GInteractions object of loops from Hi-C data} @@ -19,6 +19,10 @@ be merged} \item{dropDups}{logical determining whether duplicates should be dropped from output} } +\value{ +A GInteractions object with loop anchors merged into representative +positions. +} \description{ The representative anchor is chosen as the mode of a set of neighboring anchors. If there are no modes, the middle-most anchor of the @@ -28,6 +32,8 @@ The representative anchor is chosen as the mode of a set of break ties. } \examples{ -mergeAnchors(loops = GM12878_10KbLoops, pixelOverlap = 1) - +if (requireNamespace("loopcityData", quietly = TRUE)) { + data(GM12878_10KbLoops, package = "loopcityData") + mergeAnchors(loops = GM12878_10KbLoops, pixelOverlap = 1) +} } diff --git a/man/plotHicCommunities.Rd b/man/plotHicCommunities.Rd index 69c196d..988e0d2 100644 --- a/man/plotHicCommunities.Rd +++ b/man/plotHicCommunities.Rd @@ -44,16 +44,14 @@ If a vector of starts and ends is provided, one page will be generated per region. } \examples{ +data(GM12878_10KbLoops) hicFile <- system.file("extdata", "GM12878_chr22.hic", package = "loopcity") - -\donttest{ communities <- GM12878_10KbLoops |> mergeAnchors() |> connectLoopAnchors(overlapDist = 1e6) |> scoreInteractions(hicFile = hicFile) |> InteractionSet::interactions() |> assignCommunities() - plotHicCommunities( pdfName = tempfile(fileext = ".pdf"), communities = communities, @@ -63,4 +61,3 @@ plotHicCommunities( ends = 30e6 ) } -} diff --git a/man/scoreInteractions.Rd b/man/scoreInteractions.Rd index 599e787..f4ee89d 100644 --- a/man/scoreInteractions.Rd +++ b/man/scoreInteractions.Rd @@ -1,10 +1,25 @@ % Generated by roxygen2: do not edit by hand % Please edit documentation in R/scoreInteractions.R -\name{scoreInteractions,GInteractions,character-method} +\name{scoreInteractions} +\alias{scoreInteractions} \alias{scoreInteractions,GInteractions,character-method} \alias{scoreInteractions,InteractionArray,missing-method} \title{Score loop counts using enrichment of foreground over background} \usage{ +scoreInteractions( + x, + hicFile, + fgSize = 0, + bg = "bowtie", + bgGap = 0, + bgSize = 5, + pseudo = 5, + pruneUnder = 0, + truncate = FALSE, + resolution, + norm = "NONE" +) + \S4method{scoreInteractions}{GInteractions,character}( x, hicFile, @@ -14,7 +29,7 @@ bgSize = 5, pseudo = 5, pruneUnder = 0, - truncate = F, + truncate = FALSE, resolution, norm = "NONE" ) @@ -62,17 +77,20 @@ See also `qc_scoreInteractions`} \item{interactions}{GInteraction object of interactions from Hi-C data} } +\value{ +An object of the same class as \code{x} with a \code{score} column +added to \code{mcols}. +} \description{ Given a hi-c file or InteractionArray, pseudocounts are added to raw counts and the median of the foreground over the background is added as a score for each given interaction. } \examples{ -hicFile <- "inst/extdata/GM12878_chr22.hic" - +data(GM12878_10KbLoops) +hicFile <- system.file("extdata", "GM12878_chr22.hic", package = "loopcity") mergedLoops <- mergeAnchors(GM12878_10KbLoops, 1) connections <- connectLoopAnchors(mergedLoops, 1e6) - -scoreInteractions(connections, hicFile, mergedLoops) +scoreInteractions(connections, hicFile) } diff --git a/tests/testthat/test-assignCommunities.R b/tests/testthat/test-assignCommunities.R index 2347e01..6c1ba9b 100644 --- a/tests/testthat/test-assignCommunities.R +++ b/tests/testthat/test-assignCommunities.R @@ -1,16 +1,22 @@ ### Setup inputs +skip_if_not_installed("loopcityData") +data("GM12878_10KbLoops", package = "loopcityData", envir = environment()) + hicfile <- system.file("extdata", "GM12878_chr22.hic", - package = "loopcity") + package = "loopcity" +) loops <- GM12878_10KbLoops GenomeInfoDb::seqlevelsStyle(loops) <- "ENSEMBL" connections <- connectLoopAnchors(loops, 1e6) -scores <- scoreInteractions(x = connections, - hicFile = hicfile, - norm = "NONE", - pseudo = 5) +scores <- scoreInteractions( + x = connections, + hicFile = hicfile, + norm = "NONE", + pseudo = 5 +) loops_with_scores <- InteractionSet::interactions(scores) @@ -19,41 +25,47 @@ test_that("assignCommunities generates GInteractions", { expect_s4_class("GInteractions") }) -test_that("Various cluster types can be used",{ +test_that("Various cluster types can be used", { ## fast greedy assignCommunities(loops_with_scores, - clusterType = "fast_greedy") |> + clusterType = "fast_greedy" + ) |> expect_s4_class("GInteractions") ## walktrap assignCommunities(loops_with_scores, - clusterType = "walktrap") |> + clusterType = "walktrap" + ) |> expect_s4_class("GInteractions") ## infomap assignCommunities(loops_with_scores, - clusterType = "infomap") |> + clusterType = "infomap" + ) |> expect_s4_class("GInteractions") ## label prop assignCommunities(loops_with_scores, - clusterType = "label_prop") |> + clusterType = "label_prop" + ) |> expect_s4_class("GInteractions") ## edge betweenness assignCommunities(loops_with_scores, - clusterType = "edge_betweenness") |> + clusterType = "edge_betweenness" + ) |> expect_warning("highest modularity") ## invalid cluster type assignCommunities(loops_with_scores, - clusterType = "random_mess") |> + clusterType = "random_mess" + ) |> expect_error("not a valid clustering algorithm") }) test_that("loops are checked correctly", { ## loops must be a GInteractions object - test_df <- data.frame(x = c(1,2,3), y = c(3,6,9)) + test_df <- data.frame(x = c(1, 2, 3), y = c(3, 6, 9)) test_char <- "test" assignCommunities(test_char) |> @@ -66,7 +78,7 @@ test_that("loops are checked correctly", { test_that("scores are checked correctly", { ## `scores` must be a column in `loops` or a numeric vector ## of the same length - test_scores <- c(2,4,6,8,10) + test_scores <- c(2, 4, 6, 8, 10) assignCommunities(loops_with_scores, score = test_scores) |> expect_error("must be the same length") @@ -75,7 +87,7 @@ test_that("scores are checked correctly", { expect_error("must be a numeric vector") ## Providing scores overwrites `score` column in `loops` - ##TODO: this test isn't actually testing that the right scores are being used... + ## TODO: this test isn't actually testing that the right scores are being used... test_scores_correct <- 1:length(loops_with_scores) assignCommunities(loops_with_scores, score = test_scores_correct) |> expect_s4_class("GInteractions") @@ -90,27 +102,78 @@ test_that("scores are checked correctly", { set.seed(555) scoresInLoopsOutput <- assignCommunities(loops_with_scores) scoresSeperateOutput <- assignCommunities(loops_without_scores, - scores = as.numeric( - loops_with_scores$score)) + scores = as.numeric( + loops_with_scores$score + ) + ) + + expect_identical( + scoresInLoopsOutput$loopCommunity, + scoresSeperateOutput$loopCommunity + ) +}) + +test_that("loopCommunity is the intersection of its anchors' communities", { + result <- assignCommunities(loops_with_scores) + + anc1_comm <- InteractionSet::regions(result)[ + InteractionSet::anchorIds(result)$first + ]$anchorCommunity + + anc2_comm <- InteractionSet::regions(result)[ + InteractionSet::anchorIds(result)$second + ]$anchorCommunity + + expected <- mapply(intersect, anc1_comm, anc2_comm) + expect_identical(result$loopCommunity, expected) +}) - expect_identical(scoresInLoopsOutput$loopCommunity, - scoresSeperateOutput$loopCommunity) +test_that("border-adjustment rescues an anchor when it scores better with the neighboring community", { + ## Deterministic example, independent of real Hi-C data or clustering + ## randomness: 4 anchors on chr1, A-B assigned to community 1, C-D to + ## community 2. The B-C loop (score 10) is the only loop connecting + ## border anchor B to community 2, and scores higher than the only loop + ## inside community 1 (A-B, score 5) -- so B should be rescued into both + ## communities. This directly exercises the logic flagged by reviewer 2: + ## `comm` must be passed explicitly, and community membership must be + ## read from `anchorCommunity` (not the not-yet-assigned `loopCommunity`). + regions <- GenomicRanges::GRanges( + "chr1", + IRanges::IRanges(start = c(1, 20001, 40001, 60001), width = 10000) + ) + regions$anchorCommunity <- c(1, 1, 2, 2) + + gi <- InteractionSet::GInteractions( + anchor1 = c(1, 2, 3), + anchor2 = c(2, 3, 4), + regions = regions + ) + gi$score <- c(5, 10, 5) + + borderNode <- InteractionSet::regions(gi)[2] + + result <- .compareNodeToComm(gi, borderNode, left = FALSE, comm = 1) + + expect_identical(sort(unlist(result$anchorCommunity)), c(1, 2)) }) test_that("leiden_resolution checks work properly", { ## Give a warning if leidenResolution is provided but not used assignCommunities(loops_with_scores, - clusterType = "fast_greedy", - leidenResolution = 0.5) |> + clusterType = "fast_greedy", + leidenResolution = 0.5 + ) |> expect_warning("will not be used") ## Give a warning if more than one resolution is provided assignCommunities(loops_with_scores, - leidenResolution = c(1,2,3)) |> + leidenResolution = c(1, 2, 3) + ) |> expect_warning("first will be used") ## Error if leidenResolution is not numeric assignCommunities(loops_with_scores, - leidenResolution = "test") |> + leidenResolution = "test" + ) |> expect_error("must be numeric") }) diff --git a/tests/testthat/test-connectLoopAnchors.R b/tests/testthat/test-connectLoopAnchors.R index f0b9da9..e586a88 100644 --- a/tests/testthat/test-connectLoopAnchors.R +++ b/tests/testthat/test-connectLoopAnchors.R @@ -1,19 +1,28 @@ -regions5k <- GenomicRanges::GRanges(seqnames = "chr1", - IRanges::IRanges(seq(5e3,95e3,5e3), - width = 5001)) +regions5k <- GenomicRanges::GRanges( + seqnames = "chr1", + IRanges::IRanges(seq(5e3, 95e3, 5e3), + width = 5001 + ) +) test_that("overlaps the correct distance", { - inters <- InteractionSet::GInteractions(anchor1 = c(1,3,5), - anchor2 = c(8,10,8), - regions = regions5k) + inters <- InteractionSet::GInteractions( + anchor1 = c(1, 3, 5), + anchor2 = c(8, 10, 8), + regions = regions5k + ) - expected <- InteractionSet::GInteractions(anchor1 = c(1,3,5,8), - anchor2 = c(3,5,8,10), - regions = regions5k, - source = c("added", - "added", - "original", - "added")) |> + expected <- InteractionSet::GInteractions( + anchor1 = c(1, 3, 5, 8), + anchor2 = c(3, 5, 8, 10), + regions = regions5k, + source = c( + "added", + "added", + "original", + "added" + ) + ) |> InteractionSet::reduceRegions() connectLoopAnchors(inters, 11e3) |> diff --git a/tests/testthat/test-mergeAnchors.R b/tests/testthat/test-mergeAnchors.R index 673060b..17b7ab4 100644 --- a/tests/testthat/test-mergeAnchors.R +++ b/tests/testthat/test-mergeAnchors.R @@ -1,63 +1,80 @@ # create GRanges object to pull regions from -regions5k <- GenomicRanges::GRanges(seqnames = "chr1", - IRanges::IRanges(seq(5e3,95e3,5e3), width = 5001)) +regions5k <- GenomicRanges::GRanges( + seqnames = "chr1", + IRanges::IRanges(seq(5e3, 95e3, 5e3), width = 5001) +) -regions10k <- GenomicRanges::GRanges(seqnames = "chr1", - IRanges::IRanges(seq(10e3,90e3,10e3), width = 10001)) +regions10k <- GenomicRanges::GRanges( + seqnames = "chr1", + IRanges::IRanges(seq(10e3, 90e3, 10e3), width = 10001) +) test_that("anchors merge to correct positions", { - #create input ranges and interactions - inters <- InteractionSet::GInteractions(anchor1 = c(2,1,3), - anchor2 = c(5,7,7), - regions = regions5k) - - #create expected output ranges and interactions with duplicates - expectedDups <- InteractionSet::GInteractions(anchor1 = rep(2,3), - anchor2 = c(5,7,7), - regions = regions5k) |> + # create input ranges and interactions + inters <- InteractionSet::GInteractions( + anchor1 = c(2, 1, 3), + anchor2 = c(5, 7, 7), + regions = regions5k + ) + + # create expected output ranges and interactions with duplicates + expectedDups <- InteractionSet::GInteractions( + anchor1 = rep(2, 3), + anchor2 = c(5, 7, 7), + regions = regions5k + ) |> InteractionSet::reduceRegions() - #create expected output ranges and interactions without duplicates - expectedNoDups <- InteractionSet::GInteractions(anchor1 = c(2,2), - anchor2 = c(5,7), - regions = regions5k) |> + # create expected output ranges and interactions without duplicates + expectedNoDups <- InteractionSet::GInteractions( + anchor1 = c(2, 2), + anchor2 = c(5, 7), + regions = regions5k + ) |> InteractionSet::reduceRegions() - #generate outputs using mergeAnchors + # generate outputs using mergeAnchors outputDups <- mergeAnchors(inters, pixelOverlap = 1, dropDups = F) mergeAnchors(inters, pixelOverlap = 1) |> expect_message("Duplicates dropped.") outputNoDups <- mergeAnchors(inters, pixelOverlap = 1) - #compare expected output and function output + # compare expected output and function output expect_identical(outputDups, expectedDups) expect_identical(outputNoDups, expectedNoDups) }) test_that("there is no message if there are no duplicates", { - #create input interactions - inters <- InteractionSet::GInteractions(anchor1 = c(2,1,3), - anchor2 = c(5,7,9), - regions = regions5k) + # create input interactions + inters <- InteractionSet::GInteractions( + anchor1 = c(2, 1, 3), + anchor2 = c(5, 7, 9), + regions = regions5k + ) - #expect no message, dropDups is TRUE but there are no duplicates + # expect no message, dropDups is TRUE but there are no duplicates expect_no_message(mergeAnchors(inters, pixelOverlap = 1)) - }) test_that("mixed range widths give an error when using pixelOverlap", { - #create input interactions with mixed widths - inters <- c(InteractionSet::GInteractions(anchor1 = c(1,3), - anchor2 = c(7,9), - regions = regions5k), - InteractionSet::GInteractions(anchor1 = c(1,3), - anchor2 = c(7,9), - regions = regions10k)) - - #expect an error about different widths - expect_error(mergeAnchors(inters, pixelOverlap = 1), - "same width to use pixelOverlap") + # create input interactions with mixed widths + inters <- c( + InteractionSet::GInteractions( + anchor1 = c(1, 3), + anchor2 = c(7, 9), + regions = regions5k + ), + InteractionSet::GInteractions( + anchor1 = c(1, 3), + anchor2 = c(7, 9), + regions = regions10k + ) + ) + # expect an error about different widths + expect_error( + mergeAnchors(inters, pixelOverlap = 1), + "same width to use pixelOverlap" + ) }) - diff --git a/tests/testthat/test-scoreInteractions.R b/tests/testthat/test-scoreInteractions.R index 8521ed8..dec7482 100644 --- a/tests/testthat/test-scoreInteractions.R +++ b/tests/testthat/test-scoreInteractions.R @@ -1,39 +1,54 @@ -test_that("scores are generated", { - hicfile <- system.file("extdata", "GM12878_chr22.hic", - package = "loopcity") +skip_if_not_installed("loopcityData") +data("GM12878_10KbLoops", package = "loopcityData", envir = environment()) - mergedLoops <- - GM12878_10KbLoops |> - mariner::assignToBins(binSize = 10e3) |> - mergeAnchors(1) - GenomeInfoDb::seqlevelsStyle(mergedLoops) <- "ENSEMBL" +test_that("scores are generated", { + hicfile <- system.file("extdata", "GM12878_chr22.hic", + package = "loopcity" + ) - connections <- connectLoopAnchors(mergedLoops, 1e6) + mergedLoops <- + GM12878_10KbLoops |> + mariner::assignToBins(binSize = 10e3) |> + mergeAnchors(1) + GenomeInfoDb::seqlevelsStyle(mergedLoops) <- "ENSEMBL" - scores <- scoreInteractions(x = connections, hicFile = hicfile, norm = "NONE", - loopCalls = mergedLoops) + connections <- connectLoopAnchors(mergedLoops, 1e6) + scores <- scoreInteractions( + x = connections, hicFile = hicfile, norm = "NONE" + ) - assignCommunities(InteractionSet::interactions(scores)) + assignCommunities(InteractionSet::interactions(scores)) }) test_that("loops with multiple widths are handled correctly", { ## generate loops of various widths - regions5k <- GenomicRanges::GRanges(seqnames = "22", - IRanges::IRanges(seq(5e3,95e3,5e3), - width = 5001)) - - regions10k <- GenomicRanges::GRanges(seqnames = "22", - IRanges::IRanges(seq(10e3,90e3,10e3), - width = 10001)) - - loops <- c(InteractionSet::GInteractions(anchor1 = c(1,3), - anchor2 = c(7,9), - regions = regions5k), - InteractionSet::GInteractions(anchor1 = c(1,3), - anchor2 = c(7,9), - regions = regions10k)) + regions5k <- GenomicRanges::GRanges( + seqnames = "22", + IRanges::IRanges(seq(5e3, 95e3, 5e3), + width = 5001 + ) + ) + + regions10k <- GenomicRanges::GRanges( + seqnames = "22", + IRanges::IRanges(seq(10e3, 90e3, 10e3), + width = 10001 + ) + ) + + loops <- c( + InteractionSet::GInteractions( + anchor1 = c(1, 3), + anchor2 = c(7, 9), + regions = regions5k + ), + InteractionSet::GInteractions( + anchor1 = c(1, 3), + anchor2 = c(7, 9), + regions = regions10k + ) + ) }) - diff --git a/vignettes/img/Fig1.png b/vignettes/img/Fig1.png index 82f96b5..7334225 100644 Binary files a/vignettes/img/Fig1.png and b/vignettes/img/Fig1.png differ diff --git a/vignettes/loopcity.Rmd b/vignettes/loopcity.Rmd index 92be741..9fa4d00 100644 --- a/vignettes/loopcity.Rmd +++ b/vignettes/loopcity.Rmd @@ -4,15 +4,15 @@ author: "Sarah Parker" date: '`r Sys.Date()`' output: rmarkdown::html_vignette vignette: > - %\VignetteIndexEntry{loopcity} + %\VignetteIndexEntry{Using loopcity} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r knitr-opts, include = FALSE} knitr::opts_chunk$set( - collapse = TRUE, - comment = "#>" + collapse = TRUE, + comment = "#>" ) ``` @@ -21,19 +21,31 @@ library(loopcity) library(loopcityData) ``` +# Introduction + +`loopcity` is an R package for identifying and visualizing chromatin loop +communities from Hi-C data. Starting from a set of loop calls and a `.hic` +contact matrix, it merges redundant anchors, connects neighboring anchors into +a network, scores each interaction by Hi-C signal enrichment, and assigns loops +to communities using graph clustering. Results can be visualized as colored +arches overlaid on Hi-C maps. + # Installation -This package can be installed via Github: +`loopcity` can be installed from Bioconductor: ```{r installPackage, eval=FALSE} -if (!requireNamespace("remotes", quietly = TRUE)) - install.packages("remotes") -remotes::install_github("sarmapar/loopcity") +if (!requireNamespace("BiocManager", quietly = TRUE)) { + install.packages("BiocManager") +} +BiocManager::install("loopcity") ``` # loopcity Pipeline -While `loopcity` functions can be used in isolation, they are more powerful used in sequence. There are 4 main steps to call communities from a `.hic` file and a `GInteractions` object of loops. +While `loopcity` functions can be used in isolation, they are more powerful +used in sequence. There are 4 main steps to call communities from a `.hic` +file and a `GInteractions` object of loops. ## Starting with a loop file @@ -51,7 +63,8 @@ loops <- read.table(loopFile) |> mariner::assignToBins(binSize = 10e3) ``` -For this vignette we use the GM12878 10 kb loop calls shipped with `loopcityData`: +For this vignette we use the GM12878 10 kb loop calls shipped with +`loopcityData`: ```{r preprocessLoops} data(GM12878_10KbLoops) @@ -60,15 +73,24 @@ loops <- GM12878_10KbLoops ## Step 1: Merge Loop Anchors using `mergeAnchors` -Since loop calling can be imprecise, using loop anchors as-is can lead to multiple anchors representing one biologically relevant region. By default, any duplicate loops created via this merging are dropped, but these can be kept by setting `dropDups = F`. +Since loop calling can be imprecise, using loop anchors as-is can lead to +multiple anchors representing one biologically relevant region. By default, +any duplicate loops created via this merging are dropped, but these can be +kept by setting `dropDups = FALSE`. ```{r mergingLoops} mergedLoops <- mergeAnchors(loops) ``` -In this example, since one bin is 10Kb and `pixelOverlap` is 1 by default, all loop anchors within 1 pixel or 10Kb from each other will be grouped and merged into one representative anchor. The most common or the middle-most of these anchors is chosen as the new representative position for all neighboring anchors. Increasing `pixelOverlap` to `n` means all loop anchors within `n` pixels, or `n` \* `binSize` base pairs will be merged into one position. +In this example, since one bin is 10Kb and `pixelOverlap` is 1 by default, +all loop anchors within 1 pixel or 10Kb from each other will be grouped and +merged into one representative anchor. The most common or the middle-most of +these anchors is chosen as the new representative position for all neighboring +anchors. Increasing `pixelOverlap` to `n` means all loop anchors within `n` +pixels, or `n` \* `binSize` base pairs will be merged into one position. -The output, `mergedLoops`, is a `GInteractions` object with all the same information from `loops` but with adjusted anchor positions. +The output, `mergedLoops`, is a `GInteractions` object with all the same +information from `loops` but with adjusted anchor positions. ```{r mergedOutput} mergedLoops @@ -76,13 +98,22 @@ mergedLoops ## Step 2: Add connections using `connectLoopAnchors` -Since loop callers are limited to a binary classification of "loop" or "not a loop", we can instead create all possible loops and then create a threshold of what should be kept as a "true" loop based on a set of flexible parameters. The first step in this process is to create a set of all possible loops. The `connectLoopAnchors` function creates a new set of loops, connecting all anchors to every anchor within `overlapDist`. In this example, all anchors within 1Mb of each other are connected to create new loops. +Since loop callers are limited to a binary classification of "loop" or +"not a loop", we can instead create all possible loops and then create a +threshold of what should be kept as a "true" loop based on a set of flexible +parameters. The first step in this process is to create a set of all possible +loops. The `connectLoopAnchors` function creates a new set of loops, +connecting all anchors to every anchor within `overlapDist`. In this example, +all anchors within 1Mb of each other are connected to create new loops. ```{r connectingAnchors} connections <- connectLoopAnchors(mergedLoops, overlapDist = 1e6) ``` -This results in a new `GInteractions` object, with original metadata removed and a new metadata column named "source" where all interactions that were present in the original set of loops are labeled as "original" and added loops are labeled as "added". +This results in a new `GInteractions` object, with original metadata removed +and a new metadata column named "source" where all interactions that were +present in the original set of loops are labeled as "original" and added +loops are labeled as "added". ```{r connectionsOutput} connections @@ -90,40 +121,70 @@ connections ## Step 3: Score original and new loops using `scoreInteractions` -These added and original connections are then scored based on the enrichment of the counts of the loop pixel compared to the local background. The shape and size of this local background can be customized using the parameters `bg` for the shape, and `bgSize`, `bgGap`, and `fgSize` for the size. +These added and original connections are then scored based on the enrichment +of the counts of the loop pixel compared to the local background. The shape +and size of this local background can be customized using the parameters `bg` +for the shape, and `bgSize`, `bgGap`, and `fgSize` for the size. -This normalizes contact counts between loop sizes, since shorter loops typically have higher raw counts. To account for noise in long range loops with low counts, a pseudocount value of 5 is added to all raw counts before calculating an enrichment score. This pseudocount value can be modified by changing the `pseudo` value. +This normalizes contact counts between loop sizes, since shorter loops +typically have higher raw counts. To account for noise in long range loops +with low counts, a pseudocount value of 5 is added to all raw counts before +calculating an enrichment score. This pseudocount value can be modified by +changing the `pseudo` value. ```{r scoringLoops} hicFile <- system.file("extdata", "GM12878_chr22.hic", package = "loopcity") scores <- scoreInteractions(connections, - hicFile = hicFile) + hicFile = hicFile +) ``` -This results in a new `InteractionArray` containing the scores for all interactions and metadata such as the name of the .hic file the counts came from and the number of pseudocounts added. To get just the `GInteraction` object with the new `score` column, we can run `interactions(scores)` +This results in a new `InteractionArray` containing the scores for all +interactions and metadata such as the name of the .hic file the counts came +from and the number of pseudocounts added. To get just the `GInteraction` +object with the new `score` column, we can run `interactions(scores)` ```{r scores} -# the `scores` object is an 'InteractionArray' containing extra metadata information -scores +# the `scores` object is an 'InteractionArray' containing +# extra metadata information +scores ``` ```{r score_interactions} -# the interactions of the `InteractionArray` gives back a `GInteractions` object +# the interactions of the `InteractionArray` gives back +# a `GInteractions` object InteractionSet::interactions(scores) ``` ## Step 4: Assign anchors and loop communities using `assignCommunities` -The `GInteractions` object stored inside `score` can then be used to assign communities. This function builds a network where the nodes are each loop anchor and the edge weights are the enrichment scores calculated in the previous step. By default, all loops with a score less than the median score of original loops are removed. This threshold is printed in a message and stored in the metadata of the output object in a variable called `pruningValue`. This value can be manually set by passing a value to the `pruneUnder` parameter. - -This network is then passed through a clustering algorithm. The default algorithm is leiden, but other options include "fast_greedy", "walktrap", "infomap", "label_prop", or "edge_betweenness". This assigns each loop anchor to a community, and any nodes that are closely bordering a neighboring community are evaluated to determine if they should be included in both communities. If the average score of the border anchor and all anchors in the neighboring community is greater than the average score of all interactions between all the anchors in the current community, then the anchor is assigned to both the original and neighboring community. Loops are assigned to a community if both anchors of the loop are in the same community. +The `GInteractions` object stored inside `score` can then be used to assign +communities. This function builds a network where the nodes are each loop +anchor and the edge weights are the enrichment scores calculated in the +previous step. By default, all loops with a score less than the median score +of original loops are removed. This threshold is printed in a message and +stored in the metadata of the output object in a variable called +`pruningValue`. This value can be manually set by passing a value to the +`pruneUnder` parameter. + +This network is then passed through a clustering algorithm. The default +algorithm is leiden, but other options include "fast_greedy", "walktrap", +"infomap", "label_prop", or "edge_betweenness". This assigns each loop +anchor to a community, and any nodes that are closely bordering a neighboring +community are evaluated to determine if they should be included in both +communities. If the average score of the border anchor and all anchors in +the neighboring community is greater than the average score of all +interactions between all the anchors in the current community, then the +anchor is assigned to both the original and neighboring community. Loops are +assigned to a community if both anchors of the loop are in the same community. ```{r assignCommunities} communities <- assignCommunities(InteractionSet::interactions(scores)) ``` -To see the pruning value used to remove edges from the network, run the following code +To see the pruning value used to remove edges from the network, run the +following code ```{r pruningValue} S4Vectors::metadata(communities)$pruningValue @@ -131,7 +192,10 @@ S4Vectors::metadata(communities)$pruningValue ## Full Workflow -These four functions can be combined into one step by piping the results from one function into the next. This results in the same final `GInteractions` object which contains the merged anchor locations, enrichment scores for each loop, anchor communities, and loop communities. +These four functions can be combined into one step by piping the results from +one function into the next. This results in the same final `GInteractions` +object which contains the merged anchor locations, enrichment scores for each +loop, anchor communities, and loop communities. ```{r fullWorkflow} communities <- loops |> @@ -144,20 +208,28 @@ communities <- loops |> ## Visualizing Communities -We also provide a function to visualize these communities using the final `GInteractions` object returned from `assignCommunities` and a .hic file. This returns a pdf of the given region(s), with each region on each page. The figure contains a heatmap of hic counts with original loops annotated by black squares and added loops annotated by grey squares and arches with their height set to the calculated enrichment score colored by their assigned community. +We also provide a function to visualize these communities using the final +`GInteractions` object returned from `assignCommunities` and a .hic file. +This returns a pdf of the given region(s), with each region on each page. +The figure contains a heatmap of hic counts with original loops annotated +by black squares and added loops annotated by grey squares and arches with +their height set to the calculated enrichment score colored by their +assigned community. ```{r plot, eval = F} -plotHicCommunities(pdfName = "plot.pdf", - communities = communities, - hicFile = hicFile, - norm = "NONE", - chroms = 22, - starts = 25e6, - ends = 30e6, - zmax = 10) +plotHicCommunities( + pdfName = "plot.pdf", + communities = communities, + hicFile = hicFile, + norm = "NONE", + chroms = 22, + starts = 25e6, + ends = 30e6, + zmax = 10 +) ``` -![](img/plot.png){width="100%"} +![Hi-C map with loop communities visualized as colored arches](img/plot.png){width="100%"} ```{r sessionInfo} sessionInfo() diff --git a/vignettes/plotgardener-integration.Rmd b/vignettes/plotgardener-integration.Rmd index f41bd38..b4f717d 100644 --- a/vignettes/plotgardener-integration.Rmd +++ b/vignettes/plotgardener-integration.Rmd @@ -1,5 +1,5 @@ --- -title: "Visualizing loop communities with `plotgardener`" +title: "Integrating loop communities with genomic data tracks" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Integrating loop communities with genomic data tracks} @@ -9,8 +9,8 @@ vignette: > ```{r setup, include = FALSE} knitr::opts_chunk$set( - collapse = TRUE, - comment = "#>" + collapse = TRUE, + comment = "#>" ) ``` @@ -38,19 +38,22 @@ library(GenomeInfoDb) library(InteractionSet) library(mariner) -hic_file <- K562_hic() +hic_file <- K562_hic() loops_file <- K562_loops() loops <- data.table::fread(loops_file) |> mariner::assignToBins(binSize = 10e3) GenomeInfoDb::seqlevelsStyle(loops) <- "ENSEMBL" -merged <- mergeAnchors(loops, pixelOverlap = 1) -connected <- connectLoopAnchors(merged, overlapDist = 1e6) -scored <- scoreInteractions(connected, hicFile = hic_file, - resolution = 10e3, norm = "NONE") +merged <- mergeAnchors(loops, pixelOverlap = 1) +connected <- connectLoopAnchors(merged, overlapDist = 1e6) +scored <- scoreInteractions(connected, + hicFile = hic_file, + resolution = 10e3, norm = "NONE" +) communities <- assignCommunities(interactions(scored), - leidenResolution = 0.5) + leidenResolution = 0.5 +) ``` ## Visualize with plotgardener @@ -64,14 +67,14 @@ library(plotgardener) library(glue) ## File paths -ctcf_bw <- K562_CTCF_bw() +ctcf_bw <- K562_CTCF_bw() h3k27ac_bw <- K562_H3K27ac_bw() ## Region of interest — matches Figure 1b/c -hic_chrom <- "1" -bw_chrom <- "chr1" +hic_chrom <- "1" +bw_chrom <- "chr1" survey_chromstart <- 110e6 -survey_chromend <- 115e6 +survey_chromend <- 115e6 ## Shared parameters params <- pgParams( @@ -87,30 +90,43 @@ params <- pgParams( ## Filter communities to region survey_region <- GenomicRanges::GRanges( seqnames = hic_chrom, - ranges = IRanges::IRanges(start = survey_chromstart, - end = survey_chromend) + ranges = IRanges::IRanges( + start = survey_chromstart, + end = survey_chromend + ) ) a1_in <- GenomicRanges::countOverlaps( - anchors(communities, "first"), survey_region) > 0 + anchors(communities, "first"), survey_region +) > 0 a2_in <- GenomicRanges::countOverlaps( - anchors(communities, "second"), survey_region) > 0 + anchors(communities, "second"), survey_region +) > 0 loops_region <- communities[a1_in & a2_in] loops_region <- loops_region[ which(loops_region$score > 0 & - lengths(loops_region$loopCommunity) > 0) + lengths(loops_region$loopCommunity) > 0) ] ## Assign community colours loopcity_colors <- function(n) { - rep(c("chartreuse3", "deepskyblue3", "darkorange", - "darkorchid2", "deeppink3"), length.out = n) + rep(c( + "chartreuse3", "deepskyblue3", "darkorange", + "darkorchid2", "deeppink3" + ), length.out = n) } -comm_labels <- vapply(loops_region$loopCommunity, - function(x) if (length(x) == 0L) NA_integer_ - else as.integer(x[[1L]]), - integer(1L)) -comm_idx <- as.integer(factor(comm_labels)) -loops_region$color <- loopcity_colors(max(comm_idx, na.rm=TRUE))[comm_idx] +comm_labels <- vapply( + loops_region$loopCommunity, + function(x) { + if (length(x) == 0L) { + NA_integer_ + } else { + as.integer(x[[1L]]) + } + }, + integer(1L) +) +comm_idx <- as.integer(factor(comm_labels)) +loops_region$color <- loopcity_colors(max(comm_idx, na.rm = TRUE))[comm_idx] loops_region$color[is.na(comm_labels)] <- "grey70" loops_region$score_lg <- as.numeric(log2(loops_region$score)) @@ -118,7 +134,7 @@ loops_region$score_lg <- as.numeric(log2(loops_region$score)) pdf("figures/survey_communities.pdf", width = 9, height = 7.5) pageCreate(width = 9, height = 7.5, showGuides = FALSE) -x <- 0.5 +x <- 0.5 width <- 8.0 ## Hi-C contact map @@ -127,14 +143,17 @@ hic_plot <- plotHicRectangle( chrom = hic_chrom, zrange = c(0, 225), x = x, y = 0.4, width = width, height = 2.2 ) -annoHeatmapLegend(hic_plot, x = x - 0.25, y = 0.4, - width = 0.07, height = 0.75, - just = c("left","top"), default.units = "inches") +annoHeatmapLegend(hic_plot, + x = x - 0.25, y = 0.4, + width = 0.07, height = 0.75, + just = c("left", "top"), default.units = "inches" +) ## Community arches score_range <- range(loops_region$score_lg, na.rm = TRUE) plotPairsArches( - loops_region, params = params, chrom = hic_chrom, + loops_region, + params = params, chrom = hic_chrom, flip = TRUE, clip = TRUE, archHeight = "score_lg", range = score_range, x = x, y = 2.66, width = width, height = 0.75, @@ -142,30 +161,44 @@ plotPairsArches( ) ## CTCF signal -plotSignal(ctcf_bw, params = params, chrom = bw_chrom, - x = x, y = 3.42, width = width, height = 0.45, - fill = "#253494", linecolor = "#253494", - range = c(0, 200), scale = FALSE) -plotText("CTCF", x = x, y = 3.37, fontsize = 7, fontface = "bold", - fontcolor = "#253494", just = c("left","bottom")) +plotSignal(ctcf_bw, + params = params, chrom = bw_chrom, + x = x, y = 3.42, width = width, height = 0.45, + fill = "#253494", linecolor = "#253494", + range = c(0, 200), scale = FALSE +) +plotText("CTCF", + x = x, y = 3.37, fontsize = 7, fontface = "bold", + fontcolor = "#253494", just = c("left", "bottom") +) ## H3K27ac signal -plotSignal(h3k27ac_bw, params = params, chrom = bw_chrom, - x = x, y = 4.05, width = width, height = 0.45, - fill = "#807DBA", linecolor = "#807DBA", - range = calcSignalRange(h3k27ac_bw, chrom = bw_chrom, - chromstart = survey_chromstart, - chromend = survey_chromend, - assembly = "hg38", negData = FALSE), - scale = FALSE) -plotText("H3K27ac", x = x, y = 4.00, fontsize = 7, fontface = "bold", - fontcolor = "#807DBA", just = c("left","bottom")) +plotSignal(h3k27ac_bw, + params = params, chrom = bw_chrom, + x = x, y = 4.05, width = width, height = 0.45, + fill = "#807DBA", linecolor = "#807DBA", + range = calcSignalRange(h3k27ac_bw, + chrom = bw_chrom, + chromstart = survey_chromstart, + chromend = survey_chromend, + assembly = "hg38", negData = FALSE + ), + scale = FALSE +) +plotText("H3K27ac", + x = x, y = 4.00, fontsize = 7, fontface = "bold", + fontcolor = "#807DBA", just = c("left", "bottom") +) ## Genes and genome label -plotGenes(params = params, x = x, y = 4.54, - width = width, height = 0.5, fontsize = 6) -annoGenomeLabel(plot = hic_plot, x = x, y = 5.07, - scale = "bp", fontsize = 6) +plotGenes( + params = params, x = x, y = 4.54, + width = width, height = 0.5, fontsize = 6 +) +annoGenomeLabel( + plot = hic_plot, x = x, y = 5.07, + scale = "bp", fontsize = 6 +) dev.off() ``` @@ -177,4 +210,4 @@ object, coordinates are guaranteed to be aligned across the full figure. ```{r sessionInfo} sessionInfo() -``` \ No newline at end of file +```