inst/doc/motifbreakR-vignette.R

In motifbreakR: A Package For Predicting The Disruptiveness Of Single Nucleotide Polymorphisms On Transcription Factor Binding Sites

## ----style, echo = FALSE, results='hide', message=FALSE, warning=FALSE----------------------------
BiocStyle::markdown()
options(width=100)
knitr::opts_chunk$set(cache=TRUE, autodep=TRUE)

## ----start, cache=TRUE, message=FALSE-------------------------------------------------------------
library(motifbreakR)
pca.snps.file <- system.file("extdata", "pca.enhancer.snps", package = "motifbreakR")
pca.snps <- as.character(read.table(pca.snps.file)[,1])

## ----outline, eval=FALSE--------------------------------------------------------------------------
#  variants <- snps.from.rsid(rsid = pca.snps,
#                             dbSNP = SNPlocs.Hsapiens.dbSNP142.GRCh37,
#                             search.genome = BSgenome.Hsapiens.UCSC.hg19)
#  motifbreakr.results <- motifbreakR(snpList = variants, pwmList = MotifDb, threshold = 0.9)
#  plotMB(results = motifbreakr.results, rsid = "rs7837328", effect = "strong")

## ----whichsnps, message=FALSE, cache=TRUE---------------------------------------------------------
library(BSgenome)
available.SNPs()

## ----showbed, message=FALSE, cache=TRUE-----------------------------------------------------------
snps.file <- system.file("extdata", "snps.bed", package = "motifbreakR")
read.delim(snps.file, header = FALSE)

## ----fromrsid, message=FALSE, cache=TRUE----------------------------------------------------------
library(SNPlocs.Hsapiens.dbSNP142.GRCh37) # dbSNP142 in hg19
library(BSgenome.Hsapiens.UCSC.hg19)     # hg19 genome
head(pca.snps)
snps.mb <- snps.from.rsid(rsid = pca.snps,
                          dbSNP = SNPlocs.Hsapiens.dbSNP142.GRCh37,
                          search.genome = BSgenome.Hsapiens.UCSC.hg19)
snps.mb

## ----whichgenome, message=FALSE, cache=TRUE-------------------------------------------------------
library(BSgenome)
genomes <- available.genomes()
length(genomes)
genomes

## ----showbed2, message=TRUE, cache=TRUE-----------------------------------------------------------
snps.bed.file <- system.file("extdata", "snps.bed", package = "motifbreakR")
# see the contents
read.table(snps.bed.file, header = FALSE)

## ----getbedrs, message=FALSE, echo = TRUE, eval = FALSE-------------------------------------------
#  library(SNPlocs.Hsapiens.dbSNP142.GRCh37)
#  #import the BED file
#  snps.mb.frombed <- snps.from.file(file = snps.bed.file,
#                                    dbSNP = SNPlocs.Hsapiens.dbSNP142.GRCh37,
#                                    search.genome = BSgenome.Hsapiens.UCSC.hg19,
#                                    format = "bed")
#  snps.mb.frombed

## ----getbedrs.t, message=TRUE, echo = FALSE, eval = TRUE, cache=TRUE------------------------------
library(SNPlocs.Hsapiens.dbSNP142.GRCh37)
example.snpfrombed <- system.file("extdata", "example.snpfrombed.rda", package = "motifbreakR")
load(example.snpfrombed)
message("Warning message:
In snps.from.file(file = snps.bed.file, dbSNP = SNPlocs.Hsapiens.dbSNP142.GRCh37:
  7601289 was found as a match for chr2:12594018:G:A; using entry from dbSNP")
snps.mb.frombed

## ----getbedcust, message=FALSE, cache=TRUE--------------------------------------------------------
library(BSgenome.Drerio.UCSC.danRer7)
snps.bedfile.nors <- system.file("extdata", "danRer.bed", package = "motifbreakR")
read.table(snps.bedfile.nors, header = FALSE)
snps.mb.frombed <- snps.from.file(file = snps.bedfile.nors,
                                  search.genome = BSgenome.Drerio.UCSC.danRer7,
                                  format = "bed")
snps.mb.frombed

## ----getindel, message=FALSE, cache=TRUE----------------------------------------------------------
snps.indel.vcf <- system.file("extdata", "chek2.vcf.gz", package = "motifbreakR")
snps.indel <- variants.from.file(file = snps.indel.vcf, 
                                 search.genome = BSgenome.Hsapiens.UCSC.hg19, 
                                 format = "vcf")
snps.indel

## ----showindel, message=FALSE, cache=TRUE---------------------------------------------------------
snps.indel[nchar(snps.indel$REF) > 1 | nchar(snps.indel$ALT) > 1]

## ----motifdb, message=FALSE, cache=TRUE-----------------------------------------------------------
library(MotifDb)
MotifDb

## ----motifdbtableshow, message=FALSE, cache=TRUE, eval=FALSE--------------------------------------
#  ### Here we can see which organisms are availible under which sources
#  ### in MotifDb
#  table(mcols(MotifDb)$organism, mcols(MotifDb)$dataSource)

## ----motifdbtableres, message=FALSE, cache=TRUE, echo=FALSE---------------------------------------
knitr::kable(table(mcols(MotifDb)$organism, mcols(MotifDb)$dataSource), format = "html", table.attr="class=\"table table-striped table-hover\"")

## ----motifbreakrmot, message=FALSE, cache=TRUE----------------------------------------------------
data(motifbreakR_motif)
motifbreakR_motif

## ----hocomocomot, message=FALSE, cache=TRUE-------------------------------------------------------
data(hocomoco)
hocomoco

## ----runmotifbreakr, message=TRUE, cache=TRUE-----------------------------------------------------
results <- motifbreakR(snpList = snps.mb[1:5], filterp = TRUE,
                       pwmList = hocomoco,
                       threshold = 1e-4,
                       method = "ic",
                       bkg = c(A=0.25, C=0.25, G=0.25, T=0.25),
                       BPPARAM = BiocParallel::bpparam())

## ----getonesnp, cache=TRUE------------------------------------------------------------------------
rs1006140 <- results[names(results) %in% "rs1006140"]
rs1006140

## ----calcp.show, eval=FALSE-----------------------------------------------------------------------
#  rs1006140 <- calculatePvalue(rs1006140)
#  rs1006140

## ----calcp, cache=TRUE, echo=FALSE----------------------------------------------------------------
example.p <- system.file("extdata", "example.pvalue.rda", package = "motifbreakR")
load(example.p)
rs1006140

## ----aboutparallel, message=TRUE, cache=TRUE------------------------------------------------------
BiocParallel::registered()
BiocParallel::bpparam()

## ----plotting, cache=TRUE, fig.retina=2, fig.align='center', fig.height=8, fig.width=6------------
plotMB(results = results, rsid = "rs1006140", effect = "strong")

## ---- cache=TRUE----------------------------------------------------------------------------------
sessionInfo()