R/vcfToMAF.R

In likelet/CaMutQC: An R Package for Comprehensive Filtration and Selection of Cancer Somatic Mutations

#' vcfToMAF
#' @description Format transformation from VCF to MAF.
#'
#' @param vcfFile Directory of a VCF file, or the path to several VCF files
#' that is going to be transformed. Files should be in .vcf or .vcf.gz format.
#' @param multiVCF Logical, whether the input is a path that leads to several
#' VCFs that come from multi-region/sample/caller sequencing. Default: FALSE
#' @param inputStrelka The type of variants ('INDEL' or 'SNV') in VCF file 
#' if it is from Strelka. Default: FALSE
#' @param writeFile Whether to directly write MAF file to the disk. If FALSE,
#' a MAF data frame will be returned. If TRUE, a MAF file will be saved.
#' Default: FALSE.
#' @param MAFfile File name of the exported MAF file, if writeFile is set as
#' TRUE.
#' @param MAFdir Directory of the exported MAF file, if writeFile is set as
#' TRUE.
#' @param tumorSampleName Name of the tumor sample(s) in the VCF file(s).
#' If it is set as 'Extracted', tumorSampleName would be extracted 
#' automatically from the VCF file. Default: 'Extracted'.
#' @param normalSampleName Name the normal sample in the VCF file.
#' If it is set as 'Extracted', normalSampleName would be extracted 
#' automatically from the VCF file. Default: 'Extracted'.
#' @param ncbiBuild The reference genome used for the alignment, which will be
#' presented as value in 'NCBIbuild' column in MAF file. Default: 'GRCh38'.
#' @param MAFcenter One or more genome sequencing center reporting the variant,
#' which will be presented as value in 'Center' column in MAF. Default: '.'.
#' @param MAFstrand Genomic strand of the reported allele, which will be
#' presented as value in 'Strand' column in MAF file. Default: '+'.
#' @param filterGene Logical. Whether to filter variants without Hugo Symbol.
#' Default: FALSE
#' @param simplified Logical. Whether to extract the first thirteen columns 
#' after converting to MAF file. Default: FALSE
#'
#' @import vcfR org.Hs.eg.db clusterProfiler stringr dplyr utils
#' @importFrom stats na.omit
#' @return A detailed MAF data frame
#' @export vcfToMAF
#'
#' @examples
#' maf <- vcfToMAF(system.file("extdata", "WES_EA_T_1_mutect2.vep.vcf",
#' package="CaMutQC"))


vcfToMAF <- function(vcfFile, multiVCF = FALSE, inputStrelka = FALSE,
                     writeFile = FALSE, MAFfile = 'MAF.maf', MAFdir = './',
                     tumorSampleName = 'Extracted',
                     normalSampleName = 'Extracted', ncbiBuild = 'Extracted',
                     MAFcenter = '.', MAFstrand = '+', filterGene = FALSE,
                     simplified = FALSE){
    # if inputs are multi-sample/multi-caller VCFs
    if (multiVCF){
        filenames <- list.files(vcfFile, pattern="*.vep.vcf", full.names=TRUE)
        if (length(filenames) == 0){
          stop('No VCF file detected under the path you offered.')
        }else{
          # transform each vcf to maf separately
          dfs <- lapply(filenames, vcfhelper, tumorSampleName=tumorSampleName,
                        normalSampleName=normalSampleName, ncbiBuild=ncbiBuild,
                        MAFcenter=MAFcenter, MAFstrand=MAFstrand)
          maf <- do.call("rbind", dfs)
        }
    }else{
      # read vcf File
      ## check whether the file exists
      if (!(file_test('-f', vcfFile))){
          stop('VCF file doesn\'t exist, please check your input.')
      }
      ## check the file format
      nameSplit <- strsplit(vcfFile, split=".", fixed=TRUE)[[1]]
      suffixLast <- nameSplit[length(nameSplit)]
      suffixSecond <- nameSplit[length(nameSplit) - 1]
      if (suffixLast != 'vcf') {
          stop('Please input file in .vcf or .gz.vcf format')
      }
      maf <- vcfhelper(vcfFile, tumorSampleName=tumorSampleName,
                       normalSampleName=normalSampleName, ncbiBuild=ncbiBuild,
                       MAFcenter=MAFcenter, MAFstrand=MAFstrand,
                       inputStrelka=inputStrelka)
    }
    # simplify MAF file
    if (simplified) {
        maf <- maf[, seq_len(13)]
    }
    if (filterGene){maf <- maf[which(maf$Hugo_Symbol != ''), ]}
    if (writeFile) {
        message('The generated MAF file has been saved.')
        write.table(maf, paste0(MAFdir, MAFfile), sep="\t",
                    quote=FALSE, row.names=FALSE)
        return(maf)
    }else{return(maf)}
}

