data-raw/old-code/ReadVCFAndBAMsAndProcess_old.R

In steverozen/DBSverify: Analyze Next Generation Sequencing Data to Verify Somatic Double-BaseSubstituion (DBS) Calls

#' Determine whether sequencing reads in fact support DBSs generated by processing a VCF file.
#'
#' @param vcf.name The name of the VCF file to analyze.
#'
#' @param Nbam.name The name of the BAM file for the normal sample corresponding to \code{vcf.name}.
#'
#' @param Tbam.name The name of the BAM file for the tumor sample corresponding to \code{vcf.name}.
#'
#' @param variant.caller One of \code{"strelka"}, \code{"mutect"}, \code{"PCAWG"}, or \code{"unknown"}.
#'  If \code{variant.caller} is \code{"unknown"}, then the variant caller must have called
#'  DBSs as such. For \code{"strelka"} and {"mutect"} DBS analysis is as in \code{ICAMS};
#'  Strelka does not call DBSs, and DBSs are inferred from adjacent single base substitutions.
#'  This is hte case for \code{"PCAWG"} as well.
#'
#' @param outfile If not \code{NULL} then write the "evaluated" VCF to \code{outfile};
#'  otherwise write it to \code{paste0(vcf.name, "_evaluated.vcf")}.
#'
#' @param num.cores Number of cores to use while reading the VCF. Set this to 1 for testing
#'  and on MS Windows.
#'
#' @details Creates a new VCF file.
#'  The additional information in this VCF files is described in \code{\link{VerifyDBSVcf}}.
#'
#' @return Invisibly, a list, the first element of which is the name of the file created, the
#'   second of which is an in-memory representation of the VCF as a \code{data.table}.
#'
#' @export
#'

ReadVCFAndBAMsAndProcess_old <- function(vcf.name, Nbam.name,
                                     Tbam.name,
                                     variant.caller, num.cores = 10, outfile = NULL) {

  CheckBAM(Nbam.name)
  CheckBAM(Tbam.name)
  if (!file.exists(vcf.name)) stop("VCF file ", vcf.name, "does not exist")

  if (variant.caller %in% c("strelka", "mutect", "unknown")) {
    get.vaf.function <- NULL
  } else if (variant.caller == "PCAWG") {
    variant.caller <- "unknown"
    get.vaf.function <- GetPCAWGVAF
  }

  vcf.list <-ICAMS::ReadAndSplitVCFs(vcf.name,
                                     variant.caller = variant.caller,
                                     num.of.cores   = num.cores,
                                     get.vaf.function = get.vaf.function,
                                     max.vaf.diff   = 1)

  DBS.vcf <- vcf.list$DBS[[1]]

  evaluated.vcf <- VerifyDBSVcf(DBS.vcf, Nbam.name = Nbam.name, Tbam.name = Tbam.name)

  if (is.null(outfile)) {
    outfile <- paste0(vcf.name, "_evaluated.vcf")
  }

  cat("#", file = outfile)
  suppressWarnings(
    utils::write.table(
      evaluated.vcf,
      file = outfile,
      sep="\t",
      quote = FALSE,
      row.names = FALSE,
      append = TRUE))


  invisible(list(vcf.name = outfile, evaluated.vcf = evaluated.vcf))
}