R/ceScan.R

In CNEr: CNE Detection and Visualization

### -----------------------------------------------------------------
### The main function for scanning the axts and get the CNEs
### Not exported!
ceScanR <- function(axts, tFilter=NULL, qFilter=NULL, tSizes, qSizes, 
                    thresholds=c("49_50")){
  winSize <- as.integer(sapply(strsplit(thresholds, "_"), "[", 2))
  minScore <- as.integer(sapply(strsplit(thresholds, "_"), "[", 1))
  resFiles <- tempfile(pattern=paste(minScore, winSize, "ceScan", sep="-"), 
                       tmpdir=tempdir(), fileext="")
  ## How stupid I am...
  if(is.null(tFilter) && is.null(qFilter)){
    .Call2("myCeScan", NULL, NULL, NULL,
           NULL, NULL, NULL,
           as.character(seqnames(qSizes)), as.integer(seqlengths(qSizes)),
           as.character(seqnames(targetRanges(axts))),
           start(targetRanges(axts)), end(targetRanges(axts)),
           as.character(strand(targetRanges(axts))),
           as.character(targetSeqs(axts)),
           as.character(seqnames(queryRanges(axts))),
           start(queryRanges(axts)), end(queryRanges(axts)),
           as.character(strand(queryRanges(axts))),
           as.character(querySeqs(axts)),
           score(axts), symCount(axts), winSize, minScore,
           as.character(resFiles),
           PACKAGE="CNEr")
  }else if(is.null(tFilter) && !is.null(qFilter)){
    .Call2("myCeScan", NULL, NULL, NULL,
           as.character(seqnames(qFilter)), start(qFilter), end(qFilter),
           as.character(seqnames(qSizes)), as.integer(seqlengths(qSizes)),
           as.character(seqnames(targetRanges(axts))),
           start(targetRanges(axts)), end(targetRanges(axts)),
           as.character(strand(targetRanges(axts))),
           as.character(targetSeqs(axts)),
           as.character(seqnames(queryRanges(axts))),
           start(queryRanges(axts)), end(queryRanges(axts)),
           as.character(strand(queryRanges(axts))),
           as.character(querySeqs(axts)),
           score(axts), symCount(axts), winSize, minScore,
           as.character(resFiles),
           PACKAGE="CNEr")
  }else if(!is.null(tFilter) && is.null(qFilter)){
    .Call2("myCeScan", as.character(seqnames(tFilter)), 
           start(tFilter), end(tFilter),
           NULL, NULL, NULL,
           as.character(seqnames(qSizes)), as.integer(seqlengths(qSizes)),
           as.character(seqnames(targetRanges(axts))),
           start(targetRanges(axts)), end(targetRanges(axts)),
           as.character(strand(targetRanges(axts))),
           as.character(targetSeqs(axts)),
           as.character(seqnames(queryRanges(axts))),
           start(queryRanges(axts)), end(queryRanges(axts)),
           as.character(strand(queryRanges(axts))),
           as.character(querySeqs(axts)),
           score(axts), symCount(axts), winSize, minScore,
           as.character(resFiles),
           PACKAGE="CNEr")
  }else{
    .Call2("myCeScan", 
           as.character(seqnames(tFilter)), start(tFilter), end(tFilter),
           as.character(seqnames(qFilter)), start(qFilter), end(qFilter),
           as.character(seqnames(qSizes)), as.integer(seqlengths(qSizes)), 
           as.character(seqnames(targetRanges(axts))), 
           start(targetRanges(axts)), end(targetRanges(axts)), 
           as.character(strand(targetRanges(axts))), 
           as.character(targetSeqs(axts)),
           as.character(seqnames(queryRanges(axts))), 
           start(queryRanges(axts)), end(queryRanges(axts)), 
           as.character(strand(queryRanges(axts))), 
           as.character(querySeqs(axts)),
           score(axts), symCount(axts), winSize, minScore, 
           as.character(resFiles),
           PACKAGE="CNEr")
  }
  CNE <- lapply(resFiles, 
                function(x){
                  res <- try(read.table(x, header=FALSE, sep="\t", as.is=TRUE),
                             silent=TRUE)
                  if(class(res) == "try-error"){
                    return(GRangePairs(first=GRanges(seqinfo=tSizes), 
                                       second=GRanges(seqinfo=qSizes))
                           )
                  }
                  colnames(res) <- c("tName", "tStart", "tEnd", 
                                     "qName", "qStart", "qEnd",
                                     "strand", "score", "cigar")
                  ans <- GRangePairs(first=GRanges(seqnames=res$tName,
                                             ranges=IRanges(start=res$tStart,
                                                            end=res$tEnd),
                                             strand="+",
                                             seqinfo=tSizes),
                                     second=GRanges(seqnames=res$qName,
                                              ranges=IRanges(start=res$qStart,
                                                             end=res$qEnd),
                                              strand=res$strand,
                                              seqinfo=qSizes),
                                     score=res$score,
                                     cigar=res$cigar
                                     )
                  return(ans)
                  }
                )
  names(CNE) <-  paste(minScore, winSize, sep="_")
  unlink(resFiles)
  return(CNE)
}

