[Pruning] Seismic

R/converters.R

@@ -909,7 +909,7 @@ read.juncs = function(rafile,
  ra = readRDS(rafile)
  ## validity check written for "junctions" class
  return(Junction$new(ra))
- } else if (grepl('(.bedpe$)', rafile)){
+ } else if (grepl('bedpe(\\.gz)?$', rafile)){
  ra.path = rafile
  cols = c('chr1', 'start1', 'end1', 'chr2', 'start2', 'end2', 'name', 'score', 'str1', 'str2')
 

R/eventCallers.R

@@ -1051,9 +1051,13 @@ annotate_walks = function(walks)
 #' Shortcut to call all simple and complex event types on JaBbA graph using standard settings on all event callers. 
 #' 
 #' @param gg gGraph
+#' @param verbose
+#' @param mark
+#' @param QRP annotate quasi reciprocal events
+#' @param seismic (logical or character) if TRUE then annotates seismic amplifications using a graph-based approach, if 'Rosswog' then the Rosswog et al. algorithm is used and events are annotated as seismic.rosswog, and if 'both' then will annotate both "seismic" (graph-based) and "seismic.rosswog" (Rosswog et al. algorithm)
 #' @return gGraph with nodes and edges annotated with complex events in their node and edge metadata and in the graph meta data field $events 
 #' @export
-events = function(gg, verbose = TRUE, mark = FALSE, QRP = FALSE)
+events = function(gg, verbose = TRUE, mark = FALSE, QRP = FALSE, seismic = FALSE)
 {
  gg = gg %>% simple(mark = TRUE)
  if (verbose)
@@ -1088,7 +1092,20 @@ events = function(gg, verbose = TRUE, mark = FALSE, QRP = FALSE)
  if (verbose)
  message('Finished qrp')
  }
-
+
+ if (seismic == TRUE | seismic == 'both' | seismic == 'Rosswog'){
+ Rosswog = FALSE
+ if (seismic == 'Rosswog'){
+ Rosswog = TRUE
+ }
+ gg = gg %>% seismic(mark = TRUE, Rosswog = Rosswog)
+ if (verbose)
+ message('Finished seismic')
+ if (seismic == 'both'){
+ gg = gg %>% seismic(mark = TRUE, Rosswog = TRUE)
+ }
+ }
+
  ev = rbind(
  gg$meta$simple,
  gg$meta$chromothripsis,
@@ -1108,6 +1125,16 @@ events = function(gg, verbose = TRUE, mark = FALSE, QRP = FALSE)
  gg$meta$qrpmin,
  gg$meta$qrpmix, fill = TRUE)[, ev.id := seq_len(.N)]
  }
+ if (seismic == TRUE | seismic == 'both'){
+ ev = rbind(
+ ev,
+ gg$meta$seismic, fill = TRUE)[, ev.id := seq_len(.N)]
+ }
+ if (seismic == 'Rosswog'){
+ ev = rbind(
+ ev,
+ gg$meta$seismic.rosswog, fill = TRUE)[, ev.id := seq_len(.N)]
+ }
 
  gg$set(events = ev)
  return(gg)
@@ -2916,6 +2943,121 @@ qrp = function(gg, thresh = 1e6, max.small = 1e5,
  return(gg)
 
 }
