### R code from vignette source 'vignettes/GWASTools/inst/doc/VCF.Rnw'
### Encoding: UTF-8

###################################################
### code chunk number 1: VCF.Rnw:40-44
###################################################
library(GWASTools)
vcffile <- system.file("extdata", "sequence.vcf", package="SNPRelate")
gdsfile <- "snps.gds"
convertVcfGds(vcffile, gdsfile)


###################################################
### code chunk number 2: VCF.Rnw:50-62
###################################################
(gds <- GdsGenotypeReader(gdsfile))

getScanID(gds)
getSnpID(gds)
getChromosome(gds)
getPosition(gds)
getVariable(gds, "sample.name")
getVariable(gds, "snp.rs.id")
getVariable(gds, "snp.allele")
getGenotype(gds)

close(gds)


###################################################
### code chunk number 3: VCF.Rnw:65-66
###################################################
unlink(gdsfile)


###################################################
### code chunk number 4: VCF.Rnw:77-80
###################################################
library(VariantAnnotation)
vcffile <- system.file("extdata", "ex2.vcf", package="VariantAnnotation")
(vcf <- readVcf(vcffile, "hg18"))


###################################################
### code chunk number 5: VCF.Rnw:87-97
###################################################
# width of ref seq
rw <- width(ref(vcf))
# width of first alt seq
aw <- unlist(lapply(alt(vcf), function(x) {width(x[1])}))
# number of alternate genotypes
nalt <- elementLengths(alt(vcf))
# select only bi-allelic SNPs (monomorphic OK, so aw can be 0 or 1)
snp <- rw == 1 & aw <= 1 & nalt == 1
# subset vcf
vcf <- vcf[snp,]


###################################################
### code chunk number 6: VCF.Rnw:101-107
###################################################
chrom <- as.vector(seqnames(rowData(vcf)))
position <- start(ranges(rowData(vcf)))
alleleA <- as.character(ref(vcf))
alleleB <- as.character(unlist(alt(vcf)))
QUAL <- qual(vcf)
FILTER <- filt(vcf)


###################################################
### code chunk number 7: VCF.Rnw:111-117
###################################################
chromosome <- chrom
chromosome[chrom == "X"] <- 23
chromosome[chrom == "XY"] <- 24 # pseudoautosomal
chromosome[chrom == "Y"] <- 25
chromosome[chrom == "M"] <- 26 # mitochrondrial
chromosome <- as.integer(chromosome)


###################################################
### code chunk number 8: VCF.Rnw:126-143
###################################################
# make annotation
snp.df <- data.frame(chromosome, position, alleleA, alleleB, QUAL, FILTER,
                     row.names=rownames(vcf), stringsAsFactors=FALSE)
snp.ord <- order(snp.df$chromosome, snp.df$position)
snp.df <- snp.df[snp.ord,]
snp.df$snpID <- 1:nrow(snp.df)
cols <- c("snpID", "chromosome", "position", "alleleA", "alleleB", "QUAL", "FILTER")
hdr <- exptData(vcf)$header
desc <- c("unique integer ID = snp dimension values in NetCDF file",
  "chromosome as integer where X=23, XY=24, Y=25, M=26",
  "base position",
  "allele A (the reference allele)",
  "allele B (the alternate allele)",
  "phred-scaled quality score for the assertion made in ALT. i.e. -10log_10 prob(call in ALT is wrong)",
  paste(paste(row.names(fixed(hdr)$FILTER), fixed(hdr)$FILTER$Description, sep="="), collapse="; "))
meta <- data.frame("labelDescription"=desc, row.names=cols)
snp.df <- snp.df[,cols]


###################################################
### code chunk number 9: VCF.Rnw:151-171
###################################################
info <- info(vcf)
info.meta <- info(hdr)
vals <- list()
for (i in 1:nrow(info.meta)) {
  field <- row.names(info.meta)[i]
  number <- info.meta$Number[i]
  if (number == "A") {
    # one element per alternate allele
    # should only have one value per variant
    dat <- unlist(info[[field]])
    if (length(dat) > nrow(info)) next
    vals[[field]] <- dat
  } else if (number %in% 0:1) {
    vals[[field]] <- info[[field]]
  }
}
snp.df <- cbind(snp.df, vals)
meta.vals <- as.data.frame(info.meta[names(vals), "Description", drop=FALSE])
names(meta.vals) <- "labelDescription"
meta <- rbind(meta, meta.vals)


###################################################
### code chunk number 10: VCF.Rnw:175-177
###################################################
library(GWASTools)
snpAnnot <- SnpAnnotationDataFrame(snp.df, meta)


###################################################
### code chunk number 11: VCF.Rnw:182-186
###################################################
scan.df <- data.frame("scanID"=1:ncol(vcf), row.names=colnames(vcf))
meta.scan <- data.frame("labelDescription"="unique integer ID = sampleID in NetCDF file",
                        row.names=names(scan.df))
scanAnnot <- ScanAnnotationDataFrame(scan.df, meta.scan)


###################################################
### code chunk number 12: VCF.Rnw:192-202
###################################################
geno <- geno(vcf)$GT
geno.ab <- geno
geno.ab[,] <- NA
geno.ab[geno %in% c("0/0", "0|0")] <- 2 # AA = REF/REF
geno.ab[geno %in% c("0/1", "0|1", "1/0", "1|0")] <- 1 # AB = REF/ALT
geno.ab[geno %in% c("1/1", "1|1")] <- 0 # BB = ALT/ALT
# new snp order (by chrom and position)
geno.ab <- geno.ab[snp.ord,]
# geno.ab should be integer, not character
mode(geno.ab) <- "integer"


###################################################
### code chunk number 13: VCF.Rnw:206-207
###################################################
gq <- geno(vcf)$GQ


###################################################
### code chunk number 14: VCF.Rnw:212-220
###################################################
ncfile <- "ex2.nc"
ncdfCreate(snp.annotation=snp.df, ncdf.filename=ncfile,
           n.samples=ncol(vcf), variables=c("genotype", "quality"))
nc <- open.ncdf(ncfile, write=TRUE)
put.var.ncdf(nc, "sampleID", scanAnnot$scanID, start=1, count=-1)
put.var.ncdf(nc, "genotype", geno.ab)
put.var.ncdf(nc, "quality", gq)
close.ncdf(nc)


###################################################
### code chunk number 15: VCF.Rnw:225-232
###################################################
nc <- NcdfGenotypeReader(ncfile)
(genoData <- GenotypeData(nc, snpAnnot=snpAnnot, scanAnnot=scanAnnot))
pData(snpAnnot)
pData(scanAnnot)
getGenotype(genoData)
getVariable(genoData, "quality")
close(genoData)


###################################################
### code chunk number 16: VCF.Rnw:235-236
###################################################
unlink(ncfile)