R/BreakpointGRanges.R

Defines functions .traversable_segments findTransitiveCalls .findOverlaps_queryIns_subjectDup findInsDupOverlaps simpleEventLength simpleEventType breakpointGRangesToVCF .toVcfBreakendNotationAlt calculateReferenceHomology .constrict extractReferenceSequence extractBreakpointSequence pairs2breakpointgr .assertValidBreakpointGRanges breakpointgr2pairs countBreakpointOverlaps .distance findBreakpointOverlaps hasPartner partner

Documented in breakpointgr2pairs breakpointGRangesToVCF calculateReferenceHomology countBreakpointOverlaps extractBreakpointSequence extractReferenceSequence findBreakpointOverlaps findInsDupOverlaps findTransitiveCalls hasPartner pairs2breakpointgr partner simpleEventLength simpleEventType

#' GRanges representing the breakend coordinates of
#' structural variants
#' #@export
#setClass("BreakpointGRanges", contains="GRanges")
#' Partner breakend for each breakend.
#'
#' @details
#' All breakends must have their partner breakend included
#' in the GRanges.
#'
#' @param gr GRanges object of SV breakends
#' @param selfPartnerSingleBreakends treat single breakends as their own partner.
#' @return A GRanges object in which each entry is the partner breakend of
#' those in the input object.
#' @examples
#' #reading in a VCF file as \code{vcf}
#' vcf.file <- system.file("extdata", "gridss.vcf", package = "StructuralVariantAnnotation")
#' vcf <- VariantAnnotation::readVcf(vcf.file, "hg19")
#' #parsing \code{vcf} to GRanges object \code{gr}
#' gr <- breakpointRanges(vcf)
#' #output partner breakend of each breakend in \code{gr}
#' partner(gr)
#'@export
partner <- function(gr, selfPartnerSingleBreakends=FALSE) {
  .assertValidBreakpointGRanges(gr, allowSingleBreakends=selfPartnerSingleBreakends)
  return(gr[ifelse(selfPartnerSingleBreakends & is.na(gr$partner), names(gr), gr$partner),])
}
#' Determines whether this breakend has a valid partner in this
#' GRanges
#'
#' @param gr GRanges object of SV breakends
#' @param selfPartnerSingleBreakends treat single breakends as their own partner.
#' @return True/False for each row in the breakpoint GRanges
#' @examples
#' #Subset to chromosome 6 intra-chromosomal events \code{vcf}
#' vcf.file <- system.file("extdata", "COLO829T.purple.sv.ann.vcf.gz",
#'   package = "StructuralVariantAnnotation")
#' vcf <- VariantAnnotation::readVcf(vcf.file)
#' gr <- breakpointRanges(vcf)
#' gr <- gr[seqnames(gr) == "6"]
#' # We now need to filter out inter-chromosomal events to ensure
#' # our GRanges doesn't contain any breakpoints whose partner
#' # has already been filtered out and no longer exists in the GRanges.
#' gr <- gr[hasPartner(gr)]
#'@export
hasPartner <- function(gr, selfPartnerSingleBreakends=FALSE) {
  return((is.na(gr$partner) & selfPartnerSingleBreakends) |
    (!is.na(gr$partner) & gr$partner %in% names(gr)))
}