+
+
+#' @name seismic
+#' @description
+#' Identification of seismic amplicons according to the definitions in Rosswog et al.
+#' should be used. 
+#' 
+#' @param gg gGraph
+#' @param amp.thresh (numeric) multiply of rounded ploidy to use for detection of amplification (2). The threshold for amplification is set to >= rounded.ploidy * amp.thresh + 1 where rounded.ploidy is either 2 or 4 (determined using ploidy.thresh).
+#' @param ploidy.thresh (numeric) threshold to use above which ploidy is considered 4 (and up to it, 2). Rosswog et al. used 2 and hence this is the default value.
+#' @param min.internal (numeric) minimal number of internal junctions in a seismic amplification (14).
+#' @param mc.cores (numeric) number of cores to use (1)
+#' @param mark (logical) nodes and edges with color if they are included in a seismic amplification (TRUE)
+#' @param mark.col (character) color to use (if mark set to TRUE) to mark edges and nodes (purple).
+#' @param Rosswog (logical) run the original algorithm by Rosswog et al.
+#' @param ... additional parameters for the Rosswog et al. algorithm (see help menu for gGnome::seismic_rosswog()). This is only relevant is Rosswog = TRUE.
+#' 
+#' @return gg
+#' @export
+seismic = function(gg, amp.thresh = 2, ploidy.thresh = 2, min.internal = 14,
+ mc.cores = 1, mark = TRUE, mark.col = 'purple', Rosswog = FALSE,
+ ...)
+{
+ gg$nodes$mark(seismic = as.integer(NA))
+ gg$edges$mark(seismic = as.integer(NA))
+ gg$set(seismic = data.table())
+ if (!('cn' %in% names(mcols(gg$nodes$gr)))){
+ stop('In order to call seismic amplifications gGraph nodes must contain "cn" values')
+ }
+ if (is.null(gg$meta$ploidy)){
+ ploidy = gg$nodes$dt[!is.na(cn), sum(cn*as.numeric(width))/sum(as.numeric(width))]
+ } else {
+ ploidy = gg$meta$ploidy
+ }
+
+ if (isTRUE(Rosswog)){
+ message('Running the original algorithm by Rosswog et al.')
+ rosswog_calls = seismic_rosswog(gg, ploidy, ...)
+ gg = annotate_rosswog_calls(gg, rosswog_calls)
+ } else {
+
+ # following Rosswog et al. we use 5 as the threshold for samples with ploidy up to 2
+ # and 9 as the threshold for samples with ploidy above 3
+ amp.thresh = ifelse(ploidy <= ploidy.thresh, 2 * amp.thresh + 1, 4 * amp.thresh + 1)
+
+ amplified = (gg$nodes$dt$cn) >= amp.thresh
+
+ gg$clusters(amplified, mode = 'strong')
+
+ cids = unique(gg$nodes$dt$cluster)
+ cids = cids[!is.na(cids)]
+
+ edt = mclapply(cids, function(cid){
+ sg = gg[cluster == cid]
+ if (length(sg$edges[type == 'ALT']) < 2){
+ return(NULL)
+ }
+
+ sg.no.alt = sg$clone()
+ # mark regions with rid (region ID)
+ eids = sg.no.alt$edges$dt$type == 'REF'
+ sg.no.alt$clusters(j = eids)
+ sg$nodes$mark(rid = sg.no.alt$nodes$dt$cluster)
+
+ # for each ALT edge check if both left and right node are in the same region (i.e. it is internal)
+ internal = sg$edges[type == 'ALT']$left$dt$rid == sg$edges[type == 'ALT']$right$dt$rid
+
+ # if it is not internal then check if it falls on the edge of regions
+ ndt = sg$nodes$dt
+ left.ids = ndt[, min(node.id), by = rid]$V1
+ right.ids = ndt[, max(node.id), by = rid]$V1
+ edge.nodes = c(left.ids, right.ids)
+
+ flanking = (sg$edges[type == 'ALT']$left$dt$node.id %in% edge.nodes) & (sg$edges[type == 'ALT']$right$dt$node.id %in% edge.nodes)
+
+ other = !internal & !flanking
+
+ return(data.table(cid = cid, internal.junc = sum(internal), flanking.junc = sum(flanking), other.junc = sum(other)))
+ }, mc.cores = mc.cores)
+
+ edt = rbindlist(edt)
+
+ if (edt[,.N] > 0){
+ edt = edt[internal.junc >= min.internal]
+ cids = edt$cid
+ ev.ids = seq_along(cids)
+ if (length(cids) > 0){
+ for (ev.id in ev.ids){
+ # annotate all edges that have both breakpoints whithin the cluster
+ cid = cids[ev.id]
+ e.idx = which(gg$edges[type == 'ALT']$left$dt$cluster == cid & gg$edges[type == 'ALT']$right$dt$cluster == cid)
+ eids = gg$edges[type == 'ALT'][e.idx]$dt$edge.id
+ gg$edges[eids]$mark(seismic = ev.id)
+
+ # annotate nodes
+ gg$nodes[cluster == cid]$mark(seismic = ev.id)
+ }
+
+ names(ev.ids) = as.character(cids)
+ edt[, footprint := sapply(ev.ids[as.character(cid)], function(x){nodes2footprint.str(gg$nodes[seismic == x])})]
+ edt$type = 'seismic'
+ # change the cluster ID column to "strong" to emphesize that this is a strongly connected cluster
+ setnames(edt, 'cid', 'strong')
+ gg$set(seismic = edt)
+ }
+ }
+ }
+
+ if (mark){
+ gg$edges[!is.na(seismic)]$mark(col = mark.col)
+ gg$nodes[!is.na(seismic)]$mark(col = mark.col)
+ }
+
+ return(gg)
+}
 
 #' @name events.to.gr
 #' @title Extract event annotation as a GRanges
