VariantAnnotation.R
In VariantAnnotation: Annotation of Genetic Variants

### R code from vignette source 'VariantAnnotation.Rnw'

###################################################
### code chunk number 1: style
###################################################
BiocStyle::latex()


###################################################
### code chunk number 2: options
###################################################
options(width=72)


###################################################
### code chunk number 3: readVcf
###################################################
library(VariantAnnotation)
fl <- system.file("extdata", "chr22.vcf.gz", package="VariantAnnotation")
vcf <- readVcf(fl, "hg19")
vcf


###################################################
### code chunk number 4: readVcf_showheader
###################################################
header(vcf)


###################################################
### code chunk number 5: headeraccessors
###################################################
samples(header(vcf))
geno(header(vcf))


###################################################
### code chunk number 6: readVcf_rowRanges
###################################################
head(rowRanges(vcf), 3)


###################################################
### code chunk number 7: readVcf_fixed
###################################################
ref(vcf)[1:5]
qual(vcf)[1:5]


###################################################
### code chunk number 8: readVcf_ALT
###################################################
alt(vcf)[1:5]


###################################################
### code chunk number 9: geno_hdr
###################################################
geno(vcf)
sapply(geno(vcf), class)


###################################################
### code chunk number 10: explore_geno
###################################################
geno(header(vcf))["DS",]


###################################################
### code chunk number 11: dim_geno
###################################################
DS <-geno(vcf)$DS
dim(DS)
DS[1:3,]


###################################################
### code chunk number 12: fivenum
###################################################
fivenum(DS)


###################################################
### code chunk number 13: DS_zero
###################################################
length(which(DS==0))/length(DS)


###################################################
### code chunk number 14: DS_hist
###################################################
hist(DS[DS != 0], breaks=seq(0, 2, by=0.05), 
    main="DS non-zero values", xlab="DS")


###################################################
### code chunk number 15: info
###################################################
info(vcf)[1:4, 1:5]


###################################################
### code chunk number 16: examine_dbSNP
###################################################
library(SNPlocs.Hsapiens.dbSNP.20101109)
rd <- rowRanges(vcf)
seqlevels(rd) <- "ch22"
ch22snps <- getSNPlocs("ch22")
dbsnpchr22 <- sub("rs", "", names(rd)) %in% ch22snps$RefSNP_id
table(dbsnpchr22)


###################################################
### code chunk number 17: header_info
###################################################
info(header(vcf))[c("VT", "LDAF", "RSQ"),]


###################################################
### code chunk number 18: examine_quality
###################################################
metrics <- data.frame(QUAL=qual(vcf), inDbSNP=dbsnpchr22,
    VT=info(vcf)$VT, LDAF=info(vcf)$LDAF, RSQ=info(vcf)$RSQ)


###################################################
### code chunk number 19: examine_ggplot2
###################################################
library(ggplot2)
ggplot(metrics, aes(x=RSQ, fill=inDbSNP)) +
    geom_density(alpha=0.5) +
    scale_x_continuous(name="MaCH / Thunder Imputation Quality") +
    scale_y_continuous(name="Density") +
    theme(legend.position="top")


###################################################
### code chunk number 20: subset_ranges
###################################################
rng <- GRanges(seqnames="22", ranges=IRanges(
           start=c(50301422, 50989541), 
           end=c(50312106, 51001328),
           names=c("gene_79087", "gene_644186")))


###################################################
### code chunk number 21: subset_TabixFile
###################################################
tab <- TabixFile(fl)
vcf_rng <- readVcf(tab, "hg19", param=rng)


###################################################
### code chunk number 22: VariantAnnotation.Rnw:211-212
###################################################
head(rowRanges(vcf_rng), 3)


###################################################
### code chunk number 23: subset_scanVcfHeader
###################################################
hdr <- scanVcfHeader(fl)
## e.g., INFO and GENO fields
head(info(hdr), 3)
head(geno(hdr), 3)


###################################################
### code chunk number 24: subset_ScanVcfParam
###################################################
## Return all 'fixed' fields, "LAF" from 'info' and "GT" from 'geno'
svp <- ScanVcfParam(info="LDAF", geno="GT")
vcf1 <- readVcf(fl, "hg19", svp)
names(geno(vcf1))


