library (vcfR)

MT12VCF <- function (inputFile, outFile){

  inputDF <- read.delim (inputFile, stringsAsFactors=F)

  if ((length(unique(inputDF$normal_name)) > 1) ||
      (length(unique(inputDF$tumor_name)) > 1)){
    stop ("This function assumes there's only one normal-tumor pair in the input.")
  }

  MuTect1StdHeaders <- '##fileformat=VCFv4.1
##FILTER=<ID=PASS,Description="Accept as a confident somatic mutation">
##FILTER=<ID=REJECT,Description="Rejected as a confident somatic mutation">
##FORMAT=<ID=AD,Number=.,Type=Integer,Description="Allelic depths for the ref and alt alleles in the order listed">
##FORMAT=<ID=BQ,Number=A,Type=Float,Description="Average base quality for reads supporting alleles">
##FORMAT=<ID=DP,Number=1,Type=Integer,Description="Approximate read depth (reads with MQ=255 or with bad mates are filtered)">
##FORMAT=<ID=FA,Number=A,Type=Float,Description="Allele fraction of the alternate allele with regard to reference">
##FORMAT=<ID=GQ,Number=1,Type=Integer,Description="Genotype Quality">
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##FORMAT=<ID=PL,Number=G,Type=Integer,Description="Normalized, Phred-scaled likelihoods for genotypes as defined in the VCF specification">
##FORMAT=<ID=SS,Number=1,Type=Integer,Description="Variant status relative to non-adjacent Normal,0=wildtype,1=germline,2=somatic,3=LOH,4=post-transcriptional modification,5=unknown">
##INFO=<ID=DB,Number=0,Type=Flag,Description="dbSNP Membership">
##INFO=<ID=MQ0,Number=1,Type=Integer,Description="Total Mapping Quality Zero Reads">
##INFO=<ID=SOMATIC,Number=0,Type=Flag,Description="Somatic event">
##INFO=<ID=VT,Number=1,Type=String,Description="Variant type, can be SNP, INS or DEL">
##INFO=<ID=t_lod_fstar,Number=1,Type=Float,Description=" Log of (likelihood tumor event is real / likelihood event is sequencing error )">'

  MuTect1StdHeadersVec <- unlist(strsplit(MuTect1StdHeaders, "\n", fixed=TRUE))

  dbsnpField <- sapply(inputDF$dbsnp_site, 
                       function(x) ifelse(x=="NOVEL", "", "DB;"))
  somaticField <- sapply(inputDF$judgement, 
                         function(x) ifelse(x=="KEEP", "SOMATIC;", ""))
  filterField <- sapply(inputDF$judgement, 
                        function(x) ifelse(x=="KEEP", "PASS", "REJECT"))

  vcfFixDf <- data.frame (
    CHROM = inputDF$contig,
    POS = inputDF$position,
    ID = NA_character_,
    REF = inputDF$ref_allele,
    ALT = inputDF$alt_allele,
    QUAL = inputDF$score,
    FILTER = filterField,
    INFO = paste(sep="",
                 "MQ0=",inputDF$map_q0_reads, ";",
                 dbsnpField, somaticField,
                 "t_lod_fstar=", inputDF$t_lod_fstar,";",
                 "VT=SNP"),
    stringsAsFactors = FALSE
  )

  vcfGTdf <- data.frame (
    FORMAT = "GT:AD:BQ:DP:FA",
    NORMAL = paste (0, 
                    paste (inputDF$n_ref_count,inputDF$n_alt_count, sep=","),
                    ".",
                    inputDF$n_ref_count+inputDF$n_alt_count,
                    inputDF$n_alt_count / 
                      (inputDF$n_ref_count+inputDF$n_alt_count),
                    sep=":"
    ),
    TUMOR = paste ('0/1', 
                   paste (inputDF$t_ref_count,inputDF$t_alt_count, sep=","),
                   paste (inputDF$t_ref_sum / inputDF$t_ref_count,
                          inputDF$t_alt_sum / inputDF$t_alt_count, sep=","),
                   inputDF$t_ref_count+inputDF$t_alt_count,
                   inputDF$t_alt_count / 
                     (inputDF$t_ref_count+inputDF$t_alt_count),
                   sep=":"
    ),
    stringsAsFactors = F
  )

  vcfGTdf[inputDF$judgement=="PASS", "FORMAT"] <- paste (
    vcfGTdf[inputDF$judgement=="PASS", "FORMAT"], ":SS", sep=""
  )

  vcfGTdf[inputDF$judgement=="PASS", "NORMAL"] <- paste (
    vcfGTdf[inputDF$judgement=="PASS", "NORMAL"], ":0", sep=""
  )

  vcfGTdf[inputDF$judgement=="PASS", "TUMOR"] <- paste (
    vcfGTdf[inputDF$judgement=="PASS", "TUMOR"], ":2", sep=""
  )

  rownames(vcfGTdf)[c(2,3)] <- inputDF[1, c("normal_name", "tumor_name")]

  vcfRObj <- new(Class = "vcfR")
  vcfRObj@meta <- MuTect1StdHeadersVec
  vcfRObj@fix <- as.matrix (vcfFixDf)
  vcfRObj@gt <- as.matrix(vcfGTdf)

  write.vcf (vcfRObj, outFile)

}