@@ -2935,8 +3077,135 @@ events.to.gr = function(gg){
  warning('No events annotated in this gGraph')
  return(GRanges())
  }
- ggrl = parse.grl(gg$meta$events$footprint)
- mcols(ggrl) = gg$meta$events
+ return(events.dt.to.gr(gg$meta$events))
+}
+
+#' @name events.dt.to.gr
+#' @title convert events annotated data.table to a GRanges
+#'
+#' @author Alon Shaiber
+#' @param ev (data.table) events annotation data.table (gGraph$meta$events)
+#' @return GRanges containing ranges of annotated events along with all metadata from gg$meta$events
+#' @export
+events.dt.to.gr = function(ev){
+ ggrl = parse.grl(ev$footprint)
+ mcols(ggrl) = ev
  ggr = grl.unlist(ggrl)
  return(ggr)
 }
+
+
+#' @name seismic_rosswog
+#' @description
+#' Identification of seismic amplicons using the algorithm provided by Rosswog et al. 2021
+#' 
+#' @param gg gGraph
+#' @param ploidy ploidy value
+#' @param rosswog_dir path to a local clone of the github repository provided by Rosswog et al. (https://github.com/seismicon/SeismicAmplification)
+#' @param chrBands GRanges object (see https://github.com/seismicon/SeismicAmplification for details). If nothing is provided then the hg19 version (with no chr prefix) is used.
+#' @param minInternalSVs (numeric) see https://github.com/seismicon/SeismicAmplification (14 by default)
+#' @param cnvTol (numeric) see https://github.com/seismicon/SeismicAmplification (5e3 by default)
+#' 
+#' @return list(amplicons = GRanges, svs = data.table)
+#' @export
+seismic_rosswog = function(gg, ploidy = 2, rosswog_dir = NULL, chrBands = NULL,
+ minInternalSVs = 14, cnvTol = 5e3){
+ if (!is.null(rosswog_dir)){
+ if (dir.exists(rosswog_dir)){
+ rosswog_scripts = paste0(rosswog_dir, '/seismic_amplification_detection/seismic_amplification_detection.R')
+ if (file.exists(rosswog_scripts)){
+ source(rosswog_scripts)
+ breaks = gg$nodes$gr
+ if (!('cn' %in% names(mcols(gg$nodes$gr)))){
+ stop('In order to call seismic amplifications gGraph nodes must contain "cn" values')
+ }
+ alts = gg$junctions[type == 'ALT']
+ if (length(alts) == 0){
+ message('No junctions found in the graph so no seismic amplification')
+ return(list(amplicons = GRanges(), svs = data.table()))
+ }
+ jgrl = alts$grl
+ chrs = as.vector(unlist(seqnames(jgrl)))
+ starts = unlist(start(jgrl))
+ odd = seq(1,2*length(alts),2)
+ even = seq(2,2*length(alts),2)
+ jdt = data.table(chr1 = chrs[odd],
+ chr2 = chrs[even],
+ bp1 = starts[odd],
+ bp2 = starts[even],
+ edge.id = alts$dt$edge.id)
+ if (is.null(chrBands)){
+ message('No chrBands provided so using the default hg19 (with no chr prefix)')
+ chrBands = fread(paste0(rosswog_dir, '/seismic_amplification_detection/chromosome_bands_hg19.txt'))
+ chrBands = chrBands[, `:=`(seqnames = gsub('chr', '', chrom),
+ start = chromStart,
+ end = chromEnd)]
+ chrBands = chrBands[seqnames %in% seqlevels(breaks)]
+ bands = trim(gUtils::dt2gr(chrBands,
+ seqlengths = seqlengths(breaks)[intersect(chrBands$seqnames, seqlevels(breaks))]))
+ }
+ rosswog_calls = detect_seismic_amplification(cnv=breaks, sv=jdt, chrBands=bands,
+ ploidy = ploidy,
+ minInternalSVs = minInternalSVs,
+ cnvTol = cnvTol)
+ return(rosswog_calls)
+ }
+ }
+ }
+}
+
+
+#' @name annotate_rosswog_calls
+#' @description
+#' Takes the output of the Rosswog et al. seismic amplification detection algorithm
+#' and annotates the gGraph accordingly
+#' 
+#' @param gg gGraph
+#' @param rosswog_calls the output of seismic_rosswog()
+#' @param mark.as.rosswog (logical) if set to TRUE then seismic events are annotated as "seismic.rosswog" (TRUE), otherwise annotated as "seismic". This is meant to allow the user to annotate both graph-based "seismic" events and Rosswog et al. "seismic.rosswog" events
+#' @return gg
+annotate_rosswog_calls = function(gg, rosswog_calls, mark.as.rosswog = TRUE){
+ # mark edges and nodes
+ if (length(rosswog_calls$amplicons) > 0){
+ if (isTRUE(mark.as.rosswog)){
+ gg$edges[rosswog_calls$svs$edge.id]$mark(seismic.rosswog = rosswog_calls$svs$amplicon_id)
+ sgr = (gg$nodes$gr %*% rosswog_calls$amplicons[,'id'])
+ gg$nodes[sgr$node.id]$mark(seismic.rosswog = sgr$id)
+ } else {
+ gg$edges[rosswog_calls$svs$edge.id]$mark(seismic = rosswog_calls$svs$amplicon_id)
+ sgr = (gg$nodes$gr %*% rosswog_calls$amplicons[,'id'])
+ gg$nodes[sgr$node.id]$mark(seismic = sgr$id)
+ }
+ cols = c('id', 'nSegments', 'medianCN', 'maxCN', 'size_amplicon', 'nChrs_amplicon', 'nRegions_amplicon',
+ 'medianCN_amplicon', 'cnSpan_amplicon', 'cnStates_amplicon', 'nSegments_amplicon',
+ 'nSVs_amplicon', 'nSVsInternal_amplicon')
+ edt = gr2dt(rosswog_calls$amplicons %Q% (!duplicated(id)))[, ..cols]
+ setnames(edt, 'id', 'cid') # avoid using "id" field
+ name = ifelse(isTRUE(mark.as.rosswog), 'seismic.rosswog', 'seismic')
+ edt[, footprint := sapply(cid, function(x){nodes2footprint.str(gg$nodes[get(name) == x])})]
+ edt[, ev.id := .I]
+ edt$type = name
+ # change the cluster ID column to "strong" to emphesize that this is a strongly connected cluster
+ if (isTRUE(mark.as.rosswog)){
+ gg$set(seismic.rosswog = edt)
+ } else {
+ gg$set(seismic = edt)
+ }
+ }
+ return(gg)
+}
+
+#' @name nodes2footprint.str
+#' @description
+#' Takes a gNode and returns character footprint parsable by parse.gr
+#' 
+#' @param n gNode
+#' @param stripstrand (logical) string the strand information (TRUE)
+#' @return character
+nodes2footprint.str = function(n, stripstrand = TRUE){
+ footprint = n$footprint
+ if (stripstrand){
+ footprint = gr.stripstrand(footprint)
+ }
+ return(paste(gr.string(footprint), collapse = ';'))
+}

