第一步，利用GATK中 VariantsToTable 將vcf文件格式轉(zhuǎn)化

gatk VariantsToTable \
-V test.vcf \   #輸入文件
-R Genome.fa \
-F CHROM -F POS -F REF -F ALT -GF GT \
-O output.table     #輸出文件

第二步：使用如下R腳本將變量表轉(zhuǎn)換為hapmap。

注意：僅僅在單倍體，雙等位基因型SNP中得到測試

Rscript convert-vcf-to-hapmap.R -i output.table -o vcf_converted_to_hapmap.txt

R腳本如下

##service manual
suppressPackageStartupMessages(require(optparse))

option_list = list(
  make_option(c("-i", "--input"), type="character", default=NULL,help="Input file name. The file output from GATK'S VariantsToTable. [required]", 
              metavar="character"),
  make_option(c("-o", "--output"), type="character", default="vcf_converted_to_hapmap.txt",help="output file name [default= %default]", 
              metavar="character")
)
opt_parser = OptionParser(usage = "usage: %prog [options] args",option_list=option_list)
opt = parse_args(opt_parser)

if (is.null(opt$input)){
  print_help(opt_parser)
  #stop("At least one argument must be supplied (input file)", call.=TRUE)
  stop("The input file must be supplied")
}

vcf <- read.table(opt$input, header=TRUE, stringsAsFactors = FALSE)

rs <- paste(vcf$CHROM, vcf$POS, sep = "_")
alleles <- paste(vcf$REF, vcf$ALT, sep = "/")
chrom <- vcf$CHROM
pos <- vcf$POS
strand <- NA
assembly <- NA
center <- NA
protLSID <- NA
assayLSID <- NA
panelLSID <- NA
QCcode <- NA
genotype_df <- vcf[c(-1:-4)]


f1 <- function(x) gsub('/','',x)
f2 <- function(x) gsub('\\.\\.','NN',x)

temp <- data.frame(t(apply(vcf[c(-1:-4)],1,f1)))

genotype_df <- data.frame(t(apply(temp,1,f2)))

hapmap <- cbind(rs, alleles, chrom, pos, strand, assembly, center, protLSID, assayLSID, panelLSID, QCcode, genotype_df)


write.table(hapmap, opt$output, row.names=FALSE, quote=FALSE, sep="\t")

引自wenkaiyan-kevin/vcf-to-hapmap (github.com)

第二個腳本

java -Xmx8g -jar GenomeAnalysisTK.jar \
        -R REFERENCE_GENOME.fa \
        -T VariantsToTable \
         -V INPUT_VCF.vcf \
         -F CHROM -F POS -F ID -F QUAL -F REF -F ALT -F QUAL -F FILTER -F INFO -F FORMAT \
         -o OUTPUT.txt \
         -AMD \
         -GF GT 

    Rscript --vanilla vcf2hapmap.R OUTPUT.vcf HAPMAP.txt

#! /Library/Frameworks/R.framework/
args <- commandArgs(trailingOnly=TRUE)

if (length(args)==0) {
  stop("At least one argument must be supplied (input file).n", call.=FALSE)
} else if (length(args)==1) {
  # default output file
  args[2] <- "vcf_converted_to_hapmap.txt"
}

vcf <- read.table(args[1], header=TRUE, stringsAsFactors = FALSE)

rs <- NA
alleles <- paste(vcf$REF, vcf$ALT, sep = "/")
chrom <- vcf$CHROM
pos <- vcf$POS
strand <- "+"
assembly <- NA
center <- NA
protLSID <- NA
assayLSID <- NA
panelLSID <- NA
Qcode <- NA
genotype_df <- vcf[c(-1:-10)]

hapmap <- cbind(rs, alleles, chrom, pos, strand, assembly, center, protLSID, assayLSID, panelLSID, Qcode, genotype_df)
colnames(hapmap)[[1]]<-"rs#"
colnames(hapmap)[[6]]<-"assembly#"

write.table(hapmap, file=args[2], row.names=FALSE, quote=FALSE, sep="\t")

引自koleaz/convert-vcf-to-hapmap (github.com)