#' Finding overlapping breakpoints between two breakpoint sets
#'
#' @details
#' \code{findBreakpointOverlaps()} is an efficient adaptation of \code{findOverlaps-methods()}
#' for breakend ranges. It searches for overlaps between breakpoint objects, and return a
#' matrix including index of overlapping ranges as well as error stats.
#' All breakends must have their partner breakend included in the \code{partner}
#' field. A valid overlap requires that breakends on boths sides meets the overlapping
#' requirements.
#'
#' See GenomicRanges::findOverlaps-methods for details of overlap calculation.
#'
#' @param query,subject Both of the input objects should be GRanges objects.
#' Unlike \code{findOverlaps()}, \code{subject} cannot be ommitted. Each breakpoint
#' must be accompanied with a partner breakend, which is also in the GRanges, with the
#' partner's id recorded in the \code{partner} field.
#' See GenomicRanges::findOverlaps-methods for details.
#' @param maxgap,minoverlap Valid overlapping thresholds of a maximum gap and a minimum
#' overlapping positions between breakend intervals. Both should be scalar integers. maxgap
#' allows non-negative values, and minoverlap allows positive values.
#' See GenomicRanges::findOverlaps-methods for details.
#' @param ignore.strand Default value is FALSE. strand information is ignored when set to
#' TRUE.
#' See GenomicRanges::findOverlaps-methods for details.
#' @param sizemargin Error margin in allowable size to prevent matching of events
#' of different sizes, e.g. a 200bp event matching a 1bp event when maxgap is
#' set to 200.
#' @param restrictMarginToSizeMultiple Size restriction multiplier on event size.
#' The default value of 0.5 requires that the breakpoint positions can be off by
#' at maximum, half the event size. This ensures that small deletion do actually
#' overlap at least one base pair.
#' @examples
#' #reading in VCF files
#' query.file <- system.file("extdata", "gridss-na12878.vcf", package = "StructuralVariantAnnotation")
#' subject.file <- system.file("extdata", "gridss.vcf", package = "StructuralVariantAnnotation")
#' query.vcf <- VariantAnnotation::readVcf(query.file, "hg19")
#' subject.vcf <- VariantAnnotation::readVcf(subject.file, "hg19")
#' #parsing vcfs to GRanges objects
#' query.gr <- breakpointRanges(query.vcf)
#' subject.gr <- breakpointRanges(subject.vcf)
#' #find overlapping breakpoint intervals
#' findBreakpointOverlaps(query.gr, subject.gr)
#' findBreakpointOverlaps(query.gr, subject.gr, ignore.strand=TRUE)
#' findBreakpointOverlaps(query.gr, subject.gr, maxgap=100, sizemargin=0.5)
#' @return A dataframe containing index and error stats of overlapping breakpoints.
#'@export
findBreakpointOverlaps <- function(query, subject, maxgap=-1L, minoverlap=0L, ignore.strand=FALSE, sizemargin=NULL, restrictMarginToSizeMultiple=NULL) {
  .assertValidBreakpointGRanges(query)
  .assertValidBreakpointGRanges(subject)
  pquery = partner(query)
  squery = partner(subject)
  localhits = findOverlaps(query, subject, maxgap=maxgap, minoverlap=minoverlap, type="any", select="all", ignore.strand=ignore.strand)
  remotehits = findOverlaps(pquery, squery, maxgap=maxgap, minoverlap=minoverlap, type="any", select="all", ignore.strand=ignore.strand)
  ## duplicated() version:
  #hits = Hits(c(S4Vectors::queryHits(localhits), S4Vectors::queryHits(remotehits)), c(S4Vectors::subjectHits(localhits), S4Vectors::subjectHits(remotehits)), nLnode=nLnode(localhits), nRnode=nRnode(localhits), sort.by.query=TRUE)
  #hits = hits[duplicated(hits)]
  
  ## intersect() version:
  hits = BiocGenerics::intersect(localhits, remotehits)
  
  ## dplyr() version:
  #hits <- dplyr::bind_rows(
  #	as.data.frame(localhits, row.names=NULL),
  #	as.data.frame(remotehits, row.names=NULL))
  #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)
  
  ## dplyr() exploiting the sorted nature of the findOverlaps():
  #hits = Hits(c(S4Vectors::queryHits(localhits), S4Vectors::queryHits(remotehits)), c(S4Vectors::subjectHits(localhits), S4Vectors::subjectHits(remotehits)), nLnode=nLnode(localhits), nRnode=nRnode(localhits), sort.by.query=TRUE)
  #queryLead  = dplyr::lead(S4Vectors::queryHits(hits))
  #querySubject  = dplyr::lead(S4Vectors::queryHits(hits))
  #hits = hits[
  #	!is.na(queryLead) &d
  #	!is.na(querySubject) &
  #	queryLead == S4Vectors::queryHits(hits) &
  #	querySubject == S4Vectors::subjectHits(hits)]
  if (!is.null(sizemargin) && !is.na(sizemargin)) {
    # take into account confidence intervals when calculating event size
    callwidth <- .distance(query, pquery)
    truthwidth <- .distance(subject, squery)
    callsize <- callwidth + .replaceNa(query$insLen, 0)
    truthsize <- truthwidth + .replaceNa(subject$insLen, 0)
    sizeerror <- .distance(
      IRanges::IRanges(start=callsize$min[S4Vectors::queryHits(hits)], end=callsize$max[S4Vectors::queryHits(hits)]),
      IRanges::IRanges(start=truthsize$min[S4Vectors::subjectHits(hits)], end=truthsize$max[S4Vectors::subjectHits(hits)])
    )$min
    # event sizes must be within sizemargin
    hits <- hits[sizeerror - 1 < sizemargin * pmin(callsize$max[S4Vectors::queryHits(hits)], truthsize$max[S4Vectors::subjectHits(hits)]),]
    # further restrict breakpoint positions for small events
    localbperror <- .distance(query[S4Vectors::queryHits(hits)], subject[S4Vectors::subjectHits(hits)])$min
    remotebperror <- .distance(pquery[S4Vectors::queryHits(hits)], squery[S4Vectors::subjectHits(hits)])$min
    if (!is.null(restrictMarginToSizeMultiple)) {
      allowablePositionError <- (pmin(callsize$max[S4Vectors::queryHits(hits)], truthsize$max[S4Vectors::subjectHits(hits)]) * restrictMarginToSizeMultiple + 1)
      hits <- hits[localbperror <= allowablePositionError & remotebperror <= allowablePositionError, ]
    }
  }
  return(hits)
}
# TODO: new function to annotate a Hits object with sizeerror, localbperror, and remotebperror
#' @noRd
.distance <- function(r1, r2) {
  return(data.frame(
    min=pmax(0, pmax(start(r1), start(r2)) - pmin(end(r1), end(r2))),
    max=pmax(end(r2) - start(r1), end(r1) - start(r2))))
}
#' Counting overlapping breakpoints between two breakpoint sets
#'
#' @details
#' \code{countBreakpointOverlaps()} returns the number of overlaps between breakpoint
#' objects, based on the output of \code{findBreakpointOverlaps()}.
#' See GenomicRanges::countOverlaps-methods
#' @param querygr,subjectgr,maxgap,minoverlap,ignore.strand,sizemargin,restrictMarginToSizeMultiple
#' See \code{findBreakpointOverlaps()}.
#' @param countOnlyBest Default value set to FALSE. When set to TRUE, the result count
#' each subject breakpoint as overlaping only the best overlapping query breakpoint.
#' The best breakpoint is considered to be the one with the highest QUAL score.
#' @param breakpointScoreColumn Query column defining a score for determining which query breakpoint
#' is considered the best when countOnlyBest=TRUE.
#' @examples
#' truth_vcf = VariantAnnotation::readVcf(system.file("extdata", "na12878_chr22_Sudmunt2015.vcf", 
#' package = "StructuralVariantAnnotation"))
#' crest_vcf = VariantAnnotation::readVcf(system.file("extdata", "na12878_chr22_crest.vcf", 
#' package = "StructuralVariantAnnotation"))
#' caller_bpgr = breakpointRanges(crest_vcf)
#' caller_bpgr$true_positive = countBreakpointOverlaps(caller_bpgr, breakpointRanges(truth_vcf),
#'   maxgap=100, sizemargin=0.25, restrictMarginToSizeMultiple=0.5, countOnlyBest=TRUE)
#' @return An integer vector containing the tabulated query overlap hits.
#' @export
countBreakpointOverlaps <- function(querygr, subjectgr, countOnlyBest=FALSE,
                                    breakpointScoreColumn = "QUAL", maxgap=-1L,
                                    minoverlap=0L, ignore.strand=FALSE, sizemargin=NULL,
                                    restrictMarginToSizeMultiple=NULL) {
  hitscounts <- rep(0, length(querygr))
  hits <- as.data.frame(findBreakpointOverlaps(querygr, subjectgr, maxgap, minoverlap, ignore.strand, sizemargin=sizemargin, restrictMarginToSizeMultiple=restrictMarginToSizeMultiple))
  if (!countOnlyBest) {
    hits <- hits %>%
      dplyr::group_by(.data$queryHits) %>%
      dplyr::summarise(n=dplyr::n())
  } else {
    # assign supporting evidence to the call with the highest QUAL
    hits$QUAL <- S4Vectors::mcols(querygr)[[breakpointScoreColumn]][hits$queryHits]
    hits <- hits %>%
      dplyr::arrange(dplyr::desc(.data$QUAL), .data$queryHits) %>%
      dplyr::distinct(.data$subjectHits, .keep_all=TRUE) %>%
      dplyr::group_by(.data$queryHits) %>%
      dplyr::summarise(n=dplyr::n())
  }
  hitscounts[hits$queryHits] <- hits$n
  return(hitscounts)
}