## helper function for the transformation from VCF to MAF
vcfhelper <- function(vcfFile, tumorSampleName = 'Extracted',
                      normalSampleName = 'Extracted', ncbiBuild = 'Extracted',
                      MAFcenter = '.', MAFstrand = '+', inputStrelka = FALSE) {
    invisible(capture.output(annoVCF <- read.vcfR(vcfFile)))
    mes <- paste0(basename(vcfFile), ' loaded successfully!')
    message(mes)
    ## convert to data frames and store useful information
    vcfHeader <- as.data.frame(annoVCF@meta)
    vcfMain <- as.data.frame(annoVCF@fix)
    vcfAdditional <- as.data.frame(annoVCF@gt)
    # build data frame in MAF format
    maf <- as.data.frame(matrix(ncol=46, nrow=nrow(vcfMain)))
    colnames(maf)<-c('Hugo_Symbol', 'Entrez_Gene_Id', 'Center', 'NCBI_Build',
                    'Chromosome', 'Start_Position', 'End_Position', 'Strand',
                    'Variant_Classification', 'Variant_Type',
                    'Reference_Allele','Tumor_Seq_Allele1','Tumor_Seq_Allele2',
                    'dbSNP_RS', 'dbSNP_Val_Status', 'Tumor_Sample_Barcode',
                    'Matched_Norm_Sample_Barcode', 'Match_Norm_Seq_Allele1',
                    'Match_Norm_Seq_Allele2', 'Tumor_Validation_Allele1',
                    'Tumor_Validation_Allele2', 'Match_Norm_Validation_Allele1',
                    'Match_Norm_Validation_Allele2', 'Verification_Status',
                    'Validation_Status', 'Mutation_Status', 'Sequencing_Phase',
                    'Sequence_Source', 'Validation_Method', 'Score',
                    'BAM_File', 'Sequencer', 'Tumor_Sample_UUID',
                    'Matched_Norm_Sample_UUID', 'HGVSc', 'HGVSp', 'HGVSp_Short',
                    'Transcript_ID', 'Exon_Number', 't_depth', 't_ref_count',
                    't_alt_count', 'n_depth', 'n_ref_count', 'n_alt_count',
                    'all_effects'
    )
    # handle conditions when there is no variants in the VCF file
    if (nrow(vcfMain) == 0) {
        stop("No variants in VCF file. Program stops here.")
    }else{
        # check whether there is a paired normal sample in the VCF
        if (ncol(vcfAdditional) < 3) {
            stop("Missing columns in VCF, a paired normal sample is required.")
        }else{
            message('VCF to MAF conversion is in process...')
        }
    }
    # assign chromosome, Start_Position,Reference_Allele, Tumor_Seq_Allele1 info
    maf[, 5] <- vcfMain$CHROM
    maf[, 6] <- vcfMain$POS
    maf[, 11] <- vcfMain$REF
    maf[, 13] <- vcfMain$ALT
    # get INFO information
    infoAll <- readINFO(vcfHeader, vcfMain)
    infoFrame <- as.data.frame(infoAll[[1]])
    infoInfo <- as.data.frame(infoAll[[2]])
    # process CSQ column
    ## extract CSQ header and set as colnames
    csqHeaders <- strsplit(infoFrame[which(infoFrame$ID == 'CSQ'), 4],
                            split="Format: ")[[1]][2]
    csqHeaders <- strsplit(csqHeaders, split="\\|")[[1]]
    csqInfo <- as.data.frame(matrix(ncol=length(csqHeaders), nrow=nrow(maf)))
    colnames(csqInfo) <- csqHeaders
    ## use built functions to select proper transcripts
    csqInfo <- selectTrans(csqInfo, infoInfo)
    # fill in other information in MAF file
    ## ID conversion
    IDs <- withWarningsAndMessages(bitr(csqInfo$Gene, fromType="ENSEMBL", 
                       toType="ENTREZID", OrgDb=org.Hs.eg.db))$value
    ## center
    maf[, 3] <- MAFcenter
    ## NCBI build
    if (ncbiBuild == 'Extracted') {
      maf[, 4] <- na.omit(str_extract(vcfHeader$`annoVCF@meta`,
                                      pattern="(GRCh)[:digit:]+"))[1]
    } else {
      maf[, 4] <- ncbiBuild
    }
    ## tumor sample name and normal sample name
    if (tumorSampleName == 'Extracted'){
      # mutect format: tumor_sample=
      tumorSampleLine <- vcfHeader$`annoVCF@meta`[grep('tumor_sample=',
                                                    vcfHeader$`annoVCF@meta`)]
      if (length(tumorSampleLine) != 0){
        tumorSampleName <- strsplit(tumorSampleLine, split='=')[[1]][2]
      }else{
        # MuSE format: TUMOR=
        tumorSampleLine <- vcfHeader$`annoVCF@meta`[grep('TUMOR=',
                                                    vcfHeader$`annoVCF@meta`)]
        if (length(tumorSampleLine) != 0){
          tumorSampleName <- strsplit(strsplit(tumorSampleLine, 
                                      split=',')[[1]][1], split="=")[[1]][3]
        }else{
          tumorSampleName <- 'TUMOR'
        }
      }
    }
    if (normalSampleName == 'Extracted'){
      # mutect format: normal_sample=
      normalSampleLine <- vcfHeader$`annoVCF@meta`[grep('normal_sample=',
                                                  vcfHeader$`annoVCF@meta`)]
      if (length(normalSampleLine) != 0){
        normalSampleName <- strsplit(normalSampleLine, split='=')[[1]][2]
      }else{
        # MuSE format: NORMAL=
        normalSampleLine <- vcfHeader$`annoVCF@meta`[grep('NORMAL=',
                                        vcfHeader$`annoVCF@meta`)]
        if (length(normalSampleLine) != 0){
          normalSampleName <- strsplit(strsplit(normalSampleLine, 
                                     split=',')[[1]][1], split="=")[[1]][3]
        }else{
          normalSampleName <- 'NORMAL'
        }
      }
    }
    colnames(vcfAdditional)[colnames(vcfAdditional) == 'TUMOR'] <- tumorSampleName
    colnames(vcfAdditional)[colnames(vcfAdditional) == 'NORMAL'] <- normalSampleName
    ## strand
    maf[, 8] <- MAFstrand
    ## Hugo_Symbol
    maf[, 1] <- csqInfo[, 4]
    ## add VAF column
    maf <- cbind(maf, VAF = 0)
    maf[, c(15, 18:34, 37, 46)] <- '.'
    for (i in seq_len(nrow(maf))) {
      ## ENTREZID
      maf[i, 2] <- IDs$ENTREZID[which(IDs$ENSEMBL == csqInfo$Gene[i])][1]
      ## get correct variant Position, Variant_Type, Ref allele and Alt allele
      posType <- getVarFeature(maf[i, 6], maf[i, 11], maf[i, 13], csqInfo[i, 1])
      maf[i, 6] <- posType[[1]] # start
      maf[i, 7] <- posType[[2]] # end
      maf[i, 10] <- posType[[3]] # Variant_Type
      maf[i, 11] <- posType[[4]] # correct ref allele
      maf[i, 13] <- posType[[5]] # correct alt allele
      inframe <- posType[[6]] # inframe indicator
      ## Variant_Class
      ### select consequence in CSQ first
      cons <- strsplit(csqInfo[i, 2], split='&')[[1]]
      csqInfo[i, 2] <- cons[which.min(unlist(vapply(cons,
                                              getConsequencePriority, 10)))]
      maf[i, 9] <- getVarClass(csqInfo[i, 2], maf[i, 10], inframe)
      ## get AD and DP in INFO or FORMAT
      ## if RD field is contained in the INFO column
      if (any(strsplit(vcfAdditional[i, 1], ":")[[1]] == 'RD')) {
        RDLoc <- strsplit(vcfAdditional[i, 1], ":")[[1]] == 'RD'
        ADLoc <- strsplit(vcfAdditional[i, 1], ":")[[1]] == 'AD'
        DPLoc <- strsplit(vcfAdditional[i, 1], ":")[[1]] == 'DP'
        maf[i, "t_ref_count"] <- strsplit(strsplit(vcfAdditional[i, 
                                  tumorSampleName], ":")[[1]][RDLoc], ",")[[1]]
        maf[i, "t_alt_count"] <- strsplit(strsplit(vcfAdditional[i, 
                                  tumorSampleName], ":")[[1]][ADLoc], ",")[[1]]
        maf[i, "n_ref_count"] <- strsplit(strsplit(vcfAdditional[i, 
                                normalSampleName], ":")[[1]][RDLoc], ",")[[1]]
        maf[i, "n_alt_count"] <- strsplit(strsplit(vcfAdditional[i, 
                                normalSampleName], ":")[[1]][ADLoc], ",")[[1]]
        maf[i, "t_depth"] <- strsplit(strsplit(vcfAdditional[i, 
                                  tumorSampleName], ":")[[1]][DPLoc], ",")[[1]]
        maf[i, "n_depth"] <- strsplit(strsplit(vcfAdditional[i, 
                                normalSampleName], ":")[[1]][DPLoc], ",")[[1]]
  
        maf[i, 'VAF'] <- as.numeric(maf[i, "t_alt_count"])/
          as.numeric(maf[i, "t_depth"])
      }else{
        ADLoc <- strsplit(vcfAdditional[i, 1], ":")[[1]] == 'AD'
        ## if AD field contains information from both ref and alt
        if (any(ADLoc)){
          AD <- strsplit(strsplit(vcfAdditional[i, tumorSampleName],
                                  ":")[[1]][ADLoc], ",")[[1]][2]
          if (length(grep('DP', strsplit(vcfAdditional[i, 1], ":")[[1]]))){
            DPLoc <- strsplit(vcfAdditional[i, 1], ":")[[1]] == 'DP'
            DP <- as.numeric(strsplit(vcfAdditional[i, tumorSampleName],
                                      ":")[[1]][DPLoc])
            tDP <- as.numeric(strsplit(vcfAdditional[i, tumorSampleName],
                                       ":")[[1]][DPLoc])
            nDP <- as.numeric(strsplit(vcfAdditional[i, normalSampleName],
                                       ":")[[1]][DPLoc])
          }else if('DP' %in% colnames(infoInfo)){
            DP <- as.numeric(infoInfo[i, 'DP'])
            tDP <- DP
            nDP <- sum(as.numeric(strsplit(strsplit(vcfAdditional[i, 
                       normalSampleName],":")[[1]][ADLoc], ",")[[1]]))
          }
          ## t_depth, n_depth, t_ref_count, t_alt_count, n_ref_count, n_alt_count
          tRefAD<-as.numeric(strsplit(strsplit(vcfAdditional[i,tumorSampleName]
                                      , ":")[[1]][ADLoc], ",")[[1]][1])
          nRefAD<-as.numeric(strsplit(strsplit(vcfAdditional[i,normalSampleName]
                                      , ":")[[1]][ADLoc], ",")[[1]][1])
          nAltAD<-as.numeric(strsplit(strsplit(vcfAdditional[i,normalSampleName]
                                      , ":")[[1]][ADLoc], ",")[[1]][2])
          maf[i, "t_ref_count"] <- tRefAD
          maf[i, "t_alt_count"] <- as.numeric(AD)
          maf[i, "n_ref_count"] <- nRefAD
          maf[i, "n_alt_count"] <- nAltAD
          if (!(exists('tDP'))){
            maf[i, "t_depth"] <- tRefAD + as.numeric(AD)
            maf[i, "n_depth"] <- nRefAD + nAltAD
            maf[i, 'VAF'] <- as.numeric(AD)/maf[i, "t_depth"]
          }else{
            maf[i, "t_depth"] <- tDP
            maf[i, "n_depth"] <- nDP
            maf[i, 'VAF'] <- as.numeric(AD)/DP
          }
        }else{
          # if no AD field detected in the VCF, it should be output from Strelka
          if (inputStrelka == 'INDEL'){
            TIRLoc <- strsplit(vcfAdditional[i, 1], ":")[[1]] == 'TIR'
            TARLoc <- strsplit(vcfAdditional[i, 1], ":")[[1]] == 'TAR'
            tTIR <-strsplit(vcfAdditional[i,tumorSampleName], ":")[[1]][TIRLoc]
            tTAR <-strsplit(vcfAdditional[i,tumorSampleName], ":")[[1]][TARLoc]
            nTIR<-strsplit(vcfAdditional[i,normalSampleName], ":")[[1]][TIRLoc]
            nTAR<-strsplit(vcfAdditional[i,normalSampleName], ":")[[1]][TARLoc]
            # use tire1 information
            tAltAD <- as.numeric(strsplit(tTIR, ",")[[1]][1])
            tRefAD <- as.numeric(strsplit(tTAR, ",")[[1]][1])
            nAltAD <- as.numeric(strsplit(nTIR, ",")[[1]][1])
            nRefAD <- as.numeric(strsplit(nTAR, ",")[[1]][1])
          }else if (inputStrelka == 'SNV'){
            refAllele <- paste0(vcfMain[i, 'REF'], "U")
            altAllele <- paste0(vcfMain[i, 'ALT'], "U")
            refLoc <- strsplit(vcfAdditional[i, 1], ":")[[1]] == refAllele
            altLoc <- strsplit(vcfAdditional[i, 1], ":")[[1]] == altAllele
            tRef<-strsplit(vcfAdditional[i, tumorSampleName], ":")[[1]][refLoc]
            tAlt<-strsplit(vcfAdditional[i, tumorSampleName], ":")[[1]][altLoc]
            nRef<-strsplit(vcfAdditional[i, normalSampleName], ":")[[1]][refLoc]
            nAlt<-strsplit(vcfAdditional[i, normalSampleName], ":")[[1]][altLoc]
            tAltAD <- as.numeric(strsplit(tAlt, ",")[[1]][1])
            tRefAD <- as.numeric(strsplit(tRef, ",")[[1]][1])
            nAltAD <- as.numeric(strsplit(nAlt, ",")[[1]][1])
            nRefAD <- as.numeric(strsplit(nRef, ",")[[1]][1])
          }else{
            mes <- paste0("If your input comes from Strelka,",
                          " set 'indel' or 'SNV' for 'inputStrelka' parameter")
            stop(mes)
          }
          maf[i, "t_ref_count"] <- tRefAD
          maf[i, "t_alt_count"] <- tAltAD
          maf[i, "n_ref_count"] <- nRefAD
          maf[i, "n_alt_count"] <- nAltAD
          maf[i, "t_depth"] <- tRefAD + tAltAD
          maf[i, "n_depth"] <- nRefAD + nAltAD
          maf[i, 'VAF'] <- tAltAD/maf[i, "t_depth"]
        }
      }
      ## fill in dbSNP_RS
      if(nchar(csqInfo[i, 'Existing_variation']) == 0) {
        maf[i, 'dbSNP_RS'] <- 'novel'
      } else if (str_detect(csqInfo[i, 'Existing_variation'], 'rs')) {
        extVar <- strsplit(csqInfo[i, 'Existing_variation'], "&")[[1]]
        dbSNPLoc <- str_detect(extVar, 'rs')
        maf[i, 'dbSNP_RS'] <- paste(extVar[dbSNPLoc], collapse=",")
      } else {
        maf[i, 'dbSNP_RS'] <- '.'
      }
      ## HGVSc
      if (nchar(csqInfo[i, 'HGVSc']) != 0) {
          maf[i, 'HGVSc'] <- strsplit(csqInfo[i, 'HGVSc'], split=":")[[1]][2]
      }
      ## HGVSp
      if (nchar(csqInfo[i, 'HGVSp']) != 0) {
          maf[i, 'HGVSp'] <- strsplit(csqInfo[i, 'HGVSp'], split=":")[[1]][2]
      }
    }
    ## Transcript_ID
    maf[, 'Transcript_ID'] <- csqInfo[, 'Feature']
    ## Exon_Number
    maf[, 'Exon_Number'] <- csqInfo[, 'EXON']
    ## set Tumor_Seq_Allele1 same as ref
    maf[, 12] <- maf[, 11]
    maf[, 'Tumor_Sample_Barcode'] <- tumorSampleName
    maf[, 'Matched_Norm_Sample_Barcode'] <- normalSampleName
    maf$t_alt_count <- as.numeric(maf$t_alt_count)
    maf$t_ref_count <- as.numeric(maf$t_ref_count)
    maf$t_depth <- as.numeric(maf$t_depth)
    maf$n_alt_count <- as.numeric(maf$n_alt_count)
    maf$n_ref_count <- as.numeric(maf$n_ref_count)
    maf$n_depth <- as.numeric(maf$n_depth)
    maf1 <- jointMAF(maf[ , seq_len(46)], csqInfo, vcfMain)
    # remane NORMAL and TUMOR column so vcf from different callers can merge
    colnames(vcfAdditional)[which(colnames(vcfAdditional) 
                                   == tumorSampleName)] <- 'tumorSampleInfo'
    colnames(vcfAdditional)[which(colnames(vcfAdditional) 
                                   == normalSampleName)] <- 'normalSampleInfo'
    # handle cases when VAF is NA because t_depth = 0
    maf[is.na(maf$VAF),'VAF'] <- 0
    maf <- cbind(maf1, VAF = maf[ , 'VAF'], vcfAdditional)
    # change column type
    maf$Start_Position <- as.numeric(maf$Start_Position)
    maf$End_Position <- as.numeric(maf$End_Position)
    # change column type without displaying warnings
    maf <- maf %>% mutate_at(vars(78:86, 101:117), as.numeric)
    maf <- cbind(maf, CaTag='0')
    rownames(maf) <- seq_len(nrow(maf))
    message('VCF to MAF conversion has been done successfully!')
    return(maf)
}

# function to avoid warnings and messages generated by clusterProfiler
withWarningsAndMessages <- function(expr) {
    myWarnings <- NULL
    myMessages <- NULL
    wHandler <- function(w) {
        myWarnings <- c(myWarnings, list(w))
        invokeRestart("muffleWarning")
    }
    mHandler <- function(m) {
        myMessages <- c(myMessages, list(m))
        invokeRestart("muffleMessage")
    }
    # process the command additionally
    val <- withCallingHandlers(expr, warning=wHandler, message=mHandler)
    list(value=val, warnings=myWarnings, messages=myMessages)
}