### -----------------------------------------------------------------
### The function for ceScan: one way for axt object, two way for CNE
### Exported!
setGeneric("ceScan", function(x, tFilter=NULL, qFilter=NULL,
                              tSizes=NULL, qSizes=NULL,
                              window=50L, identity=50L)
  standardGeneric("ceScan"))

setMethod("ceScan", "Axt", function(x, tFilter=NULL, qFilter=NULL,
                                    tSizes=NULL, qSizes=NULL,
                                    window=50L, identity=50L){
  ceScanAxt(x, tFilter=tFilter, qFilter=qFilter, tSizes=tSizes, qSizes=qSizes, 
            window=window, identity=identity)
})

setMethod("ceScan", "CNE", function(x, tFilter=NULL, qFilter=NULL,
                                    tSizes=NULL, qSizes=NULL,
                                    window=50L, identity=50L){
  axt12 <- readAxt(x@axt12Fn)
  axt21 <- readAxt(x@axt21Fn)
  cne12 <- ceScan(axt12, tFilter, qFilter,
                  tSizes=seqinfo(TwoBitFile(x@assembly1Fn)),
                  qSizes=seqinfo(TwoBitFile(x@assembly2Fn)),
                  window=window, identity=identity)
  cne21 <- ceScan(axt21, qFilter, tFilter,
                  tSizes=seqinfo(TwoBitFile(x@assembly2Fn)),
                  qSizes=seqinfo(TwoBitFile(x@assembly1Fn)),
                  window=window, identity=identity)
  ans <- list()
  for(i in 1:length(cne12)){
    ans[[names(cne12)[i]]] <- BiocGenerics:::replaceSlots(x, 
                                                          CNE12=cne12[[i]],
                                                          CNE21=cne21[[i]],
                                                          window=window[i],
                                                          identity=identity[i])
  }
  return(ans)
})

ceScanAxt <- function(axts, tFilter=NULL, qFilter=NULL,
                      tSizes=NULL, qSizes=NULL,
                      window=50L, identity=50L){
  # Prepare tSizes
  if(is.null(tSizes)){
    ## tSizes is NULL
    tSizes <- seqinfo(targetRanges(axts))
    if(any(is.na(seqlengths(tSizes)))){
      stop("tSizes must be provided or seqinfo is available in Axt object!")
    }
  }else if(!is(tSizes, "Seqinfo")){
    ## tSizes is integer vector
    tSizes <- Seqinfo(seqnames=names(tSizes), seqlengths=tSizes)
  }
  ## Check tFilter, tFilter's seqnames must %in% tSizes
  if(!is.null(tFilter)){
    if(!is(tFilter, "GRanges")){
      stop("tFilter must be NULL or a GRanges object!")
    }else{
      if(!all(seqnames(tFilter) %in% seqnames(tSizes))){
        stop("All the chromosomes in tFilter must exist in tSizes!")
      }
    }
  }
  ## Check axts target chromosomes
  if(!all(seqnames(targetRanges(axts)) %in% seqnames(tSizes))){
    stop("All the chromosomes in targetRanges of Axt object must exist in tSizes!")
  }
  
  # Prepare qSizes
  if(is.null(qSizes)){
    ## qSizes is NULL
    qSizes <- seqinfo(queryRanges(axts))
    if(any(is.na(seqlengths(qSizes)))){
      stop("qSizes must be provided or seqinfo is available in Axt object!")
    }
  }else if(!is(qSizes, "Seqinfo")){
    ## qSizes is integer vector
    qSizes <- Seqinfo(seqnames=names(qSizes), seqlengths=qSizes)
  }
  ## Check qFilter, qFilter's seqnames must %in% qSizes
  if(!is.null(qFilter)){
    if(!is(qFilter, "GRanges")){
      stop("qFilter must be NULL or a GRanges object!")
    }else{
      if(!all(seqnames(qFilter) %in% seqnames(qSizes))){
        stop("All the chromosomes in qFilter must exist in qSizes!")
      }
    }
  }
  ## Check axts query chromosomes
  if(!all(seqnames(queryRanges(axts)) %in% seqnames(qSizes))){
    stop("All the chromosomes in queryRanges of Axt object must exist in qSizes!")
  }
  
  # Check identity and window size
  if(any(window < identity)){
    stop("The scanning window size must be equal or larger than identity!")
  }
  if(length(identity) != length(window)){
    stop("The length of identity and window are different!")
  }
  thresholds <- paste(identity, window, sep="_")
  ans <- ceScanR(axts, tFilter=tFilter, qFilter=qFilter, 
                 tSizes=tSizes, qSizes=qSizes, thresholds=thresholds)
}