#' Converts a breakpoint GRanges object to a Pairs object
#' @param bpgr breakpoint GRanges object
#' @param writeQualAsScore write the breakpoint GRanges QUAL field as the score
#' fields for compatibility with BEDPE rtracklayer export
#' @param writeName write the breakpoint GRanges QUAL field as the score
#' fields for compatibility with BEDPE rtracklayer export
#' @param bedpeName function that returns the name to use for the breakpoint.
#' Defaults to the sourceId, name column, or row names (in that priority) of
#' the first breakend of each pair.
#' @param firstInPair function that returns TRUE for breakends that are considered
#' the first in the pair, and FALSE for the second in pair breakend. By default,
#' the first in the pair is the breakend with the lower ordinal in the breakpoint
#' GRanges object.
#' @examples
#' vcf.file <- system.file("extdata", "gridss.vcf", package = "StructuralVariantAnnotation")
#' bpgr <- breakpointRanges(VariantAnnotation::readVcf(vcf.file))
#' pairgr <- breakpointgr2pairs(bpgr)
#' #rtracklayer::export(pairgr, con="example.bedpe")
#' @return Pairs GRanges object suitable for export to BEDPE by rtracklayer
#' @rdname pairs2breakpointgr
#' @export
breakpointgr2pairs <- function(
  bpgr,
  writeQualAsScore=TRUE,
  writeName=TRUE,
  bedpeName = NULL,
  firstInPair = NULL) {
  .assertValidBreakpointGRanges(bpgr, "Cannot convert breakpoint GRanges to Pairs: ", allowSingleBreakends=FALSE)
  
  if (is.null(bedpeName)) {
    bedpeName = function(gr) { .replaceNull(.replaceNull(gr$sourceId, gr$name), names(gr)) }
  }
  if (is.null(firstInPair)) {
    firstInPair = function(gr) { seq_along(gr) < match(gr$partner, names(gr)) }
  }
  isFirst = firstInPair(bpgr)
  pairgr = S4Vectors::Pairs(bpgr[isFirst], partner(bpgr)[isFirst])
  if (writeName) {
    S4Vectors::mcols(pairgr)$name = bedpeName(S4Vectors::first(pairgr))
  }
  if (writeQualAsScore) {
    S4Vectors::mcols(pairgr)$score = S4Vectors::first(pairgr)$QUAL
  }
  return(pairgr)
}
#' @noRd
.assertValidBreakpointGRanges <- function(bpgr, friendlyErrorMessage="", allowSingleBreakends=TRUE) {
  if (is.null(names(bpgr))) {
    stop(paste0(friendlyErrorMessage, "Breakpoint GRanges require names"))
  }
  if (any(is.na(names(bpgr)))) {
    stop(paste0(friendlyErrorMessage, "Breakpoint GRanges names cannot be NA"))
  }
  if (any(duplicated(names(bpgr)))) {
    stop(paste0(friendlyErrorMessage, "Breakpoint GRanges names cannot duplicated"))
  }
  if (!allowSingleBreakends & any(is.na(bpgr$partner))) {
    stop(paste0(friendlyErrorMessage, "Breakpoint GRanges contains single breakends"))
  }
  if (any(duplicated(bpgr$partner) & !is.na(bpgr$partner))) {
    stop(paste0(friendlyErrorMessage,
                "Multiple breakends with the sample partner identified. ",
                "Breakends with multiple partners not currently supported by Breakpoint GRanges."))
  }
  else if (!all(is.na(bpgr$partner) | (bpgr$partner %in% names(bpgr) & names(bpgr) %in% bpgr$partner))) {
    stop(paste0(friendlyErrorMessage,
                "Unpartnered breakpoint found. ",
                "All breakpoints must contain a partner in the breakpoint GRanges."))
  }
}
#' Converts a BEDPE Pairs containing pairs of GRanges loaded using to a breakpoint GRanges object.
#' @details
#' Breakpoint-level column names will override breakend-level column names.
#' @param pairs a Pairs object consisting of two parallel genomic loci.
#' @param placeholderName prefix to use to ensure each entry has a unique ID.
#' @param firstSuffix first in pair name suffix to ensure breakend name uniqueness
#' @param secondSuffix second in pair name suffix to ensure breakend name uniqueness
#' @param nameField Fallback field for row names if the Pairs object does not contain any names.
#' BEDPE files loaded using rtracklayer use the "name" field.
#' @param renameScoreToQUAL renames the 'score' column to 'QUAL'.
#' Performing this rename results in a consistent variant quality score column
#' name for variant loaded from BEDPE and VCF.
#' @examples
#' bedpe.file <- system.file("extdata", "gridss.bedpe", package = "StructuralVariantAnnotation")
#' bedpe.pairs <- rtracklayer::import(bedpe.file)
#' bedpe.bpgr <- pairs2breakpointgr(bedpe.pairs)
#' @return Breakpoint GRanges object.
#' @export
pairs2breakpointgr <- function(
		pairs,
		placeholderName="bedpe",
		firstSuffix="_1", secondSuffix="_2",
		nameField="name",
		renameScoreToQUAL=TRUE) {
	n <- names(pairs)
	if (is.null(n)) {
		# BEDPE uses the "name" field
		if (nameField %in% names(S4Vectors::mcols(pairs))) {
			n <- S4Vectors::mcols(pairs)[[nameField]]
			mcols(pairs)$sourceId <- n
		} else {
			n <- rep(NA_character_, length(pairs))
		}
	}
	# ensure row names are unique
	n <- ifelse(is.na(n) | n == "" | n =="." | duplicated(n), paste0(placeholderName, seq_along(n)), n)
	#
	gr <- c(S4Vectors::first(pairs), S4Vectors::second(pairs))
	names(gr) <- c(paste0(n, firstSuffix), paste0(n, secondSuffix))
	gr$partner <- c(paste0(n, secondSuffix), paste0(n, firstSuffix))
	for (col in names(S4Vectors::mcols(pairs))) {
		if (col %in% nameField) {
			# drop columns we have processed
		} else {
			S4Vectors::mcols(gr)[[col]] <- S4Vectors::mcols(pairs)[[col]]
		}
	}
	if (renameScoreToQUAL) {
		names(mcols(gr))[which(names(mcols(gr)) == "score")] <- "QUAL"
		}
	return(gr)
}