###################################################
### code chunk number 25: subset_ScanVcfParam_new
###################################################
svp_all <- ScanVcfParam(info="LDAF", geno="GT", which=rng) 
svp_all


###################################################
### code chunk number 26: locate_rename_seqlevels
###################################################
library(TxDb.Hsapiens.UCSC.hg19.knownGene)
txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene
seqlevels(vcf) <- "chr22"
rd <- rowRanges(vcf)
loc <- locateVariants(rd, txdb, CodingVariants())
head(loc, 3)


###################################################
### code chunk number 27: AllVariants (eval = FALSE)
###################################################
## allvar <- locateVariants(rd, txdb, AllVariants())


###################################################
### code chunk number 28: locate_gene_centric
###################################################
## Did any coding variants match more than one gene?
splt <- split(mcols(loc)$GENEID, mcols(loc)$QUERYID) 
table(sapply(splt, function(x) length(unique(x)) > 1))

## Summarize the number of coding variants by gene ID.
splt <- split(mcols(loc)$QUERYID, mcols(loc)$GENEID)
head(sapply(splt, function(x) length(unique(x))), 3)


###################################################
### code chunk number 29: predictCoding
###################################################
library(BSgenome.Hsapiens.UCSC.hg19)
coding <- predictCoding(vcf, txdb, seqSource=Hsapiens)
coding[5:7]


###################################################
### code chunk number 30: predictCoding_frameshift
###################################################
## CONSEQUENCE is 'frameshift' where translation is not possible
coding[mcols(coding)$CONSEQUENCE == "frameshift"]


###################################################
### code chunk number 31: nonsynonymous
###################################################
nms <- names(coding)
idx <- mcols(coding)$CONSEQUENCE == "nonsynonymous"
nonsyn <- coding[idx]
names(nonsyn) <- nms[idx]
rsids <- unique(names(nonsyn)[grep("rs", names(nonsyn), fixed=TRUE)])


###################################################
### code chunk number 32: polyphen
###################################################
library(PolyPhen.Hsapiens.dbSNP131)

pp <- select(PolyPhen.Hsapiens.dbSNP131, keys=rsids,
          cols=c("TRAININGSET", "PREDICTION", "PPH2PROB"))
head(pp[!is.na(pp$PREDICTION), ])


###################################################
### code chunk number 33: snpMatrix
###################################################
res <- genotypeToSnpMatrix(vcf)
res


###################################################
### code chunk number 34: snpMatrix_ALT
###################################################
allele2 <- res$map[["allele.2"]]
## number of alternate alleles per variant
unique(elementNROWS(allele2))


###################################################
### code chunk number 35: VariantAnnotation.Rnw:449-464
###################################################
fl.gl <- system.file("extdata", "gl_chr1.vcf", package="VariantAnnotation") 
vcf.gl <- readVcf(fl.gl, "hg19")
geno(vcf.gl)

## Convert the "GL" FORMAT field to a SnpMatrix
res <- genotypeToSnpMatrix(vcf.gl, uncertain=TRUE)
res
t(as(res$genotype, "character"))[c(1,3,7), 1:5]

## Compare to a SnpMatrix created from the "GT" field
res.gt <- genotypeToSnpMatrix(vcf.gl, uncertain=FALSE)
t(as(res.gt$genotype, "character"))[c(1,3,7), 1:5]

## What are the original likelihoods for rs58108140?
geno(vcf.gl)$GL["rs58108140", 1:5]


###################################################
### code chunk number 36: writeVcf
###################################################
fl <- system.file("extdata", "ex2.vcf", package="VariantAnnotation")

out1.vcf <- tempfile()
out2.vcf <- tempfile()
in1 <- readVcf(fl, "hg19")
writeVcf(in1, out1.vcf)
in2 <- readVcf(out1.vcf, "hg19")
writeVcf(in2, out2.vcf)
in3 <- readVcf(out2.vcf, "hg19")

identical(rowRanges(in1), rowRanges(in3))
identical(geno(in1), geno(in2))


###################################################
### code chunk number 37: sessionInfo
###################################################
sessionInfo()