### -----------------------------------------------------------------
### Merge two side cnes (GRangePairs object), mcols will be discarded.
### strand information is no longer important.
### Exported!
setGeneric("cneMerge", function(cne12, cne21) standardGeneric("cneMerge"))

setMethod("cneMerge", signature(cne12="CNE", cne21="missing"),
          function(cne12, cne21){
            cneMerged <- cneMergeGRangePairs(cne12@CNE12, cne12@CNE21)
            BiocGenerics:::replaceSlots(cne12, CNEMerged=cneMerged)
          })

setMethod("cneMerge", signature(cne12="GRangePairs", cne21="GRangePairs"),
          function(cne12, cne21){
            cneMergeGRangePairs(cne12, cne21)
          })

cneMergeGRangePairs <- function(cne12, cne21){
  if(!is(cne12, "GRangePairs") || !is(cne21, "GRangePairs")){
    stop("cne12 and cne21 must be a GRangePairs object!")
  }
  strand(cne12@first) <- strand(cne12@second) <- strand(cne21@first) <- 
    strand(cne21@second) <- "+"
  cne <- c(cne12, swap(cne21))
  
  # 1. using reduce: the problem: it won't deal with {1,2,3}, {1,2} case
  # firstReduce <- reduce(first(cne), with.revmap=TRUE)
  # lastReduce <- reduce(last(cne), with.revmap=TRUE)
  # overlapFirstLast <- IntegerList(intersect(firstReduce$revmap,
  #                                           lastReduce$revmap))
  # 
  # ## First deal with the merged ranges
  # overlapFirstLast <- overlapFirstLast[elementNROWS(overlapFirstLast) > 1L]
  # firstIndex <- match(paste(overlapFirstLast, collapse="-"),
  #                     paste(firstReduce$revmap, collapse="-"))
  # lastIndex <- match(paste(overlapFirstLast, collapse="-"),
  #                    paste(lastReduce$revmap, collapse="-"))
  # ansFirst <- firstReduce[firstIndex]
  # mcols(ansFirst) <- NULL
  # ansLast <- lastReduce[lastIndex]
  # mcols(ansLast) <- NULL
  # 
  # ## Then deal with the unmerged ranges
  # ansFirst <- c(ansFirst, first(cne)[-sort(unlist(overlapFirstLast))])
  # ansLast <- c(ansLast, last(cne)[-sort(unlist(overlapFirstLast))])
  # return(GRangePairs(first=ansFirst, last=ansLast))
  
  # 2. using the findOverlaps:
  firstHits <- findOverlaps(first(cne), type="within",
                            drop.self=TRUE, drop.redundant=TRUE)
  lastHist <- findOverlaps(second(cne), type="within",
                           drop.self=TRUE, drop.redundant=TRUE)
  redundance <- IRanges::intersect(firstHits, lastHist)
  if(length(redundance) == 0L){
    return(cne)
  }else{
    return(cne[-queryHits(redundance)])
  }
}