#' Extracts the breakpoint sequence.
#'
#' @details
#' The sequence is the sequenced traversed from the reference anchor bases
#' to the breakpoint. For backward (-) breakpoints, this corresponds to the
#' reverse compliment of the reference sequence bases.
#'
#' @param gr breakpoint GRanges
#' @param ref Reference BSgenome
#' @param anchoredBases Number of bases leading into breakpoint to extract
#' @param remoteBases Number of bases from other side of breakpoint to extract
#' @return Breakpoint sequence around the variant position.
#' @export
extractBreakpointSequence <- function(gr, ref, anchoredBases, remoteBases=anchoredBases) {
	localSeq <- extractReferenceSequence(gr, ref, anchoredBases, 0)
	insSeq <- ifelse(strand(gr) == "-",
					 as.character(Biostrings::reverseComplement(DNAStringSet(.replaceNa(gr$insSeq,"")))),
					 .replaceNa(gr$insSeq, ""))
	remoteSeq <- as.character(Biostrings::reverseComplement(DNAStringSet(
		extractReferenceSequence(partner(gr), ref, remoteBases, 0))))
	return(paste0(localSeq, insSeq, remoteSeq))
}
#' Returns the reference sequence around the breakpoint position
#'
#' @details
#' The sequence is the sequenced traversed from the reference anchor bases
#' to the breakpoint. For backward (-) breakpoints, this corresponds to the
#' reverse compliment of the reference sequence bases.
#'
#' @param gr breakpoint GRanges
#' @param ref Reference BSgenome
#' @param anchoredBases Number of bases leading into breakpoint to extract
#' @param followingBases Number of reference bases past breakpoint to extract
#' @return Reference sequence around the breakpoint position.
#' @export
extractReferenceSequence <- function(gr, ref, anchoredBases, followingBases=anchoredBases) {
	assertthat::assert_that(is(gr, "GRanges"))
	assertthat::assert_that(is(ref, "BSgenome"))
	gr <- .constrict(gr)
	seqgr <- GRanges(seqnames=GenomeInfoDb::seqnames(gr), ranges=IRanges::IRanges(
		start=start(gr) - ifelse(strand(gr) == "-", followingBases, anchoredBases - 1),
		end=end(gr) + ifelse(strand(gr) == "-", anchoredBases - 1, followingBases)))
	startPad <- pmax(0, 1 - start(seqgr))
	endPad <- pmax(0, end(seqgr) - GenomeInfoDb::seqlengths(ref)[as.character(GenomeInfoDb::seqnames(seqgr))])
	GenomicRanges::ranges(seqgr) <- IRanges::IRanges(start=start(seqgr) + startPad, end=end(seqgr) - endPad)
	seq <- Biostrings::getSeq(ref, seqgr)
	seq <- paste0(stringr::str_pad("", startPad, pad="N"), as.character(seq), stringr::str_pad("", endPad, pad="N"))
	# DNAStringSet doesn't like out of bounds subsetting
	seq <- ifelse(strand(gr) == "-", as.character(Biostrings::reverseComplement(DNAStringSet(seq))), seq)
	return(seq)
}
#' constrict
#' @param gr GRanges object
#' @param ref reference 
#' @param position only 'middle' position is accepted.
#' @return A constricted GRanges object.
#' @noRd
.constrict <- function(gr, ref=NULL,position="middle") {
	isLower <- start(gr) < start(partner(gr))
	# Want to call a valid breakpoint
	#  123 456
	#
	#  =>   <= + -
	#  >   <== f f
	#
	#  =>  =>  + +
	#  >   ==> f c
	roundDown <- isLower | strand(gr) == "-"
	if (position == "middle") {
		pos <- (start(gr) + end(gr)) / 2
		GenomicRanges::ranges(gr) <- IRanges::IRanges(
			start=ifelse(roundDown,floor(pos), ceiling(pos)),
			width=1, names=names(gr))

	} else {
		stop(paste("Unrecognised position", position))
	}
	if (!is.null(ref)) {
		GenomicRanges::ranges(gr) <- IRanges::IRanges(start=pmin(pmax(1, start(gr)), GenomeInfoDb::seqlengths(ref)[as.character(GenomeInfoDb::seqnames(gr))]), width=1)
	}
	return(gr)
}

