fix.R

junction.support = function (reads, junctions = NULL, bwa = NULL, ref = NULL, pad = 500, 
    pad.ref = pad * 20, realign = TRUE, walks = NULL, verbose = TRUE, 
    ...) 
{
    if (!inherits(reads, "GRanges") || is.null(reads$qname) || 
        is.null(reads$cigar) || is.null(reads$seq) || is.null(reads$flag)) 
        stop("read input must be GRanges with fields $qname, $cigar, $seq, $flag and optionally $AS")
    if (!is.null(junctions)) 
        walks = jJ(junctions$grl)$gw(pad = pad)
    if (is.null(walks)) 
        stop("Either walks or junctions must be provided")
    if (!realign) {
        if (is.null(junctions)) 
            junctions = walks$edges$junctions
        reads = gr.flipstrand(reads)
        reads$R1 = bamUtils::bamflag(reads$flag)[, "isFirstMateRead"] > 
            0
        r1 = reads %Q% (R1 == TRUE) %>% as.data.table
        r2 = reads %Q% (R1 == FALSE) %>% as.data.table
        ov = merge(r1, r2, by = "qname")
        sl = seqlengths(reads)
        grl = grl.pivot(GRangesList(dt2gr(ov[, .(seqnames = seqnames.x, 
            start = start.x, end = end.x, strand = strand.x)], 
            seqlengths = sl), dt2gr(ov[, .(seqnames = seqnames.y, 
            start = start.y, end = end.y, strand = strand.y)], 
            seqlengths = sl)))
        values(grl)$qname = ov$qname
        jn = merge(jJ(grl), junctions, cartesian = TRUE, pad = pad)
        if (!length(jn)) 
            return(reads[c()])
        ## out = merge(as.data.table(reads), unique(jn$dt[, .(qname, 
        ##     junction.id = subject.id)]), by = "qname") %>% dt2gr(seqlengths = sl)
        out = dt2gr(setDT(dplyr::inner_join(as.data.frame(reads),
                                            as.data.frame(unique(jn$dt[, .(qname, junction.id = subject.id)])), by = "qname")), seqlengths = sl)
        return(out)
    }
    if (inherits(bwa, "character") && file.exists(bwa)) {
        if (verbose) 
            message("Loading BWA index")
        bwa = BWA(bwa)
    }
    if (!inherits(ref, "DNAStringSet")) {
        if (verbose) 
            message("Loading genome reference as DNAStringSet")
        ref = rtracklayer::import(bwa@reference)
    }
    names(ref) = strsplit(names(ref), "\\s+") %>% sapply("[", 
        1)
    if (length(setdiff(seqnames(walks$nodes$gr), seqlevels(ref)))) 
        stop("seqlevels mismatch between junctions / walks and reference, please check reference (e.g. chr issues)")
    if (length(setdiff(seqnames(walks$nodes$gr), seqlevels(bwa)))) 
        stop("seqlevels mismatch between junctions / walks and BWA reference, please check reference (e.g. chr issues)")
    if (verbose) 
        message("Building and mapping derivative contigs")
    contig = bwa[ref[gr.fix(walks$grl, ref, drop = TRUE)]]
    if (verbose) 
        message("Building reference contigs flanking junctions")
    contigref = ref[gr.fix(walks$edges$junctions$footprint + 
        pad.ref, ref, drop = TRUE)]
    if (verbose) 
        message("Making gChain mapping contigs to reference")
    cg.contig = gChain::cgChain(contig)
    if (verbose) 
        message("Running contig support")
    reads = contig.support(reads, contig, ref = contigref, cg.contig = cg.contig, 
        ...)
    reads$junction.id = reads$contig.id
    return(reads)
}

overwritefun("junction.support", "junction.support", "skitools")

junction.support = skitools::junction.support