### R code from vignette source 'vignettes/VariantAnnotation/inst/doc/VariantAnnotation.Rnw' ################################################### ### code chunk number 1: options ################################################### options(width=72) ################################################### ### code chunk number 2: readVcf ################################################### library(VariantAnnotation) fl <- system.file("extdata", "chr22.vcf.gz", package="VariantAnnotation") vcf <- readVcf(fl, "hg19") vcf ################################################### ### code chunk number 3: readVcf_showheader ################################################### hdr <- exptData(vcf)[["header"]] hdr ################################################### ### code chunk number 4: headeraccessors ################################################### fixed(hdr) head(info(hdr), 3) ################################################### ### code chunk number 5: readVcf_rowData ################################################### head(rowData(vcf)) ################################################### ### code chunk number 6: readVcf_fixed ################################################### head(fixed(vcf), 3) qual(vcf)[1:3,] ################################################### ### code chunk number 7: readVcf_ALT ################################################### alternate <- values(alt(vcf))[["ALT"]] alternate ## number of ALT values per variant unique(elementLengths(alternate)) head(unlist(alternate)) ################################################### ### code chunk number 8: readVcf_info_header ################################################### head(info(hdr), 4) ################################################### ### code chunk number 9: readVcf_info_data ################################################### info(vcf)[1:3, 2:3] ################################################### ### code chunk number 10: VariantAnnotation.Rnw:137-138 ################################################### geno(hdr) ################################################### ### code chunk number 11: readVCF_geno ################################################### geno(vcf) geno(vcf)$GT[1:3,1:5] geno(vcf)$DS[1:3,1:5] ################################################### ### code chunk number 12: subset_ranges ################################################### rng <- GRanges(seqnames="22", ranges=IRanges(c(50301422, 50989541), c(50312106, 51001328))) names(rng) <- c("gene_79087", "gene_644186") ################################################### ### code chunk number 13: subset_TabixFile ################################################### tab <- TabixFile(fl) tab ################################################### ### code chunk number 14: subset_call ################################################### vcf_rng <- readVcf(tab, "hg19", rng) vcf_rng ################################################### ### code chunk number 15: VariantAnnotation.Rnw:189-190 ################################################### rowData(vcf_rng) ################################################### ### code chunk number 16: subset_scanVcfHeader ################################################### hdr <- scanVcfHeader(fl) hdr ################################################### ### code chunk number 17: subset_infoformat ################################################### ## INFO fields info_DF <- info(hdr) rownames(info_DF) ## FORMAT fields geno_DF <- geno(hdr) rownames(geno_DF) ################################################### ### code chunk number 18: subset_elements ################################################### info_DF[rownames(info_DF) == "LDAF", ] geno_DF[rownames(geno_DF) == "GT", ] ################################################### ### code chunk number 19: subset_ScanVcfParam ################################################### ## Return "ALT" from 'fixed', "LAF" from 'info' and "GT" from 'geno' svp <- ScanVcfParam(fixed="ALT", info="LDAF", geno="GT") ## Return all 'fixed' fields, "LAF" from 'info' and "GT" from 'geno' svp <- ScanVcfParam(info="LDAF", geno="GT") svp ################################################### ### code chunk number 20: subset_ScanVcfParam2 ################################################### vcf_flds <- readVcf(fl, "hg19", svp) geno(vcf_flds) head(info(vcf_flds), 3) ################################################### ### code chunk number 21: subset_ScanVcfParam_new ################################################### svp_all <- ScanVcfParam(info="LDAF", geno="GT", which=rng) svp_all ################################################### ### code chunk number 22: subset_both ################################################### readVcf(tab, "hg19", svp_all) ################################################### ### code chunk number 23: seqlevels_rd ################################################### rowdat <- rowData(vcf) seqlevels(rowdat) ################################################### ### code chunk number 24: seqlevels_TxDb ################################################### library(TxDb.Hsapiens.UCSC.hg19.knownGene) txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene head(seqlevels(txdb)) ################################################### ### code chunk number 25: seqlevels_rename ################################################### ## rename variant seqlevels in the VCF object vcf <- renameSeqlevels(vcf, c("22"="chr22")) ## extract the rowData with modified seqlevels rd <- rowData(vcf) ## confirm seqlevels are the same intersect(seqlevels(rd), seqlevels(txdb)) ################################################### ### code chunk number 26: CodingVariants ################################################### loc <- locateVariants(rd, txdb, CodingVariants()) head(loc, 4) ################################################### ### code chunk number 27: SpliceSiteVariants ################################################### head(locateVariants(rd, txdb, SpliceSiteVariants()), 4) ################################################### ### code chunk number 28: AllVariants (eval = FALSE) ################################################### ## allvar <- locateVariants(rd, txdb, AllVariants()) ################################################### ### code chunk number 29: locatVariants_example ################################################### ## Did any coding variants match more than one gene? table(sapply(split(values(loc)[["geneID"]], values(loc)[["queryID"]]), function(x) length(unique(x)) > 1)) ## Summarize the number of coding variants by gene ID idx <- sapply(split(values(loc)[["queryID"]], values(loc)[["geneID"]]), unique) sapply(idx, length) ################################################### ### code chunk number 30: predictCoding ################################################### library(BSgenome.Hsapiens.UCSC.hg19) coding <- predictCoding(vcf, txdb, seqSource=Hsapiens) coding[5:9] ################################################### ### code chunk number 31: predictCoding_frameshift ################################################### ## consequence is 'frameshift' where translation is not possible coding[values(coding)[["consequence"]] == "frameshift"] ################################################### ### code chunk number 32: nonsynonymous ################################################### nms <- names(coding) idx <- values(coding)[["consequence"]] == "nonsynonymous" nonsyn <- coding[idx] names(nonsyn) <- nms[idx] rsids <- unique(names(nonsyn)[grep("rs", names(nonsyn), fixed=TRUE)]) ################################################### ### code chunk number 33: sift ################################################### library(SIFT.Hsapiens.dbSNP132) ## rsids in the package head(keys(SIFT.Hsapiens.dbSNP132)) ## list available columns cols(SIFT.Hsapiens.dbSNP132) ## select a subset of columns ## a warning is thrown when a key is not found in the database subst <- c("RSID", "PREDICTION", "SCORE", "AACHANGE", "PROTEINID") sift <- select(SIFT.Hsapiens.dbSNP132, keys=rsids, cols=subst) head(sift) ################################################### ### code chunk number 34: polyphen ################################################### library(PolyPhen.Hsapiens.dbSNP131) inSIFT <- unique(sift$RSID[!is.na(sift$PREDICTION)]) pp <- select(PolyPhen.Hsapiens.dbSNP131, keys=inSIFT, cols=c("TRAININGSET", "PREDICTION", "PPH2PROB")) head(pp[!is.na(pp$PREDICTION), ]) ################################################### ### code chunk number 35: snpMatrix ################################################### calls <- geno(vcf)$GT a0 <- values(ref(vcf))[["REF"]] a1 <- values(alt(vcf))[["ALT"]] res <- MatrixToSnpMatrix(calls, a0, a1) res ################################################### ### code chunk number 36: snpMatrix_ALT ################################################### allele2 <- res$map[["allele.2"]] ## number of alternate alleles per variant unique(elementLengths(allele2)) unlist(allele2) ################################################### ### code chunk number 37: longform_header ################################################### rownames(info_DF) ################################################### ### code chunk number 38: expand_info ################################################### param <- ScanVcfParam(fixed="ALT", info="HOMSEQ") gr <- readVcfLongForm(fl, "hg19", param) head(gr) ################################################### ### code chunk number 39: 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(in2, in3) ################################################### ### code chunk number 40: sessionInfo ################################################### sessionInfo()