#' Calculates the length of inexact homology between the breakpoint sequence
#' and the reference
#'
#' @param gr reakpoint GRanges
#' @param ref reference BSgenome
#' @param anchorLength Number of bases to consider for homology
#' @param margin Number of additional reference bases include. This allows
#'		for inexact homology to be detected even in the presence of indels.
#' @param mismatch see Biostrings::pairwiseAlignment
#' @param gapOpening see Biostrings::pairwiseAlignment
#' @param gapExtension see Biostrings::pairwiseAlignment
#' @param match see Biostrings::pairwiseAlignment
#' @return A dataframe containing the length of inexact homology between the 
#' breakpoint sequence and the reference.
#' @export
calculateReferenceHomology <- function(gr, ref,
									   anchorLength=300,
									   margin=5,
									   match=2, mismatch=-6, gapOpening=5, gapExtension=3 # bwa
									   #match = 1, mismatch = -4, gapOpening = 6, gapExtension = 1, # bowtie2
) {
	# shrink anchor for small events to prevent spanning alignment
	aLength <- .replaceNa(pmin(anchorLength, abs(gr$svLen) + 1), anchorLength)
	anchorSeq <- extractReferenceSequence(gr, ref, aLength, 0)
	anchorSeq <- sub(".*N", "", anchorSeq)
	# shrink anchor with Ns
	aLength <- nchar(anchorSeq)
	varseq <- extractBreakpointSequence(gr, ref, aLength)
	varseq <- sub("N.*", "", varseq)
	bpLength <- nchar(varseq) - aLength
	nonbpseq <- extractReferenceSequence(gr, ref, 0, bpLength + margin)
	nonbpseq <- sub("N.*", "", nonbpseq)
	refseq <- paste0(anchorSeq, nonbpseq)

	partnerIndex <- match(gr$partner, names(gr))

	if (all(refseq=="") && all(varseq=="")) {
		# Workaround of Biostrings::pairwiseAlignment bug
		return(data.frame(
			exacthomlen=rep(NA, length(gr)),
			inexacthomlen=rep(NA, length(gr)),
			inexactscore=rep(NA, length(gr))))
	}

	aln <- Biostrings::pairwiseAlignment(varseq, refseq, type="local",
										 substitutionMatrix=nucleotideSubstitutionMatrix(match, mismatch, FALSE, "DNA"),
										 gapOpening=gapOpening, gapExtension=gapExtension, scoreOnly=FALSE)
	ihomlen <- Biostrings::nchar(aln) - aLength - deletion(nindel(aln))[,2] - insertion(nindel(aln))[,2]
	ibphomlen <- ihomlen + ihomlen[partnerIndex]
	ibpscore <- score(aln) + score(aln)[partnerIndex] - 2 * aLength * match

	# TODO: replace this with an efficient longest common substring function
	# instead of S/W with a massive mismatch/gap penalty
	penalty <- anchorLength * match
	matchLength <- Biostrings::pairwiseAlignment(varseq, refseq, type="local",
												 substitutionMatrix=nucleotideSubstitutionMatrix(match, -penalty, FALSE, "DNA"),
												 gapOpening=penalty, gapExtension=0, scoreOnly=TRUE) / match
	ehomlen <- matchLength - aLength
	ebphomlen <- ehomlen + ehomlen[partnerIndex]

	ebphomlen[aLength == 0] <- NA
	ibphomlen[aLength == 0] <- NA
	ibpscore[aLength == 0] <- NA
	return(data.frame(
		exacthomlen=ebphomlen,
		inexacthomlen=ibphomlen,
		inexactscore=ibpscore))
}


#' Converts to breakend notation
#' @param gr GRanges object.
#' @param insSeq insert sequence of the GRanges.
#' @param ref reference sequence of the GRanges.
#' @return breakendAlt or breakpointAlt depending on whether the variant is partnered.
#' @noRd
.toVcfBreakendNotationAlt = function(gr, insSeq=gr$insSeq, ref=gr$REF) {
	assertthat::assert_that(all(width(gr) == 1))
	assertthat::assert_that(!is.null(insSeq))
	assertthat::assert_that(all(insSeq != ""))
	assertthat::assert_that(!is.null(gr$partner))
	isBreakpoint = !is.na(gr$partner)
	breakendAlt = ifelse(as.character(strand(gr)) == "+", paste0(gr$insSeq, "."), paste0(".", gr$insSeq))
	gr$partner[isBreakpoint] = names(gr)[isBreakpoint] # self partner to prevent errors
	partnergr = gr[gr$partner]
	partnerDirectionChar = ifelse(strand(partnergr) == "+", "]", "[")
	breakpointAlt = ifelse(as.character(strand(gr)) == "+",
						   paste0(ref, insSeq, partnerDirectionChar, GenomeInfoDb::seqnames(partnergr), ":", start(partnergr), partnerDirectionChar),
						   paste0(partnerDirectionChar, GenomeInfoDb::seqnames(partnergr), ":", start(partnergr), partnerDirectionChar, insSeq, ref))
	return (ifelse(isBreakpoint, breakpointAlt, breakendAlt))
}