tests/testthat/test_gGnome_event_callers.R

@@ -5,11 +5,11 @@ library(gTrack)
 gg.jabba = gG(jabba = system.file('extdata/hcc1954', 'jabba.rds', package="gGnome"))
 
 test_that('event caller', {
- ## run events caller with qrp
-############# not running this currently since it is erroring out (See: https://app.travis-ci.com/github/mskilab/gGnome/builds/237955222 )
- gg.events = events(gg.jabba, QRP = TRUE)
+ ## run events caller with qrp and with seismic
+ gg.events = events(gg.jabba, QRP = TRUE, seismic = TRUE)
  expect_is(gg.events$meta$events, 'data.table')
  expect_true(any(grepl('qrp', gg.events$meta$events$type)))
+ expect_true(any(grepl('seismic', gg.events$meta$events$type)))
 
  ## run events caller without qrp
  gg.events.no.qrp = events(gg.jabba, QRP = FALSE)
@@ -109,3 +109,12 @@ test_that('microhomology', {
  expect_error(microhomology(gg, fa))
 
 })
+
+test_that('seismic', {
+ # test the Rosswog et al. algorithm version
+ # clone the repo
+ seismic_dir = paste0(tempdir(), '/', 'SeismicAmplification')
+ system(paste0('git clone https://github.com/ShaiberAlon/SeismicAmplification.git ', seismic_dir))
+ s = seismic(gg.jabba, Rosswog = TRUE, rosswog_dir = seismic_dir)
+ expect_error(seismic(gG()))
+})