### R code from vignette source 'D6_RegulatoryVariants.Rnw' ################################################### ### code chunk number 1: setup ################################################### library(Morgan2013) library(org.Hs.eg.db) library(TxDb.Hsapiens.UCSC.hg19.knownGene) library(BSgenome.Hsapiens.UCSC.hg19) library(seqnames.db) library(gmapR) library(VariantTools) library(lattice) library(MotifDb) library(parallel) ################################################### ### code chunk number 2: parallel ################################################### library(parallel); options(mc.cores=detectCores()) ################################################### ### code chunk number 3: annotation-libraries ################################################### library(org.Hs.eg.db) library(TxDb.Hsapiens.UCSC.hg19.knownGene) library(BSgenome.Hsapiens.UCSC.hg19) ################################################### ### code chunk number 4: path-to-bam ################################################### path <- "~/BigData/TERT/bam" ################################################### ### code chunk number 5: sanity-check ################################################### stopifnot(length(dir(path, "bam$")) == 11) stopifnot(length(dir(path, "bai$")) == 11) ################################################### ### code chunk number 6: consistent-names ################################################### seqnameStyle(BSgenome.Hsapiens.UCSC.hg19) <- seqnameStyle(TxDb.Hsapiens.UCSC.hg19.knownGene) <- "NCBI" ################################################### ### code chunk number 7: anno-symbols ################################################### bsgenome <- BSgenome.Hsapiens.UCSC.hg19 txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene geneid <- select(org.Hs.eg.db, "TERT", "ENTREZID", "SYMBOL") ################################################### ### code chunk number 8: VariantTools-setup ################################################### VARIANTTOOLS_OK <- local({ sys <- !Sys.info()["sysname"] %in% c("Darwin", "Windows") sys && require(VariantTools) && require(gmapR) }) stopifnot(VARIANTTOOLS_OK) ################################################### ### code chunk number 9: mini-genome ################################################### roi <- GRanges("5", IRanges(1, seqlengths(bsgenome)[["5"]])) chr5seq <- getSeq(bsgenome, roi) names(chr5seq) <- "5" ################################################### ### code chunk number 10: make-mini ################################################### library(gmapR) genome.5 <- GmapGenome(chr5seq, name="hg19_5", create=TRUE) ################################################### ### code chunk number 11: TERT-promoters ################################################### rng <- transcriptsBy(txdb, "gene")[[geneid$ENTREZID]] ################################################### ### code chunk number 12: TERT-starts ################################################### unique(start(rng)) ################################################### ### code chunk number 13: TERT-promoter ################################################### unique(start(rng)) pregion <- promoters(range(rng), upstream=330, downstream=0) seqlevels(pregion) <- "5" ################################################### ### code chunk number 14: setup-variant-calling ################################################### library(VariantTools) vtparam <- VariantTallyParam(genome.5, readlen=101L, which=pregion, indels=TRUE) ################################################### ### code chunk number 15: call-variants ################################################### fls <- dir(path, "bam$", full=TRUE) called <- lapply(fls, callVariants, tally.param=vtparam) ################################################### ### code chunk number 16: call-variants-tidy ################################################### len <- elementLengths(called) called <- do.call(c, called) id <- sub(".TERT.bam", "", basename(fls)) called$id <- factor(rep(id, len), levels=unique(id)) ################################################### ### code chunk number 17: load-called ################################################### data(called, package="Morgan2013") ################################################### ### code chunk number 18: explore-variants ################################################### library(lattice) altCounts <- xtabs(cycleCount.10.91 ~ start(called) + id, mcols(called)) plt <- levelplot(altCounts, xlab="Position", ylab=NULL, scales=list(x=list(rot=45)), aspect="fill", col.regions=rev(gray.colors(100, 0, 1))) ################################################### ### code chunk number 19: variant-levelplot ################################################### pdf("variant-levelplot.pdf", width=10, height=3) print(plt) xx <- dev.off() ################################################### ### code chunk number 20: snp-def ################################################### snp <- GRanges("5", IRanges(1295228, width=1)) snp <- flank(snp, 10, both=TRUE) ################################################### ### code chunk number 21: ref-seq ################################################### refSeq <- getSeq(bsgenome, snp) ################################################### ### code chunk number 22: variantSequences ################################################### library(Morgan2013) variantSequences altSeq <- variantSequences(called, snp, refSeq) ################################################### ### code chunk number 23: alt-seq ################################################### altConsensus <- DNAStringSet(consensusString(altSeq)) ################################################### ### code chunk number 24: complement ################################################### complement(c(refSeq, altConsensus)) ################################################### ### code chunk number 25: pwms ################################################### library(MotifDb) idx <- with(mcols(MotifDb), organism=="Hsapiens" & dataSource == "JASPAR_CORE") jasparHumanPWMs <- MotifDb[idx] ################################################### ### code chunk number 26: pwm-match ################################################### ## matchPWMs: helper in Morgan2013 minScore = "90%" # high min. matching score: 1 base change (refHits <- matchPWMs(jasparHumanPWMs, reverseComplement(refSeq)[[1]], minScore)) (altHits <- matchPWMs(jasparHumanPWMs, reverseComplement(altConsensus)[[1]], minScore))