#' Converts the given breakpoint GRanges object to VCF format in breakend
#' notation.
#'
#' @param gr breakpoint GRanges object. Can contain both breakpoint and single 
#' breakend SV records.
#' @param ... For cbind and rbind a list of VCF objects. For all other methods 
#' ... are additional arguments passed to methods. See VCF class in 
#' VariantAnnotation for more details.
#' @return A VCF object.
breakpointGRangesToVCF <- function(gr, ...) {
	if (is.null(gr$insSeq)) {
		gr$insSeq = rep("", length(gr))
	}
	nominalgr = GRanges(seqnames=GenomeInfoDb::seqnames(gr), 
	                    ranges=IRanges::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)
	
	VariantAnnotation::VCF(rowRanges = GRanges(), colData = S4Vectors::DataFrame(), 
	    exptData = list(header = VCFHeader()), fixed = S4Vectors::DataFrame(), 
	    info = S4Vectors::DataFrame(), geno = S4Vectors::SimpleList(), ..., collapsed=FALSE, 
	    verbose = FALSE)

}
#' Type of simplest explanation of event. Possible types are:
#' | Type | Description |
#' | BND | Single breakend |
#' | CTX | Interchromosomal translocation |
#' | INV | Inversion. |
#' | DUP | Tandem duplication |
#' | INS | Insertion |
#' | DEL | Deletion |
#' 
#' Note that both ++ and -- breakpoint will be classified as inversions regardless of whether both breakpoint that consistitute an actual inversion exists or not
#' 
#' @param gr breakpoint GRanges object
#' @param insertionLengthThreshold portion of inserted bases compared to total event size to be classified as an insertion.
#' For example, a 5bp deletion with 5 inserted bases will be classified as an INS event.
#' @return Type of simplest explanation of event
#' @export
simpleEventType <- function(gr, insertionLengthThreshold=0.5) {
	if (is.null(gr$partner)) {
		gr$partner = rep(NA_character_, length(gr))
	}
	pgr = partner(gr, selfPartnerSingleBreakends=TRUE)
	return(
		ifelse(is.na(gr$partner), "BND", 
			ifelse(seqnames(gr) != seqnames(pgr), "CTX", # inter-chromosomosal
				ifelse(strand(gr) == strand(pgr), "INV",
					ifelse(gr$insLen >= abs(simpleEventLength(gr)) * insertionLengthThreshold, "INS", # TODO: improve classification of complex events
						ifelse(xor(start(gr) < start(pgr), strand(gr) == "-"), "DEL",
							"DUP"))))))
}
#' Length of event if interpreted as an isolated breakpoint.
#' @param gr breakpoint GRanges object
#' @return Length of the simplest explanation of this breakpoint/breakend.
#' @export
simpleEventLength <- function(gr) {
	if (is.null(gr$partner)) {
		gr$partner = rep(NA_character_, length(gr))
	}
	pgr = partner(gr, selfPartnerSingleBreakends=TRUE)
	return(
		ifelse(seqnames(gr) != seqnames(pgr) | as.logical(strand(gr) == strand(pgr) | is.na(gr$partner)), NA_integer_,
			gr$insLen + 1 + ifelse(as.logical(strand(gr) == "+"), start(gr) - start(pgr), start(pgr) - start(gr))))
}
#' Finds duplication events that are reported as inserts.
#' As sequence alignment algorithms do no allow backtracking, long read-based
#' variant callers will frequently report small duplication as insertion events.
#' Whilst both the duplication and insertion representations result in the same
#' sequence, this representational difference is problematic when comparing
#' variant call sets.
#' 
#' WARNING: this method does not check that the inserted sequence actually matched the duplicated sequence.
#' @param query a breakpoint GRanges object
#' @param subject a breakpoint GRanges object
#' @param maxgap maximum distance between the insertion position and the duplication
#' @param maxsizedifference maximum size difference between the duplication and insertion.
#' @return Hits object containing the ordinals of the matching breakends
#' in the query and subject 
#' @export
findInsDupOverlaps <- function(query, subject, maxgap=-1L, maxsizedifference=0L) {
	.assertValidBreakpointGRanges(query)
	.assertValidBreakpointGRanges(subject)
	query$ordinal = seq_len(length(query))
	subject$ordinal = seq_len(length(subject))
	query$set = simpleEventType(query)
	query$sel = simpleEventLength(query)
	subject$set = simpleEventType(subject)
	subject$sel = simpleEventLength(subject)
	pquery = partner(query)
	psubject = partner(subject)
	query$isLowBreakend = start(query) < start(pquery) | (start(query) == start(pquery) & query$ordinal < pquery$ordinal)
	subject$isLowBreakend = start(subject) < start(psubject) | (start(subject) == start(psubject) & subject$ordinal < psubject$ordinal)
	
	qins_to_sdup = .findOverlaps_queryIns_subjectDup(query, subject, psubject, maxgap=maxgap, maxsizedifference=maxsizedifference)
	sins_to_qdup = .findOverlaps_queryIns_subjectDup(subject, query, pquery, maxgap=maxgap, maxsizedifference=maxsizedifference)
	lowhits = data.frame(
		qhits=c(qins_to_sdup$queryHits, sins_to_qdup$subjectHits),
		shits=c(qins_to_sdup$subjectHits, sins_to_qdup$queryHits))
	# add upper to upper match
	bothhits = S4Vectors::Hits(
		from=c(lowhits$qhits, pquery$ordinal[lowhits$qhits]),
		to=c(lowhits$shits, psubject$ordinal[lowhits$shits]),
		nLnode=length(query),
		nRnode=length(subject))
	return(bothhits)
}
#' @noRd
.findOverlaps_queryIns_subjectDup <- function(query, subject, psubject , maxgap=-1L, maxsizedifference=0L) {
	subject$HighEndPosition = end(psubject)
	subject = subject[subject$set == "DUP" & subject$isLowBreakend]
	end(subject) = subject$HighEndPosition
	# expand by one since insertion can preceed, succeed, or be in the middle of the dup
	start(subject) = start(subject) - 1
	query = query[query$set == "INS" & query$isLowBreakend]
	hits = findOverlaps(query, subject, maxgap=maxgap, ignore.strand=TRUE)
	hits = hits[abs(query$sel[S4Vectors::queryHits(hits)] - subject$sel[S4Vectors::subjectHits(hits)]) <= maxsizedifference]
	# TODO: filter by 
	# Translate back to ordinals of what was passed in to us
	return(data.frame(
		queryHits=query$ordinal[S4Vectors::queryHits(hits)],
		subjectHits=subject$ordinal[S4Vectors::subjectHits(hits)]))
}
#' Identifies potential transitive imprecise calls that can be explained by
#' traversing multiple breakpoints.
#' 
#' Transitive calls are imprecise breakpoints or breakpoints with inserted sequence
#' that can be explained by a sequence of breakpoints.
#' That is, A-C calls in which additional sequence may be between A and C that
#' can be explained by A-B-C.
#' 
#' @param transitiveGr a breakpoint GRanges object containing imprecise calls
#' @param subjectGr breakpoints to traverse
#' @param maximumImpreciseInsertSize
#' Expected number of bases to traverse imprecise calls.
#' @param minimumTraversedBreakpoints
#' Minimum number of traversed breakpoints to consider a transitive
#' @param maximumTraversedBreakpoints
#' Maximum number of breakpoints to traverse when looking for an explanation of the transitive calls
#' @param positionalMargin 
#' Allowable margin of error when matching call positional overlaps.
#' A non-zero margin allows for matching of breakpoint with imperfect homology.
#' @param insertionLengthMargin
#' Allowable difference in length between the inserted sequence and the traversed
#' path length.
#' Defaults to 50bp to allow for long read indel errors.
#' @param insLen
#' Integer vector of same length as `transitiveGr` indicating the number
#' of bases inserted at the breakpoint.
#' 
#' Defaults to transitiveGr$insLen which will be present if the GRanges
#' was loaded from a VCF using breakpointRanges()
#' @param impreciseTransitiveCalls
#' Boolean vector of same length as `transitiveGr` indicating which calls
#' are imprecise calls. Defaults to calls with a non-zero interval size
#' that have no homology.
#' @param impreciseSubjectCalls
#' Boolean vector of same length as `subjectGr` indicating which calls
#' are imprecise calls. Defaults to calls with a non-zero interval size
#' that have no homology.
#' @param allowImprecise Allow traversal of imprecise calls.
#' Defaults to FALSE as to prevent spurious results which skip
#' some breakpoints when traversing multiple breakpoints
#' E.g. An A-D transitive from an underlying A-B-C-D rearrangement
#' will include A-B-D and A-C-D results if allowImprecise=TRUE.
#' @return `DataFrame` containing the transitive calls traversed with the following columns:
#' | column | meaning | 
#' | ------ | ------- | 
#' | transitive_breakpoint_name | Name of the transitive breakpoint a path was found for |
#' | total_distance | Total length (in bp) of the path |
#' | traversed_breakpoint_names | `CharacterList` of names of breakpoint traversed in the path |
#' | distance_to_traversed_breakpoint | `IntegerList` of distances from start of path to end of traversing breakpoint |
#' @export
findTransitiveCalls <- function(
		transitiveGr,
		subjectGr,
		maximumImpreciseInsertSize=700,
		minimumTraversedBreakpoints=2,
		maximumTraversedBreakpoints=6,
		positionalMargin=8,
		insertionLengthMargin=50,
		insLen=transitiveGr$insLen,
		impreciseTransitiveCalls=(transitiveGr$HOMLEN == 0 | is.null(transitiveGr$HOMLEN)) & start(transitiveGr) != end(transitiveGr),
		impreciseSubjectCalls=(subjectGr$HOMLEN == 0 | is.null(subjectGr$HOMLEN)) & start(subjectGr) != end(subjectGr),
		allowImprecise=FALSE) {
	if (is.null(insLen)) {
		stop("Missing insLen")
	}
	transitiveGr$.isImprecise = impreciseTransitiveCalls
	transitiveGr$insLen = insLen
	transitiveGr = transitiveGr[hasPartner(transitiveGr)]
	transitiveGr$.isImprecise = transitiveGr$.isImprecise | partner(transitiveGr)$.isImprecise
	transitiveGr$minimumTransitiveLength = ifelse(insLen > 0, pmax(0, insLen - insertionLengthMargin), 0)
	transitiveGr$maximumTransitiveLength = ifelse(insLen > 0, insLen + insertionLengthMargin, ifelse(impreciseTransitiveCalls, maximumImpreciseInsertSize, 0))
	transitiveGr = transitiveGr[transitiveGr$maximumTransitiveLength > 0]
	transitiveGr = transitiveGr[hasPartner(transitiveGr)]
	transitiveGr$ordinal = seq_len(length(transitiveGr))
	transitiveGr$partnerOrdinal = partner(transitiveGr)$ordinal
	if (!allowImprecise) {
		subjectGr = subjectGr[!impreciseSubjectCalls]
	}
	subjectGr = subjectGr[hasPartner(subjectGr)]
	# centre-align subject intervals to simplify the traversal logic
	start(subjectGr) = (start(subjectGr) + end(subjectGr)) / 2
	if (is.null(subjectGr$insLen)) {
		warning("insLen field missing. Assuming all traversed breakpoints have no inserted sequence")
		subjectGr$insLen = 0
	}
	subjectGr$insLen = .replaceNa(subjectGr$insLen, 0)
	subjectGr$ordinal = seq_len(length(subjectGr))
	subjectGr$partnerOrdinal = partner(subjectGr)$ordinal
	# transitive breakpoint must occur within the confidence interval bounds
	terminal_matches = as.data.frame(GenomicRanges::findOverlaps(transitiveGr, subjectGr, maxgap=positionalMargin, ignore.strand=FALSE)) %>%
		dplyr::select(
			terminalStartOrdinal=.data$queryHits,
			transitiveOrdinal=.data$subjectHits) %>%
		dplyr::mutate(
			terminalEndOrdinal=transitiveGr$partnerOrdinal[.data$terminalStartOrdinal])
	current_traversals = dplyr::inner_join(terminal_matches, terminal_matches, by=c("terminalStartOrdinal"="terminalEndOrdinal"), suffix=c(".start", ".end")) %>%
		dplyr::select(
			terminalStartOrdinal=.data$terminalStartOrdinal,
			currentTraverseInOrdinal=.data$transitiveOrdinal.start,
			endingTraverseOutOrdinal=.data$transitiveOrdinal.end,
			terminalEndOrdinal=.data$terminalEndOrdinal) %>%
		dplyr::mutate(
			currentTraverseOutOrdinal=subjectGr$partnerOrdinal[.data$currentTraverseInOrdinal],
			distance=subjectGr$insLen[.data$currentTraverseInOrdinal],
			# TAB is used as a placeholder as it's a disallowed character in VCF and causes a parsing error in BEDPE
			breakpointsTraversed=paste0(names(subjectGr)[.data$currentTraverseInOrdinal], "	"),
			traversedDistances=paste0(.data$distance, "	"),
			traversedBreakpoints=1,
			minimumTransitiveLength=transitiveGr$minimumTransitiveLength[.data$terminalStartOrdinal],
			maximumTransitiveLength=transitiveGr$maximumTransitiveLength[.data$terminalStartOrdinal])
	
	resultdf = data.frame(
		terminalStartOrdinal=integer(0),
		currentTraverseInOrdinal=integer(0),
		endingTraverseOutOrdinal=integer(0),
		terminalEndOrdinal=integer(0),
		currentTraverseOutOrdinal=integer(0),
		distance=integer(0),
		ordinalsTraversed=character(0),
		traversedDistances=character(0),
		traversedBreakpoints=integer(0))
	if (nrow(current_traversals > 0)) {
		traversable_segments = .traversable_segments(subjectGr, max(current_traversals$maximumTransitiveLength))
		# now we traverse
		while (nrow(current_traversals) > 0) {
			# Terminal paths
			current_traversals = current_traversals %>%
				dplyr::filter(.data$distance <= .data$maximumTransitiveLength & .data$traversedBreakpoints <= maximumTraversedBreakpoints) %>%
				dplyr::mutate(is_complete_path=.data$currentTraverseOutOrdinal == .data$endingTraverseOutOrdinal)
			resultdf = dplyr::bind_rows(
				resultdf,
				current_traversals %>% dplyr::filter(
					.data$is_complete_path &
					.data$distance >= .data$minimumTransitiveLength &
					.data$traversedBreakpoints >= minimumTraversedBreakpoints))
			# traverse
			current_traversals = current_traversals %>%
				dplyr::filter(!.data$is_complete_path) %>%
				dplyr::inner_join(traversable_segments, by=c("currentTraverseInOrdinal"="segmentStartExternalOrdinal")) %>%
				dplyr::mutate(
					currentTraverseInOrdinal=.data$segmentEndInternalOrdinal,
					currentTraverseOutOrdinal=.data$segmentEndExternalOrdinal,
					distance=.data$distance + .data$segmentLength + .data$segmentEndAdditionalLength,
					traversedBreakpoints=.data$traversedBreakpoints + 1,
					traversedDistances=paste0(.data$traversedDistances, .data$distance, "	"),
					breakpointsTraversed=paste0(.data$breakpointsTraversed, names(subjectGr)[.data$segmentEndInternalOrdinal], "	"))
			# How can drop columns without a "findTransitiveCalls: no visible binding for global variable 'segmentStartAdditionalLength'" NOTE in R check?
			#	dplyr::select(
			#		-segmentStartInternalOrdinal,
			#		-segmentLength,
			#		-segmentEndInternalOrdinal,
			#		-segmentEndExternalOrdinal,
			#		-segmentStartAdditionalLength,
			#		-segmentEndAdditionalLength)
			current_traversals$segmentStartInternalOrdinal = NULL
			current_traversals$segmentLength = NULL
			current_traversals$segmentEndInternalOrdinal = NULL
			current_traversals$segmentEndExternalOrdinal = NULL
			current_traversals$segmentStartAdditionalLength = NULL
			current_traversals$segmentEndAdditionalLength = NULL
		}
	}
	# transform results into long form
	resultdf = resultdf %>%
		dplyr::mutate(
			transitive_breakpoint_name=names(transitiveGr)[.data$terminalStartOrdinal],
			total_distance=as.integer(.data$distance),
			# trim trailing tab
			traversed_breakpoint_names=stringr::str_sub(.data$breakpointsTraversed, end=stringr::str_length(.data$breakpointsTraversed) - 1),
			distance_to_traversed_breakpoint=stringr::str_sub(.data$traversedDistances, end=stringr::str_length(.data$traversedDistances) - 1)) %>%
		dplyr::select(
			transitive_breakpoint_name=.data$transitive_breakpoint_name,
			total_distance=.data$total_distance,
			traversed_breakpoint_names=.data$traversed_breakpoint_names,
			distance_to_traversed_breakpoint=.data$distance_to_traversed_breakpoint) %>%
		as.data.frame() %>%
		S4Vectors::DataFrame()
	resultdf$traversed_breakpoint_names = IRanges::CharacterList(stringr::str_split(resultdf$traversed_breakpoint_names, stringr::fixed("	")))
	resultdf$distance_to_traversed_breakpoint = IRanges::IntegerList(stringr::str_split(resultdf$distance_to_traversed_breakpoint, stringr::fixed("	")))
	return(resultdf)
}
#' @noRd
.traversable_segments = function(gr, maxgap) {
	as.data.frame(GenomicRanges::findOverlaps(gr, gr, maxgap=maxgap, ignore.strand=TRUE)) %>%
		dplyr::filter(
			(as.logical(strand(gr)[.data$queryHits] == "-") & as.logical(strand(gr)[.data$subjectHits] == "+") & start(gr)[.data$queryHits] <= start(gr)[.data$subjectHits]) |
				(as.logical(strand(gr)[.data$queryHits] == "+") & as.logical(strand(gr)[.data$subjectHits] == "-") & start(gr)[.data$queryHits] >= start(gr)[.data$subjectHits])) %>%
		dplyr::select(
			segmentStartInternalOrdinal=.data$queryHits,
			segmentEndInternalOrdinal=.data$subjectHits) %>%
		dplyr::mutate(
			segmentStartExternalOrdinal=gr$partnerOrdinal[.data$segmentStartInternalOrdinal],
			segmentEndExternalOrdinal=gr$partnerOrdinal[.data$segmentEndInternalOrdinal],
			segmentLength=1 + abs(start(gr)[.data$segmentStartInternalOrdinal] - start(gr)[.data$segmentEndInternalOrdinal]),
			segmentStartAdditionalLength=gr$insLen[.data$segmentStartInternalOrdinal],
			segmentEndAdditionalLength=gr$insLen[.data$segmentEndInternalOrdinal]) %>%
		dplyr::select(
			segmentStartExternalOrdinal=.data$segmentStartExternalOrdinal,
			segmentStartInternalOrdinal=.data$segmentStartInternalOrdinal,
			segmentStartAdditionalLength=.data$segmentStartAdditionalLength,
			segmentLength=.data$segmentLength,
			segmentEndAdditionalLength=.data$segmentEndAdditionalLength,
			segmentEndInternalOrdinal=.data$segmentEndInternalOrdinal,
			segmentEndExternalOrdinal=.data$segmentEndExternalOrdinal)
}
d-cameron/StructuralVariantAnnotation documentation built on Dec. 26, 2021, 6:14 p.m.