### -----------------------------------------------------------------
### blatCNE: 
### Exported!
blatCNE <- function(cne, blatOptions=NULL, cutIdentity=90){
  if(!is(cne, "CNE")){
    stop("CNE must be a CNE class object!")
  }
  if(!file.exists(cne@assembly1Fn)){
    stop("The assembly1 twoBit file must exit!")
  }
  if(!file.exists(cne@assembly2Fn)){
    stop("The assembly2 twoBit file must exit!")
  }
  blatOptionsALL <- list("DEF_BLAT_OPT_WSLO"=
                         "-tileSize=9 -minScore=24 -repMatch=16384",
                         "DEF_BLAT_OPT_WSMID"=
                         "-tileSize=10 -minScore=28 -repMatch=4096",
                         "DEF_BLAT_OPT_WSHI"=
                         "-tileSize=11 -minScore=30 -repMatch=1024")
  winSize <- cne@window
  if(is.null(blatOptions)){
    if(winSize > 45L)
      blatOptions <- blatOptionsALL[["DEF_BLAT_OPT_WSHI"]]
    else if(winSize > 35L)
      blatOptions <- blatOptionsALL[["DEF_BLAT_OPT_WSMID"]]
    else
      blatOptions <- blatOptionsALL[["DEF_BLAT_OPT_WSLO"]]
  }
  if(cutIdentity > 100  || cutIdentity < 0)
    stop("cutIdentity must be between 0 and 100!")
  cutoffs1 <- cne@cutoffs1
  cutoffs2 <- cne@cutoffs2
  
  .run_blat <- function(cne, cutIdentity, whichAssembly=c("first", "last"),
                        blatOptions){
    whichAssembly <- match.arg(whichAssembly)
    temp_cne <- tempfile(pattern="cne-")
    temp_psl <- tempfile(pattern="psl-")
    # For Blat, the start is 0-based and end is 1-based. 
    # So make cne's coordinates to comply with it.
    if(whichAssembly == "first"){
      assemblyTwobit <- cne@assembly1Fn
      cneDataFrame <- paste0(assemblyTwobit, ":", 
                    as.character(seqnames(first(cne@CNEMerged))),
                    ":", format(start(first(cne@CNEMerged))-1,
                                trim=TRUE, scientific=FALSE),
                    "-", format(end(first(cne@CNEMerged)),
                                trim=TRUE, scientific=FALSE))
    }else{
      assemblyTwobit <- cne@assembly2Fn
      cneDataFrame <- paste0(assemblyTwobit, ":", 
                    as.character(seqnames(second(cne@CNEMerged))), 
                    ":", format(start(second(cne@CNEMerged))-1,
                                trim=TRUE, scientific=FALSE),
                    "-", format(end(second(cne@CNEMerged)),
                                trim=TRUE, scientific=FALSE))
    }
    cneDataFrame <- unique(cneDataFrame)
    writeLines(cneDataFrame, con=temp_cne)
    cmd <- paste0(cne@aligner, " ", blatOptions, " ",
                  "-minIdentity=", cutIdentity,
                  " ", assemblyTwobit, " ", temp_cne, " ", temp_psl)
    my.system(cmd)
    unlink(temp_cne)
    return(temp_psl)
  }
  psl1Fn <- .run_blat(cne, cutIdentity, "first", blatOptions)
  psl2Fn <- .run_blat(cne, cutIdentity, "last", blatOptions)
  psl1 <- read.table(psl1Fn, header=FALSE, sep="\t", skip=5, as.is=TRUE)
  psl2 <- read.table(psl2Fn, header=FALSE, sep="\t", skip=5, as.is=TRUE)
  colnames(psl1) <- colnames(psl2) <- 
    c("matches", "misMatches", "repMatches", "nCount", 
      "qNumInsert", "qBaseInsert", "tNumInsert", "tBaseInsert", 
      "strand", "qName", "qSize", "qStart", "qEnd", "tName", 
      "tSize", "tStart", "tEnd", "blockCount", "blockSizes", 
      "qStarts", "tStarts")
  ##psl1 = subset(psl1, matches/qSize >= cutIdentity/100)
  ##for some reason, subset will cause error duing the Bioconductor build
  psl1 <- psl1[psl1$matches / psl1$qSize >= cutIdentity/100, ]
  ##psl2 = subset(psl2, matches/qSize >= cutIdentity/100)
  psl2 <- psl2[psl2$matches / psl2$qSize >= cutIdentity/100, ]
  psl1 <- table(psl1[ , "qName"])
  psl2 <- table(psl2[ , "qName"])
  CNEtNameIndex <- unname(psl1[paste0(as.character(
                                        seqnames(first(cne@CNEMerged))),
                                      ":", format(start(first(cne@CNEMerged))-1,
                                                  trim=TRUE, scientific=FALSE),
                                      "-", format(end(first(cne@CNEMerged)),
                                                  trim=TRUE, scientific=FALSE))
                               ]
                          )
  CNEtNameIndex[is.na(CNEtNameIndex)] <- 0
  CNEqNameIndex <- unname(psl2[paste0(as.character(
                                        seqnames(second(cne@CNEMerged))), 
                                      ":", format(start(second(cne@CNEMerged))-1,
                                                  trim=TRUE, scientific=FALSE),
                                      "-", format(end(second(cne@CNEMerged)),
                                                  trim=TRUE, scientific=FALSE))
                               ]
                          )
  CNEqNameIndex[is.na(CNEqNameIndex)] <- 0
  cneFinal <- cne@CNEMerged[as.integer(CNEtNameIndex) <= cutoffs1 &
                            as.integer(CNEqNameIndex) <= cutoffs2]
  BiocGenerics:::replaceSlots(cne, CNEFinal=cneFinal)
}