R/GroupHeatSpot_internal.R

In tangchao7498/ICAS:

#' @title ReadBam
#' @importFrom Rsamtools scanBamFlag
#' @importFrom Rsamtools ScanBamParam
#' @importFrom GenomicRanges reduce
#' @importFrom BiocGenerics start
#' @importFrom BiocGenerics end
#' @importFrom Rsamtools BamFile
#' @importFrom GenomicAlignments readGAlignments
#' @importFrom GenomeInfoDb seqlevels
#' @importFrom GenomicAlignments coverage
#' @importFrom GenomicAlignments strand
#' @importFrom BiocGenerics start end

ReadBam <- function(BamFile, region, ignore.strand = T, singleEnd = F, antisense = T, strandMode = 2, mapqFilter = 0, ...) {
  flag <- Rsamtools::scanBamFlag(isSecondaryAlignment = FALSE, isNotPassingQualityControls = FALSE, isUnmappedQuery = FALSE, ...)
  sbp <- Rsamtools::ScanBamParam(flag = flag, mapqFilter = mapqFilter, which = region)

  if(ignore.strand) {
    if(singleEnd) {
      galp0 <- GenomicAlignments::readGAlignments(Rsamtools::BamFile(BamFile), use.names = FALSE, param = sbp)
    } else {
      galp0 <- GenomicAlignments::readGAlignmentPairs(Rsamtools::BamFile(BamFile), use.names = FALSE, param = sbp)
    }
  } else {
    if(singleEnd) {
      if(antisense) {
        galp0 <- GenomicAlignments::readGAlignments(Rsamtools::BamFile(BamFile), use.names = FALSE, param = sbp)
        GenomicAlignments::strand(galp0) <- ifelse(GenomicAlignments::strand(galp0) == "*", "*", ifelse(GenomicAlignments::strand(galp0) == "+", "-", "+"))
      } else {
        galp0 <- GenomicAlignments::readGAlignments(Rsamtools::BamFile(BamFile), use.names = FALSE, param = sbp)
      }
    } else {
      galp0 <- GenomicAlignments::readGAlignmentPairs(Rsamtools::BamFile(BamFile), use.names = FALSE, param = sbp, strandMode = strandMode)
    }
  }

  galp0 <- galp0[!is.na(seqnames(galp0))]
  GenomeInfoDb::seqlevels(galp0) <- GenomeInfoDb::seqlevels(region)

  if(ignore.strand) {
    galp0 <- GenomicAlignments::coverage(galp0)[[1]]
  } else {
    galp_plus <- GenomicAlignments::coverage(galp0[GenomicAlignments::strand(galp0) == "+"])[[1]]
    galp_minus <- GenomicAlignments::coverage(galp0[GenomicAlignments::strand(galp0) == "-"])[[1]]
    galp0 <- GenomicAlignments::coverage(galp0)[[1]]
  }

  if(ignore.strand) {
    BiocGenerics::strand(region) <- "*"
  }

  if(S4Vectors::runValue(BiocGenerics::strand(region)) == "*") {
    depth <- galp0[BiocGenerics::start(region):BiocGenerics::end(region)]
  } else {
    if(S4Vectors::runValue(BiocGenerics::strand(region)) == "+") {
      depth <- galp_plus[BiocGenerics::start(region):BiocGenerics::end(region)]
    } else {
      depth <- galp_minus[BiocGenerics::start(region):BiocGenerics::end(region)]
    }
  }
  return(depth)
}

#' @title CumPos

CumPos <- function(x) {
  cumsum(c(0, x[-length(x)])) + c(x/2)
}


#' @title PosReadTest
#' @importFrom lmtest waldtest
#'

PosReadTest <- function(depth, group) {
  model1 <- tryCatch(aov(depth ~ group), error = function(e) NULL)
  return(tryCatch(lmtest::waldtest(model1, test = "F")$`Pr(>F)`[2], error = function(e) 0))
}