improved findBreakpointOverlaps() performance when both query and subject are large GRanges

Author: Daniel Cameron

DESCRIPTION

@@ -1,8 +1,8 @@
 Package: StructuralVariantAnnotation
 Type: Package
 Title: VariantionAnnoration for Structural Variants
-Version: 0.7.1
-Date: 2018-08-23
+Version: 0.7.2
+Date: 2018-01-25
 Author: Daniel Cameron
 Maintainer: Daniel Cameron <daniel.l.cameron@gmail.com>
 Description: StructuralVariantAnnotation contains useful helper
@@ -20,7 +20,8 @@ Depends:
 	stats,
     R (>= 3.2.2),
     S4Vectors (>= 0.10.0),
-    VariantAnnotation (>= 1.16.3)
+    VariantAnnotation (>= 1.16.3),
+    dplyr
 Imports:
     assertthat,
     Biostrings (>= 2.40.0),

R/BreakpointGRanges.R

@@ -31,32 +31,15 @@ partner <- function(gr) {
 #'
 #'@export
 findBreakpointOverlaps <- function(query, subject, maxgap=-1L, minoverlap=0L, ignore.strand=FALSE, sizemargin=0.25, restrictMarginToSizeMultiple=0.5) {
-	hitdf <- as.data.frame(findOverlaps(query, subject, maxgap=maxgap, minoverlap=minoverlap, type="any", select="all", ignore.strand=ignore.strand), row.names=NULL)
-	# instead of running findOverlaps(partner(query), partner(subject), ...
-	# we can reduce our runtime cost by just performing partner index lookups
-	# partner lookups
-	subjectPartnerIndexLookup <- seq_along(names(subject))
-	names(subjectPartnerIndexLookup) <- names(subject)
-	queryPartnerIndexLookup <- seq_along(names(query))
-	names(queryPartnerIndexLookup) <- names(query)
-	phitdf <- data.frame(
-		queryHits=queryPartnerIndexLookup[query$partner[hitdf$queryHits]],
-		subjectHits=subjectPartnerIndexLookup[subject$partner[hitdf$subjectHits]])
-	hits <- rbind(hitdf, phitdf, make.row.names=FALSE)
+	hits <- dplyr::bind_rows(
+		as.data.frame(findOverlaps(query, subject, maxgap=maxgap, minoverlap=minoverlap, type="any", select="all", ignore.strand=ignore.strand), row.names=NULL),
+		as.data.frame(findOverlaps(query, subject, maxgap=maxgap, minoverlap=minoverlap, type="any", select="all", ignore.strand=ignore.strand), row.names=NULL))
 	# we now want to do:
 	# hits <- hits[duplicated(hits),] # both breakends match
 	# but for large hit sets (such as focal false positive loci) we run out of memory (>32GB)
-	# instead, we sort then check that we match the previous record
-	hits <- hits[base::order(hits$queryHits, hits$subjectHits), ]
-	lg <- function(x) {
-		if (length(x) == 0) {
-			return(x)
-		} else {
-			return(c(-1, x[1:(length(x)-1)])) # -1 to ensure FALSE match instead of NA match
-		}
-	}
-	isDup <- hits$queryHits == lg(hits$queryHits) & hits$subjectHits == lg(hits$subjectHits)
-	hits <- hits[isDup,]
+	# instead, we sort then check that we match the next record
+	hits = hits %>% dplyr::arrange(queryHits, subjectHits) %>%
+		dplyr::filter(!is.na(dplyr::lead(.$queryHits)) & !is.na(dplyr::lead(.$subjectHits)) & dplyr::lead(.$queryHits) == .$queryHits & dplyr::lead(.$subjectHits) == .$subjectHits)
 	if (!is.null(sizemargin) && !is.na(sizemargin)) {
 		# take into account confidence intervals when calculating event size
 		callwidth <- .distance(query, partner(query))
@@ -343,29 +326,28 @@ calculateBlastHomology <- function(gr, ref, db, anchorLength=150) {
 #' @param gr breakpoint GRanges object. Can contain both breakpoint and single breakend SV records
 #'
 #'@export
-breakpointGRangesToVCF <- function(gr) {
-	if (is.null(gr$insSeq)) {
-		gr$insSeq = rep("", length(gr))
-	}
-	nominalgr = GRanges(seqnames=seqnames(gr), ranges=IRanges(start=(end(gr) + start(gr)) / 2, width=1))
-	if (is.null(gr$REF)) {
-		gr$REF = rep("N", length(gr))
-	}
-	gr$ALT[is.na(gr$ALT)] = ""
-	if (is.null(gr$ALT)) {
-		gr$ALT = rep("", length(gr))
-	}
-	gr$ALT[is.na(gr$ALT)] = ""
-	gr$ALT[gr$ALT == ""] = .toVcfBreakendNotationAlt(gr)[gr$ALT == ""]
-	ciposstart = start(gr) - start(nominalgr)
-	ciposend = end(gr) - end(nominalgr)
-	vcf = VCF(rowRanges=nominalgr, collapsed=FALSE)
-	fixeddf = data.frame(
-		ALT=gr$ALT,
-		REF=gr$REF,
-		QUAL=gr$QUAL,
-		FILTER=gr$FILTER)
-
-	VCF(rowRanges = GRanges(), colData = DataFrame(), exptData = list(header = VCFHeader()), fixed = DataFrame(), info = DataFrame(), geno = SimpleList(), ..., collapsed=FALSE, verbose = FALSE
-}
+# breakpointGRangesToVCF <- function(gr) {
+# 	if (is.null(gr$insSeq)) {
+# 		gr$insSeq = rep("", length(gr))
+# 	}
+# 	nominalgr = GRanges(seqnames=seqnames(gr), ranges=IRanges(start=(end(gr) + start(gr)) / 2, width=1))
+# 	if (is.null(gr$REF)) {
+# 		gr$REF = rep("N", length(gr))
+# 	}
+# 	gr$ALT[is.na(gr$ALT)] = ""
+# 	if (is.null(gr$ALT)) {
+# 		gr$ALT = rep("", length(gr))
+# 	}
+# 	gr$ALT[is.na(gr$ALT)] = ""
+# 	gr$ALT[gr$ALT == ""] = .toVcfBreakendNotationAlt(gr)[gr$ALT == ""]
+# 	ciposstart = start(gr) - start(nominalgr)
+# 	ciposend = end(gr) - end(nominalgr)
+# 	vcf = VCF(rowRanges=nominalgr, collapsed=FALSE)
+# 	fixeddf = data.frame(
+# 		ALT=gr$ALT,
+# 		REF=gr$REF,
+# 		QUAL=gr$QUAL,
+# 		FILTER=gr$FILTER)
+# 	VCF(rowRanges = GRanges(), colData = DataFrame(), exptData = list(header = VCFHeader()), fixed = DataFrame(), info = DataFrame(), geno = SimpleList(), ..., collapsed=FALSE, verbose = FALSE)
+# }
 

tests/testthat/test-BreakpointGRanges.R

@@ -211,6 +211,18 @@ test_that("calculateBlastHomology", {
 	#bh <- calculateBlastHomology(gr, hg19, "~/blastdb/16SMicrobial")
 
 })
+test_that("performance_test_partner", {
+	n = 10000
+	gr = GRanges(
+		seqnames="1",
+		ranges=IRanges(start=1:(2*n), width=1),
+		partner=c(paste0(1:n, "o"), paste0(1:n, "h")))
+	names(gr)=c(paste0(1:n, "h"), paste0(1:n, "o"))
+	tictoc::tic(paste0("Start", n))
+	pgr = partner(gr)
+	tictoc::toc()
+})
+