R/disjoinPafAlignments.R

In daewoooo/SVbyEye: Visualization of genomic structural variants

#' Disjoin overlapping PAF alignments.
#'
#' This function takes loaded PAF alignments using \code{\link{readPaf}} function and identify overlapping
#' genomic ranges either in target or query coordinates and then split PAF alignments at the positions of
#' these overlaps into a disjoined set of genomic ranges. Alternatively, user defined set of genomic ranges
#' in parameter `disjoin.gr` can be used to disjoin PAF alignments.
#'
#' @param min.overlap A minimum number of overlapping base pairs to disjoin PAF alignments.
#' @param coordinates PAF coordinates to disjoin, either 'target' or 'query'.
#' @param disjoin.gr A \code{\link{GRanges-class}} object containing genomic ranges where PAF alignments will be disjoined
#' (must be in defined `coordinates`).
#' @inheritParams breakPaf
#' @return A \code{tibble} of disjoined PAF alignments.
#' @import GenomicRanges
#' @importFrom GenomicAlignments GAlignments cigarNarrow explodeCigarOpLengths
#' @importFrom GenomeInfoDb seqnames
#' @importFrom S4Vectors subjectHits queryHits
#' @importFrom IRanges IRanges
#' @importFrom methods is
#' @importFrom tibble tibble
#' @author David Porubsky
#' @export
#' @examples
#' ## Get PAF to plot
#' paf.file <- system.file("extdata", "test1.paf", package = "SVbyEye")
#' ## Read in PAF
#' paf.table <- readPaf(paf.file = paf.file, include.paf.tags = TRUE, restrict.paf.tags = "cg")
#' ## Disjoin PAF at target coordinates
#' disj.paf.table <- disjoinPafAlignments(paf.table = paf.table, coordinates = 'target')
#' ## Disjoin PAF at user defined coordinates
#' disj.gr <- GenomicRanges::GRanges(seqnames = 'query.region',
#'                                   ranges = IRanges::IRanges(start = 16300000, end = 16350000))
#' disj.paf.table <- disjoinPafAlignments(paf.table = paf.table,
#'                                        coordinates = 'query', disjoin.gr = disj.gr)
#'
disjoinPafAlignments <- function(paf.table, min.overlap = 1000, coordinates = 'target', disjoin.gr = NULL) {
  ptm <- startTimedMessage(paste0("[disjoinPafAlignments] Disjoining ", coordinates, " alignments"))

  ## Check user input ##
  ## Make sure PAF has at least 12 mandatory fields
  if (ncol(paf.table) >= 12) {
    paf <- paf.table
  } else {
    stop("Submitted PAF alignments do not contain a minimum of 12 mandatory fields, see PAF file format definition !!!")
  }

  ## Process target ranges ##
  if (coordinates == 'target') {
    target.gr <- GenomicRanges::makeGRangesFromDataFrame(paf, seqnames.field = "t.name", start.field = "t.start", end.field = "t.end")

    if (methods::is(disjoin.gr, "GRanges")) {
      hits <- GenomicRanges::findOverlaps(disjoin.gr, target.gr)
      if (length(hits) > 0) {
        todisj.gr <- disjoin.gr
      } else {
        message("No genomic ranges in 'disjoin.gr' overlaps target coordinates")
      }
    } else {
      ## Get regions covered more than once
      cov.gr <- as(GenomicRanges::coverage(target.gr), 'GRanges')
      cov.gr <- cov.gr[cov.gr$score > 1]
      ## Filter by minimum overlapping regions
      todisj.gr <- cov.gr[GenomicRanges::width(cov.gr) >= min.overlap, 0]
    }

    ## Get disjoined set of query and target ranges
    if (length(todisj.gr) > 0) {
      ## Remove regions to disjoin from target
      subtr.gr <- suppressWarnings( unlist(GenomicRanges::subtract(target.gr, todisj.gr, ignore.strand=TRUE)) )
      ## Get intersection between regions to disjoin and target
      inters.gr <- suppressWarnings( GenomicRanges::intersect(target.gr, todisj.gr, ignore.strand=TRUE) )
      target.gr <- suppressWarnings( sort(c(subtr.gr, inters.gr), ignore.strand=TRUE) )
      ## Lift target ranges to query coordinates
      query.gr <- suppressMessages( liftRangesToAlignment(paf.table = paf, gr = target.gr, direction = 'target2query', report.cigar.str = TRUE) )
      ## Expand target ranges to match mapped query ranges
      target.gr <- target.gr[query.gr$idx]
      target.gr$alignmentsHits <- query.gr$alignmentsHits
      GenomicRanges::strand(target.gr) <- GenomicRanges::strand(query.gr)
      cigars.region <- query.gr$cg
      processCIGAR <- TRUE
    } else {
      processCIGAR <- FALSE
    }
    ## Process query ranges ##
  } else if (coordinates == 'query') {
    query.gr <- GenomicRanges::makeGRangesFromDataFrame(paf, seqnames.field = "q.name", start.field = "q.start", end.field = "q.end")

    if (methods::is(disjoin.gr, "GRanges")) {
      hits <- GenomicRanges::findOverlaps(disjoin.gr, query.gr)
      if (length(hits) > 0) {
        todisj.gr <- disjoin.gr
      } else {
        message("No genomic ranges in 'disjoin.gr' overlaps target coordinates")
      }
    } else {
      ## Get regions covered more than once
      cov.gr <- as(GenomicRanges::coverage(query.gr), 'GRanges')
      cov.gr <- cov.gr[cov.gr$score > 1]
      ## Filter by minimum overlapping regions
      todisj.gr <- cov.gr[GenomicRanges::width(cov.gr) >= min.overlap, 0]
    }

    ## Get disjoined set of query and target ranges
    if (length(todisj.gr) > 0) {
      ## Remove regions to disjoin from target
      subtr.gr <- suppressWarnings( unlist(GenomicRanges::subtract(query.gr, todisj.gr, ignore.strand=TRUE)) )
      ## Get intersection between regions to disjoin and target
      inters.gr <- suppressWarnings( GenomicRanges::intersect(query.gr, todisj.gr, ignore.strand=TRUE) )
      query.gr <- suppressWarnings( sort(c(subtr.gr, inters.gr), ignore.strand=TRUE) )
      ## Lift query ranges to target coordinates
      target.gr <- suppressMessages( liftRangesToAlignment(paf.table = paf, gr = query.gr, direction = 'query2target', report.cigar.str = TRUE) )
      ## Expand query ranges to match mapped target ranges
      query.gr <- query.gr[target.gr$idx]
      GenomicRanges::strand(query.gr) <- GenomicRanges::strand(target.gr)
      cigars.region <- target.gr$cg
      processCIGAR <- TRUE
    } else {
      processCIGAR <- FALSE
    }
  }

  if (processCIGAR) {
  #   ## Subset the cigar string per target region ##
    # ## Define PAF alignment object
    # alignments <- GenomicAlignments::GAlignments(
    #   seqnames = paf$t.name,
    #   pos = as.integer(paf$t.start) + 1L,
    #   #pos = as.integer(paf$t.start),
    #   cigar = paf$cg,
    #   strand = GenomicRanges::strand(paf$strand),
    #   names = paf$t.name
    # )
    # ## Get overlaps between target ranges and alignments
    # #hits <- findOverlaps(target.gr, alignments, minoverlap = min.overlap)
    # hits <- GenomicRanges::findOverlaps(target.gr, alignments)
  #   end(target.gr[which.max(end(target.gr))]) <- max(end(alignments))
  #   hits <- GenomicRanges::findOverlaps(target.gr, alignments, type = 'within')
  #
  #   ## Keep only longest overlaps for each target range
  #   intersect.gr <- GenomicRanges::pintersect(target.gr[S4Vectors::queryHits(hits)], as(alignments[S4Vectors::subjectHits(hits)], 'GRanges'))
  #   overlap.width <- GenomicRanges::width(intersect.gr)
  #   max.overlap <- sapply(split(overlap.width, S4Vectors::queryHits(hits)), max)
  #   hits <- hits[overlap.width %in% max.overlap]
  #   #hits <- hits[which(unique(overlap.width) %in% max.overlap)]
  #
  #   ## Get local alignment coordinates [in parallel fashion]
  #   #lifted.gr <- GenomicAlignments::pmapToAlignments(x = target.gr, alignments = alignments[S4Vectors::subjectHits(hits)])
  #   #names(lifted.gr) <- NULL
  #
  #   ## Set target range to local alignment coordinates
  #   suppressWarnings( target.local.gr <- GenomicRanges::shift(target.gr, shift = -(start(alignments)[S4Vectors::subjectHits(hits)] - 1)) )
  #   suppressWarnings( target.local.gr <- GenomicRanges::shift(target.gr, shift = -(start(alignments)[target.gr$alignmentsHits] - 1)) )
  #   start(target.local.gr[start(target.local.gr) == 0]) <- 1
  #   #target.local.gr <- GenomicRanges::trim(target.local.gr)
  #
  #   ## Define start and end position for CIGAR extraction
  #   #starts <- pmax(GenomicRanges::start(lifted.gr), 1)
  #   #width <- GenomicRanges::width(target.gr) - 1
  #   #ends <- pmin(GenomicRanges::end(lifted.gr), GenomicRanges::width(alignments)[S4Vectors::subjectHits(hits)])
  #   starts <- GenomicRanges::start(target.local.gr)
  #   width <- GenomicRanges::width(target.local.gr)
  #
  #   cigars <- alignments@cigar[S4Vectors::subjectHits(hits)]
  #   cigars.region <- vapply(seq_along(starts), function(i) {
  #     tryCatch(
  #       {
  #         GenomicAlignments::cigarNarrow(cigar = cigars[i], start = starts[i], width = width[i])[1]
  #       },
  #       error = function(e) {
  #         return("1=")
  #       }
  #     )
  #   }, FUN.VALUE = character(1))

    ## Get alignment length from regional cigars
    aln.len <- vapply(cigars.region, function(cg) sum(GenomicAlignments::explodeCigarOpLengths(cigar = cg)[[1]]), USE.NAMES = FALSE, FUN.VALUE = numeric(1))
    ## Get matched bases from regional cigars
    n.match <- vapply(cigars.region, function(cg) sum(GenomicAlignments::explodeCigarOpLengths(cigar = cg, ops = c("=", "M"))[[1]]), USE.NAMES = FALSE, FUN.VALUE = numeric(1))
    ## Get single-base mismatches from regional cigars
    n.mismatch <- vapply(cigars.region, function(cg) sum(GenomicAlignments::explodeCigarOpLengths(cigar = cg, ops = c("X"))[[1]]), USE.NAMES = FALSE, FUN.VALUE = numeric(1))

    ## Prepare modified PAF for export ##
    # return.paf <- tibble::tibble("q.name" = as.character(GenomeInfoDb::seqnames(query.gr))
    #                           ,"q.len" = paf$q.len[S4Vectors::subjectHits(hits)]
    #                           ,"q.start" = GenomicRanges::start(query.gr)
    #                           ,"q.end" = GenomicRanges::end(query.gr)
    #                           ,"strand" = as.character( GenomicRanges::strand(query.gr) )
    #                           ,"t.name" = as.character(GenomeInfoDb::seqnames(target.gr))
    #                           ,"t.len" = paf$t.len[S4Vectors::subjectHits(hits)]
    #                           ,"t.start" = GenomicRanges::start(target.gr)
    #                           ,"t.end" = GenomicRanges::end(target.gr)
    #                           ,"n.match" = n.match
    #                           ,"aln.len" = aln.len
    #                           ,"mapq" = paf$mapq[S4Vectors::subjectHits(hits)]
    #                           ,"cg" = cigars.region)

    return.paf <- tibble::tibble("q.name" = as.character(GenomeInfoDb::seqnames(query.gr))
                                 ,"q.len" = paf$q.len[target.gr$alignmentsHits]
                                 ,"q.start" = GenomicRanges::start(query.gr)
                                 ,"q.end" = GenomicRanges::end(query.gr)
                                 ,"strand" = as.character( GenomicRanges::strand(query.gr) )
                                 ,"t.name" = as.character(GenomeInfoDb::seqnames(target.gr))
                                 ,"t.len" = paf$t.len[target.gr$alignmentsHits]
                                 ,"t.start" = GenomicRanges::start(target.gr)
                                 ,"t.end" = GenomicRanges::end(target.gr)
                                 ,"n.match" = n.match
                                 ,"aln.len" = aln.len
                                 ,"mapq" = paf$mapq[target.gr$alignmentsHits]
                                 ,"cg" = cigars.region)

  } else {
    return.paf <- paf
    message('    No genomic ranges to disjoin were found !!!')
  }

  stopTimedMessage(ptm)
  ## Return PAF
  return(return.paf)
}