motifbreakR: an Introduction
Simon G. Coetzee, Gerhard A. Coetzee and Dennis J. Hazelett
2019-05-02
1 Introduction
This document offers an introduction and overview of motifBreakR, which allows the biologist to judge whether the sequence surrounding a polymorphism or mutation is a good match to known transcription factor binding sites, and how much information is gained or lost in one allele of the polymorphism relative to another or mutation vs. wildtype. motifBreakR is flexible, giving a choice of algorithms for interrogation of genomes with motifs from public sources that users can choose from; these are 1) a weighted-sum, 2) log-probabilities, and 3) relative entropy. motifBreakR can predict effects for novel or previously described variants in public databases, making it suitable for tasks beyond the scope of its original design. Lastly, it can be used to interrogate any genome curated within Bioconductor.
motifBreakR works with position probability matrices (PPM). PPM are derived as the fractional occurrence of nucleotides A,C,G, and T at each position of a position frequency matrix (PFM). PFM are simply the tally of each nucleotide at each position across a set of aligned sequences. With a PPM, one can generate probabilities based on the genome, or more practically, create any number of position specific scoring matrices (PSSM) based on the principle that the PPM contains information about the likelihood of observing a particular nucleotide at a particular position of a true transcription factor binding site.
This guide includes a brief overview of the processing flow, an example focusing more in depth on the practical aspect of using motifBreakR, and finally a detailed section on the scoring methods employed by the package.
2 Processing overview
motifBreakR may be used to interrogate SNPs or SNVs for their potential effect on transcription factor binding by examining how the two alleles of the variant effect the binding score of a motif. The basic process is outlined in the figure below.
2.1 Outline of process
motifbreakR workflow: How inputs (trapezoids) are generated from R functions (rectangles). Diamonds represent decisions of the user.
This straightforward process allows the interrogation of SNPs and SNVs in the context of the different species represented by BSgenome packages (at least 22 different species) and allows the use of the full MotifDb data set that includes over 4200 motifs across 8 studies and 22 organisms that we have supplemented with over 2800 additional motifs across four additional studies in Humans see data(encodemotif)1, data(factorbook)2, data(hocomoco)3 and data(homer)4 for the additional studies that we have included.
Practically motifBreakR has involves three phases.
- Read in Single Nucleotide Variants: The first step is to read in the list of variants. Variants can be a list of rsIDs if your SNPs are represented in one of the SNPlocs packages, may be included as a .bed file with a specifically formatted name field, or may be provided as a .vcf file. We then transform these input such that it may be read my motifBreakR.
- Find Broken Motifs: Next we present motifBreakR with the input generated in the previous step, along with a set of motifs formatted as class
MotifList, and your preferred scoring method. - Visualize SNPs and motifs: Finally we can visualize which motifs are broken by any individual SNP using our plotting function.
3 How To Use motifbreakR: A Practical Example
This section offers an example of how to use motifBreakR to identify potentially disrupted transcription factor binding sites due to 701 SNPs output from a FunciSNP analysis of Prostate Cancer (PCa) genome wide association studies (GWAS) risk loci. The SNPs are included in this package here:
library(motifbreakR)
pca.snps.file <- system.file("extdata", "pca.enhancer.snps", package = "motifbreakR")
pca.snps <- as.character(read.table(pca.snps.file)[,1])The simplest form of a motifBreakR analysis is summarized as follows:
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")Lets look at these steps more closely and see how we can customize our analysis.
3.1 Step 1 | Read in Single Nucleotide Variants
Variants can be input either as a list of rsIDs or as a .bed file. The main factor determining which you will use is if your variants have rsIDs that are included in one of the Bioconductor SNPlocs packages. The present selection is seen here:
library(BSgenome)
available.SNPs()## [1] "SNPlocs.Hsapiens.dbSNP.20101109" "SNPlocs.Hsapiens.dbSNP.20120608"
## [3] "SNPlocs.Hsapiens.dbSNP141.GRCh38" "SNPlocs.Hsapiens.dbSNP142.GRCh37"
## [5] "SNPlocs.Hsapiens.dbSNP144.GRCh37" "SNPlocs.Hsapiens.dbSNP144.GRCh38"
## [7] "SNPlocs.Hsapiens.dbSNP149.GRCh38" "SNPlocs.Hsapiens.dbSNP150.GRCh38"
## [9] "SNPlocs.Hsapiens.dbSNP151.GRCh38" "XtraSNPlocs.Hsapiens.dbSNP141.GRCh38"
## [11] "XtraSNPlocs.Hsapiens.dbSNP144.GRCh37" "XtraSNPlocs.Hsapiens.dbSNP144.GRCh38"For cases where your rsIDs are not available in a SNPlocs package, or you have novel variants that are not cataloged at all, variants may be entered in BED format as seen below:
snps.file <- system.file("extdata", "snps.bed", package = "motifbreakR")
read.delim(snps.file, header = FALSE)## V1 V2 V3 V4 V5 V6
## 1 chr2 12581137 12581138 rs10170896 0 +
## 2 chr2 12594017 12594018 chr2:12594018:G:A 0 +
## 3 chr3 192388677 192388678 rs13068005 0 +
## 4 chr4 122361479 122361480 rs12644995 0 +
## 5 chr6 44503245 44503246 chr6:44503246:A:T 0 +
## 6 chr6 44503247 44503248 chr6:44503248:G:C 0 +
## 7 chr6 85232897 85232898 rs4510639 0 +
## 8 chr6 44501872 44501873 rs932680 0 +Our requirements for the BED file are that it must include chromosome, start, end, name, score and strand fields – where the name field is required to be in one of two formats, either an rsID that is present in a SNPlocs package, or in the form chromosome:position:referenceAllele:alternateAllele e.g., chr2:12594018:G:A. It is also essential that the fields are TAB separated, not a mixture of tabs and spaces.
More to the point here are the two methods for reading in the variants.
3.1.1 SNPs from rsID:
We use the SNPlocs.Hsapiens.dbSNP142.GRCh37 which is the SNP locations and alleles defined in dbSNP142 as a source for looking up our rsIDs and BSgenome.Hsapiens.UCSC.hg19 which holds the reference sequence for UCSC genome build hg19. Additional SNPlocs packages are availble from Bioconductor.
library(SNPlocs.Hsapiens.dbSNP142.GRCh37) # dbSNP142 in hg19
library(BSgenome.Hsapiens.UCSC.hg19) # hg19 genome
head(pca.snps)## [1] "rs1551515" "rs1551513" "rs17762938" "rs4473999" "rs7823297" "rs9656964"snps.mb <- snps.from.rsid(rsid = pca.snps,
dbSNP = SNPlocs.Hsapiens.dbSNP142.GRCh37,
search.genome = BSgenome.Hsapiens.UCSC.hg19)## Warning in rowids2rowidx(user_rowids, ids, x_rowids, ifnotfound): SNP ids not found: rs78914317, rs75425437, rs114099824, rs79509278, rs74738513
##
## They were dropped.snps.mb## GRanges object with 700 ranges and 4 metadata columns:
## seqnames ranges strand | SNP_id alleles_as_ambig REF ALT
## <Rle> <IRanges> <Rle> | <character> <DNAStringSet> <DNAStringSet> <DNAStringSet>
## rs10007915 chr4 106065308 * | rs10007915 S C G
## rs10015716 chr4 95548550 * | rs10015716 R G A
## rs10034824 chr4 95524838 * | rs10034824 K G T
## rs10056823 chr5 115609454 * | rs10056823 S C G
## rs1006140 chr19 38778915 * | rs1006140 R A G
## ... ... ... ... . ... ... ... ...
## rs9901746 chr17 36103149 * | rs9901746 R G A
## rs9908087 chr17 69106937 * | rs9908087 K T G
## rs991429 chr17 69109773 * | rs991429 R G A
## rs9973650 chr2 238380266 * | rs9973650 R G A
## rs998071 chr4 95591976 * | rs998071 S C G
## -------
## seqinfo: 93 sequences (1 circular) from hg19 genomeA far greater variety of variants may be read into motifBreakR via the snps.from.file function. In fact motifBreakR will work with any organism present as a Bioconductor BSgenome package. This includes 76 genomes representing 22 species.
library(BSgenome)
genomes <- available.genomes()
length(genomes)## [1] 93genomes## [1] "BSgenome.Alyrata.JGI.v1" "BSgenome.Amellifera.BeeBase.assembly4"
## [3] "BSgenome.Amellifera.UCSC.apiMel2" "BSgenome.Amellifera.UCSC.apiMel2.masked"
## [5] "BSgenome.Aofficinalis.NCBI.V1" "BSgenome.Athaliana.TAIR.04232008"
## [7] "BSgenome.Athaliana.TAIR.TAIR9" "BSgenome.Btaurus.UCSC.bosTau3"
## [9] "BSgenome.Btaurus.UCSC.bosTau3.masked" "BSgenome.Btaurus.UCSC.bosTau4"
## [11] "BSgenome.Btaurus.UCSC.bosTau4.masked" "BSgenome.Btaurus.UCSC.bosTau6"
## [13] "BSgenome.Btaurus.UCSC.bosTau6.masked" "BSgenome.Btaurus.UCSC.bosTau8"
## [15] "BSgenome.Carietinum.NCBI.v1" "BSgenome.Celegans.UCSC.ce10"
## [17] "BSgenome.Celegans.UCSC.ce11" "BSgenome.Celegans.UCSC.ce2"
## [19] "BSgenome.Celegans.UCSC.ce6" "BSgenome.Cfamiliaris.UCSC.canFam2"
## [21] "BSgenome.Cfamiliaris.UCSC.canFam2.masked" "BSgenome.Cfamiliaris.UCSC.canFam3"
## [23] "BSgenome.Cfamiliaris.UCSC.canFam3.masked" "BSgenome.Cjacchus.UCSC.calJac3"
## [25] "BSgenome.Dmelanogaster.UCSC.dm2" "BSgenome.Dmelanogaster.UCSC.dm2.masked"
## [27] "BSgenome.Dmelanogaster.UCSC.dm3" "BSgenome.Dmelanogaster.UCSC.dm3.masked"
## [29] "BSgenome.Dmelanogaster.UCSC.dm6" "BSgenome.Drerio.UCSC.danRer10"
## [31] "BSgenome.Drerio.UCSC.danRer5" "BSgenome.Drerio.UCSC.danRer5.masked"
## [33] "BSgenome.Drerio.UCSC.danRer6" "BSgenome.Drerio.UCSC.danRer6.masked"
## [35] "BSgenome.Drerio.UCSC.danRer7" "BSgenome.Drerio.UCSC.danRer7.masked"
## [37] "BSgenome.Ecoli.NCBI.20080805" "BSgenome.Gaculeatus.UCSC.gasAcu1"
## [39] "BSgenome.Gaculeatus.UCSC.gasAcu1.masked" "BSgenome.Ggallus.UCSC.galGal3"
## [41] "BSgenome.Ggallus.UCSC.galGal3.masked" "BSgenome.Ggallus.UCSC.galGal4"
## [43] "BSgenome.Ggallus.UCSC.galGal4.masked" "BSgenome.Ggallus.UCSC.galGal5"
## [45] "BSgenome.Hsapiens.1000genomes.hs37d5" "BSgenome.Hsapiens.NCBI.GRCh38"
## [47] "BSgenome.Hsapiens.UCSC.hg17" "BSgenome.Hsapiens.UCSC.hg17.masked"
## [49] "BSgenome.Hsapiens.UCSC.hg18" "BSgenome.Hsapiens.UCSC.hg18.masked"
## [51] "BSgenome.Hsapiens.UCSC.hg19" "BSgenome.Hsapiens.UCSC.hg19.masked"
## [53] "BSgenome.Hsapiens.UCSC.hg38" "BSgenome.Hsapiens.UCSC.hg38.masked"
## [55] "BSgenome.Mdomestica.UCSC.monDom5" "BSgenome.Mfascicularis.NCBI.5.0"
## [57] "BSgenome.Mfuro.UCSC.musFur1" "BSgenome.Mmulatta.UCSC.rheMac2"
## [59] "BSgenome.Mmulatta.UCSC.rheMac2.masked" "BSgenome.Mmulatta.UCSC.rheMac3"
## [61] "BSgenome.Mmulatta.UCSC.rheMac3.masked" "BSgenome.Mmulatta.UCSC.rheMac8"
## [63] "BSgenome.Mmusculus.UCSC.mm10" "BSgenome.Mmusculus.UCSC.mm10.masked"
## [65] "BSgenome.Mmusculus.UCSC.mm8" "BSgenome.Mmusculus.UCSC.mm8.masked"
## [67] "BSgenome.Mmusculus.UCSC.mm9" "BSgenome.Mmusculus.UCSC.mm9.masked"
## [69] "BSgenome.Osativa.MSU.MSU7" "BSgenome.Ptroglodytes.UCSC.panTro2"
## [71] "BSgenome.Ptroglodytes.UCSC.panTro2.masked" "BSgenome.Ptroglodytes.UCSC.panTro3"
## [73] "BSgenome.Ptroglodytes.UCSC.panTro3.masked" "BSgenome.Ptroglodytes.UCSC.panTro5"
## [75] "BSgenome.Ptroglodytes.UCSC.panTro6" "BSgenome.Rnorvegicus.UCSC.rn4"
## [77] "BSgenome.Rnorvegicus.UCSC.rn4.masked" "BSgenome.Rnorvegicus.UCSC.rn5"
## [79] "BSgenome.Rnorvegicus.UCSC.rn5.masked" "BSgenome.Rnorvegicus.UCSC.rn6"
## [81] "BSgenome.Scerevisiae.UCSC.sacCer1" "BSgenome.Scerevisiae.UCSC.sacCer2"
## [83] "BSgenome.Scerevisiae.UCSC.sacCer3" "BSgenome.Sscrofa.UCSC.susScr11"
## [85] "BSgenome.Sscrofa.UCSC.susScr3" "BSgenome.Sscrofa.UCSC.susScr3.masked"
## [87] "BSgenome.Tgondii.ToxoDB.7.0" "BSgenome.Tguttata.UCSC.taeGut1"
## [89] "BSgenome.Tguttata.UCSC.taeGut1.masked" "BSgenome.Tguttata.UCSC.taeGut2"
## [91] "BSgenome.Vvinifera.URGI.IGGP12Xv0" "BSgenome.Vvinifera.URGI.IGGP12Xv2"
## [93] "BSgenome.Vvinifera.URGI.IGGP8X"3.1.2 SNPs from BED formatted file:
Here we examine two possibilities. In one case we have a mixture of rsIDs and our naming scheme that allows for arbitrary variants. Second we have a list of variants for the zebrafish Danio rerio that does not have a SNPlocs package, but does have it’s genome present among the availible.genomes().
snps.bed.file <- system.file("extdata", "snps.bed", package = "motifbreakR")
# see the contents
read.table(snps.bed.file, header = FALSE)## V1 V2 V3 V4 V5 V6
## 1 chr2 12581137 12581138 rs10170896 0 +
## 2 chr2 12594017 12594018 chr2:12594018:G:A 0 +
## 3 chr3 192388677 192388678 rs13068005 0 +
## 4 chr4 122361479 122361480 rs12644995 0 +
## 5 chr6 44503245 44503246 chr6:44503246:A:T 0 +
## 6 chr6 44503247 44503248 chr6:44503248:G:C 0 +
## 7 chr6 85232897 85232898 rs4510639 0 +
## 8 chr6 44501872 44501873 rs932680 0 +Seeing as we have some SNPs listed by their rsIDs we can query those by including a SNPlocs object as an argument to snps.from.file
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## 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## GRanges object with 8 ranges and 4 metadata columns:
## seqnames ranges strand | SNP_id alleles_as_ambig REF
## <Rle> <IRanges> <Rle> | <character> <DNAStringSet> <DNAStringSet>
## rs10170896 chr2 12581138 + | rs10170896 R G
## rs12644995 chr4 122361480 + | rs12644995 M C
## rs13068005 chr3 192388678 + | rs13068005 R G
## rs4510639 chr6 85232898 + | rs4510639 Y C
## rs932680 chr6 44501873 + | rs932680 K G
## 7601289 chr2 12594018 + | 7601289 R G
## chr6:44503246:A:T chr6 44503246 + | chr6:44503246:A:T W A
## chr6:44503248:G:C chr6 44503248 + | chr6:44503248:G:C S G
## ALT
## <DNAStringSet>
## rs10170896 A
## rs12644995 A
## rs13068005 A
## rs4510639 T
## rs932680 T
## 7601289 A
## chr6:44503246:A:T T
## chr6:44503248:G:C C
## -------
## seqinfo: 93 sequences (1 circular) from hg19 genomeWe see also that one of our custom variants chr2:12594018:G:A was actually already included in dbSNP, and was therefor annotated in the output as rs7601289
If our BED file includes no rsIDs, then we may omit the dbSNP argument from the function. This example uses variants from Danio rerio
library(BSgenome.Drerio.UCSC.danRer7)
snps.bedfile.nors <- system.file("extdata", "danRer.bed", package = "motifbreakR")
read.table(snps.bedfile.nors, header = FALSE)## V1 V2 V3 V4 V5 V6
## 1 chr18 13030932 13030933 chr18:13030933:G:A 0 +
## 2 chr18 30445455 30445456 chr18:30445456:T:A 0 +
## 3 chr5 22065023 22065024 chr5:22065024:A:T 0 +
## 4 chr14 36140941 36140942 chr14:36140942:T:A 0 +
## 5 chr3 16701576 16701577 chr3:16701577:T:A 0 +
## 6 chr14 20887995 20887996 chr14:20887996:G:A 0 +
## 7 chr7 25195449 25195450 chr7:25195450:G:T 0 +
## 8 chr2 59181852 59181853 chr2:59181853:A:G 0 +
## 9 chr3 58162674 58162675 chr3:58162675:C:T 0 +
## 10 chr22 18708824 18708825 chr22:18708825:T:A 0 +snps.mb.frombed <- snps.from.file(file = snps.bedfile.nors,
search.genome = BSgenome.Drerio.UCSC.danRer7,
format = "bed")
snps.mb.frombed## GRanges object with 10 ranges and 4 metadata columns:
## seqnames ranges strand | SNP_id alleles_as_ambig REF
## <Rle> <IRanges> <Rle> | <character> <DNAStringSet> <DNAStringSet>
## chr18:13030933:G:A chr18 13030933 + | chr18:13030933:G:A R G
## chr18:30445456:T:A chr18 30445456 + | chr18:30445456:T:A W T
## chr5:22065024:A:T chr5 22065024 + | chr5:22065024:A:T W A
## chr14:36140942:T:A chr14 36140942 + | chr14:36140942:T:A W T
## chr3:16701577:T:A chr3 16701577 + | chr3:16701577:T:A W T
## chr14:20887996:G:A chr14 20887996 + | chr14:20887996:G:A R G
## chr7:25195450:G:T chr7 25195450 + | chr7:25195450:G:T K G
## chr2:59181853:A:G chr2 59181853 + | chr2:59181853:A:G R A
## chr3:58162675:C:T chr3 58162675 + | chr3:58162675:C:T Y C
## chr22:18708825:T:A chr22 18708825 + | chr22:18708825:T:A W T
## ALT
## <DNAStringSet>
## chr18:13030933:G:A A
## chr18:30445456:T:A A
## chr5:22065024:A:T T
## chr14:36140942:T:A A
## chr3:16701577:T:A A
## chr14:20887996:G:A A
## chr7:25195450:G:T T
## chr2:59181853:A:G G
## chr3:58162675:C:T T
## chr22:18708825:T:A A
## -------
## seqinfo: 26 sequences (1 circular) from danRer7 genomesnps.from.file also can take as input a vcf file with SNVs, by using format = "vcf".
3.2 Step 2 | Find Broken Motifs
Now that we have our data in the required format, we may continue to the task at hand, and determine which variants modify potential transcription factor binding. An important element of this task is identifying a set of transcription factor binding motifs that we wish to query. Fortunately MotifDb includes a large selection of motifs across multiple species that we can see here:
library(MotifDb)
MotifDb## MotifDb object of length 9933
## | Created from downloaded public sources: 2013-Aug-30
## | 9933 position frequency matrices from 14 sources:
## | FlyFactorSurvey: 614
## | HOCOMOCOv10: 1066
## | HOMER: 332
## | JASPAR_2014: 592
## | JASPAR_CORE: 459
## | ScerTF: 196
## | SwissRegulon: 684
## | UniPROBE: 380
## | cisbp_1.02: 874
## | hPDI: 437
## | jaspar2016: 1209
## | jaspar2018: 1564
## | jolma2013: 843
## | stamlab: 683
## | 61 organism/s
## | Hsapiens: 4616
## | Mmusculus: 1411
## | Dmelanogaster: 1287
## | Scerevisiae: 1051
## | Athaliana: 803
## | Celegans: 90
## | other: 675
## Scerevisiae-ScerTF-ABF2-badis
## Scerevisiae-ScerTF-CAT8-badis
## Scerevisiae-ScerTF-CST6-badis
## Scerevisiae-ScerTF-ECM23-badis
## Scerevisiae-ScerTF-EDS1-badis
## ...
## Mmusculus-UniPROBE-Zfp740.UP00022
## Mmusculus-UniPROBE-Zic1.UP00102
## Mmusculus-UniPROBE-Zic2.UP00057
## Mmusculus-UniPROBE-Zic3.UP00006
## Mmusculus-UniPROBE-Zscan4.UP00026### Here we can see which organisms are availible under which sources
### in MotifDb
table(mcols(MotifDb)$organism, mcols(MotifDb)$dataSource)| FlyFactorSurvey | HOCOMOCOv10 | HOMER | JASPAR_2014 | JASPAR_CORE | ScerTF | SwissRegulon | UniPROBE | cisbp_1.02 | hPDI | jaspar2016 | jaspar2018 | jolma2013 | stamlab | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Acarolinensis | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 0 | 0 | 0 |
| Amajus | 0 | 0 | 0 | 3 | 3 | 0 | 0 | 0 | 0 | 0 | 3 | 3 | 0 | 0 |
| Anidulans | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 8 | 0 | 0 | 0 | 0 | 0 |
| Apisum | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 0 | 0 | 0 |
| Aterreus | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 0 | 0 | 0 |
| Athaliana | 0 | 0 | 0 | 48 | 5 | 0 | 0 | 0 | 107 | 0 | 191 | 452 | 0 | 0 |
| Bdistachyon | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 2 | 0 | 0 | 0 | 0 | 0 |
| Celegans | 0 | 0 | 0 | 15 | 5 | 0 | 0 | 2 | 22 | 0 | 23 | 23 | 0 | 0 |
| Cparvum | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 1 | 0 | 0 | 0 | 0 | 0 |
| Csativa | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 2 | 0 | 0 | 0 | 0 | 0 |
| Ddiscoideum | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 9 | 0 | 0 | 0 | 0 | 0 |
| Dmelanogaster | 614 | 0 | 0 | 131 | 125 | 0 | 0 | 0 | 138 | 0 | 139 | 140 | 0 | 0 |
| Drerio | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 3 | 0 | 0 | 0 | 0 | 0 |
| Gaculeatus | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 0 | 0 | 0 |
| Gallus | 0 | 0 | 0 | 1 | 2 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| Ggallus | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 4 | 4 | 0 | 0 |
| Hcapsulatum | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 0 | 0 | 0 |
| Hroretzi | 0 | 0 | 0 | 1 | 1 | 0 | 0 | 0 | 0 | 0 | 1 | 1 | 0 | 0 |
| Hsapiens | 0 | 640 | 0 | 117 | 66 | 0 | 684 | 2 | 313 | 437 | 442 | 522 | 710 | 683 |
| Hvulgare | 0 | 0 | 0 | 1 | 1 | 0 | 0 | 0 | 0 | 0 | 1 | 1 | 0 | 0 |
| Mdomestica | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 0 | 0 | 0 |
| Mgallopavo | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 0 | 0 | 0 |
| Mmurinus | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 0 | 0 | 0 |
| Mmusculus | 0 | 426 | 0 | 66 | 47 | 0 | 0 | 282 | 132 | 0 | 165 | 160 | 133 | 0 |
| Mmusculus;Hsapiens | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 3 | 0 | 0 |
| Mmusculus;Rnorvegicus | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 2 | 0 | 0 |
| Mmusculus;Rnorvegicus;Hsapiens | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 6 | 0 | 0 |
| Mmusculus;Rnorvegicus;Hsapiens;Ocuniculus | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0 |
| Mmusculus;Rnorvegicus;Omykiss;Ggallus;Hsapiens | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0 |
| Mmusculus;Rnorvegicus;Xlaevis;Stropicalis;Ggallus;Hsapiens;Btaurus;Ocuniculus | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 2 | 0 | 0 |
| Mmusculus;Rrattus;Hsapiens;Ocuniculus | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0 |
| Mtruncatula | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 0 | 0 | 0 |
| NA | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 40 | 0 | 0 |
| Ncrassa | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 15 | 0 | 0 | 0 | 0 | 0 |
| Ngruberi | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 0 | 0 | 0 |
| Nhaematococca | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 0 | 0 | 0 |
| Nsp. | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 1 | 0 | 0 |
| Nsylvestris | 0 | 0 | 0 | 1 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| Nvectensis | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 0 | 0 | 0 |
| Ocuniculus | 0 | 0 | 0 | 1 | 1 | 0 | 0 | 0 | 0 | 0 | 1 | 1 | 0 | 0 |
| Osativa | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 5 | 0 | 0 | 0 | 0 | 0 |
| Otauri | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 2 | 0 | 0 | 0 | 0 | 0 |
| Pcapensis | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 0 | 0 | 0 |
| Pfalciparum | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 2 | 26 | 0 | 0 | 0 | 0 | 0 |
| Phybrida | 0 | 0 | 0 | 1 | 1 | 0 | 0 | 0 | 0 | 0 | 1 | 1 | 0 | 0 |
| Ppatens | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 7 | 0 | 0 | 0 | 0 | 0 |
| Ppygmaeus | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 0 | 0 | 0 |
| Psativum | 0 | 0 | 0 | 3 | 3 | 0 | 0 | 0 | 1 | 0 | 3 | 3 | 0 | 0 |
| Ptetraurelia | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 0 | 0 | 0 |
| Rnorvegicus | 0 | 0 | 0 | 6 | 8 | 0 | 0 | 0 | 2 | 0 | 12 | 7 | 0 | 0 |
| Rnorvegicus;Hsapiens | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0 |
| Rrattus | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 2 | 1 | 0 | 0 |
| Scerevisiae | 0 | 0 | 0 | 177 | 177 | 196 | 0 | 91 | 60 | 0 | 175 | 175 | 0 | 0 |
| Taestivam | 0 | 0 | 0 | 1 | 1 | 0 | 0 | 0 | 0 | 0 | 1 | 1 | 0 | 0 |
| Tthermophila | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 0 | 0 | 0 |
| Vertebrata | 0 | 0 | 0 | 12 | 4 | 0 | 0 | 0 | 0 | 0 | 1 | 1 | 0 | 0 |
| Vvinifera | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 0 | 0 | 0 |
| Xlaevis | 0 | 0 | 0 | 1 | 1 | 0 | 0 | 0 | 0 | 0 | 3 | 2 | 0 | 0 |
| Xtropicalis | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 0 | 0 | 0 |
| Zmays | 0 | 0 | 0 | 6 | 6 | 0 | 0 | 0 | 1 | 0 | 6 | 8 | 0 | 0 |
We have leveraged the MotifList introduced by MotifDb to include an additional set of motifs that have been gathered to this package:
data(motifbreakR_motif)
motifbreakR_motif## MotifDb object of length 2817
## | Created from downloaded public sources: 2013-Aug-30
## | 2817 position frequency matrices from 4 sources:
## | ENCODE-motif: 2065
## | FactorBook: 79
## | HOCOMOCO: 426
## | HOMER: 247
## | 1 organism/s
## | Hsapiens: 2817
## Hsapiens-ENCODE-motifs-SIX5_disc1
## Hsapiens-ENCODE-motifs-MYC_disc1
## Hsapiens-ENCODE-motifs-SRF_disc1
## Hsapiens-ENCODE-motifs-AP1_disc1
## Hsapiens-ENCODE-motifs-SIX5_disc2
## ...
## Hsapiens-HOMER-yy1.motif
## Hsapiens-HOMER-zbtb33.motif
## Hsapiens-HOMER-zfx.motif
## Hsapiens-HOMER-znf263.motif
## Hsapiens-HOMER-znf711.motifThe different studies included in this data set may be called individually; for example:
data(hocomoco)
hocomoco## MotifDb object of length 426
## | Created from downloaded public sources: 2013-Aug-30
## | 426 position frequency matrices from 1 source:
## | HOCOMOCO: 426
## | 1 organism/s
## | Hsapiens: 426
## Hsapiens-HOCOMOCO-AHR_si
## Hsapiens-HOCOMOCO-AIRE_f2
## Hsapiens-HOCOMOCO-ALX1_si
## Hsapiens-HOCOMOCO-ANDR_do
## Hsapiens-HOCOMOCO-AP2A_f2
## ...
## Hsapiens-HOCOMOCO-ZN333_f1
## Hsapiens-HOCOMOCO-ZN350_f1
## Hsapiens-HOCOMOCO-ZN384_f1
## Hsapiens-HOCOMOCO-ZN423_f1
## Hsapiens-HOCOMOCO-ZN589_f1See ?motifbreakR_motif for more information and citations.
Some of our data sets include a sequenceCount. These include FlyFactorSurvey, hPDI, JASPAR_2014, JASPAR_CORE, and jolma2013 from MotifDb and HOCOMOCO from the set of motifbreakR_motif. Using these we calculate a pseudocount to allow us to calculate the logarithms in the case where we have a 0 in a pwm. The calculation for incorporating pseudocounts is ppm <- (ppm * sequenceCount + 0.25)/(sequenceCount + 1). If the sequenceCount for a particular ppm is NA we use 20 as a default sequenceCount.
Now that we have the three necessary components to run motifBreakR:
- A
BSgenomeobject for our organism, in this caseBSgenome.Hsapiens.UCSC.hg19 - A
MotifListobject containing our motifs, in this casehocomoco, - And our input
GRangesobject generated bysnps.from.rsid, in this casesnps.mb
We get to the task of actually running the function motifbreakR().
We have several options that we may pass to the function, the main ones that will dictate how long the function will run with a given set of variants and motifs are the threshold we pass and the method we use to score.
Here we specify the snpList, pwmList, threshold that we declare as the cutoff for reporting results, filterp set to true declares that we are filtering by p-value, the method, and bkg the relative nucleotide frequency of A, C, G, and T.
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())The results reveal which variants disrupt which motifs, and to which degree. If we want to examine a single variant, we can select one like this:
rs1006140 <- results[names(results) %in% "rs1006140"]
rs1006140## GRanges object with 11 ranges and 18 metadata columns:
## seqnames ranges strand | REF ALT snpPos motifPos
## <Rle> <IRanges> <Rle> | <DNAStringSet> <DNAStringSet> <integer> <integer>
## rs1006140 chr19 38778913-38778923 - | A G 38778915 9
## rs1006140 chr19 38778913-38778923 - | A G 38778915 9
## rs1006140 chr19 38778915-38778925 - | A G 38778915 11
## rs1006140 chr19 38778912-38778924 - | A G 38778915 10
## rs1006140 chr19 38778910-38778920 - | A G 38778915 6
## rs1006140 chr19 38778905-38778926 - | A G 38778915 12
## rs1006140 chr19 38778907-38778925 - | A G 38778915 11
## rs1006140 chr19 38778910-38778920 - | A G 38778915 6
## rs1006140 chr19 38778907-38778923 - | A G 38778915 9
## rs1006140 chr19 38778905-38778924 - | A G 38778915 10
## rs1006140 chr19 38778912-38778928 - | A G 38778915 14
## geneSymbol dataSource providerName providerId seqMatch
## <character> <character> <character> <character> <character>
## rs1006140 EGR1 HOCOMOCO EGR1_f2 EGR1_HUMAN ccAcccaccct
## rs1006140 EGR2 HOCOMOCO EGR2_si EGR2_HUMAN ccAcccaccct
## rs1006140 EGR4 HOCOMOCO EGR4_f1 EGR4_HUMAN Acccaccctcc
## rs1006140 EPAS1 HOCOMOCO EPAS1_si EPAS1_HUMAN cccAcccaccctc
## rs1006140 KLF1 HOCOMOCO KLF1_f1 KLF1_HUMAN tccccAcccac
## rs1006140 RREB1 HOCOMOCO RREB1_si RREB1_HUMAN cccactccccAcccaccctcct
## rs1006140 SP1 HOCOMOCO SP1_f2 SP1_HUMAN cactccccAcccaccctcc
## rs1006140 SP1 HOCOMOCO SP1_f1 SP1_HUMAN tccccAcccac
## rs1006140 SP2 HOCOMOCO SP2_si SP2_HUMAN cactccccAcccaccct
## rs1006140 SP3 HOCOMOCO SP3_f1 SP3_HUMAN cccactccccAcccaccctc
## rs1006140 ZBTB4 HOCOMOCO ZBTB4_si ZBTB4_HUMAN cccAcccaccctcctgt
## pctRef pctAlt scoreRef scoreAlt Refpvalue
## <numeric> <numeric> <numeric> <numeric> <logical>
## rs1006140 0.872576362913107 0.919856011083053 9.92842685969322 10.466388405248 <NA>
## rs1006140 0.960333809386458 0.996981439092749 10.0425988795408 10.4258379589411 <NA>
## rs1006140 0.904466067177726 0.915572992725868 7.42410931066732 7.5152780480774 <NA>
## rs1006140 0.918315516234323 0.79870400254378 7.09686011612808 6.17248699389117 <NA>
## rs1006140 0.962133797324208 0.808767865451189 8.99202486538536 7.5586792363881 <NA>
## rs1006140 0.825601779803683 0.716284210979446 10.0684130278547 8.73525888375256 <NA>
## rs1006140 0.900129993865763 0.927682614187614 10.092293503227 10.4012145845887 <NA>
## rs1006140 0.907157686157437 0.940596439155853 9.55675262677591 9.90902422787688 <NA>
## rs1006140 0.850346403843652 0.905003499448131 6.21074088981548 6.60994415222112 <NA>
## rs1006140 0.917522805258346 0.948382954812404 8.96608373907422 9.26765083202962 <NA>
## rs1006140 0.803054830846981 0.764173920551312 12.7442409323818 12.127212468764 <NA>
## Altpvalue alleleRef alleleAlt effect
## <logical> <numeric> <numeric> <character>
## rs1006140 <NA> 0.155509026321165 0.748436131852517 weak
## rs1006140 <NA> 0.233333333333333 0.716666666666667 weak
## rs1006140 <NA> 0.100612612234726 0.530474476237945 weak
## rs1006140 <NA> 0.869570911818731 0 strong
## rs1006140 <NA> 1 0 strong
## rs1006140 <NA> 0.935589712108114 0 strong
## rs1006140 <NA> 0.148364998792684 0.647582503097021 weak
## rs1006140 <NA> 0.106401451920105 0.644070946496345 weak
## rs1006140 <NA> 0 0.619047619047619 weak
## rs1006140 <NA> 0.0873493975903614 0.627560240963857 weak
## rs1006140 <NA> 0.799044082390768 0.200955917609232 weak
## -------
## seqinfo: 93 sequences (1 circular) from hg19 genomeHere we can see that SNP rs1006140 disrupts multiple motifs. We can then check what the pvalue for each allele is with regard to each motif, using calculatePvalue.
rs1006140 <- calculatePvalue(rs1006140)
rs1006140## GRanges object with 11 ranges and 18 metadata columns:
## seqnames ranges strand | REF ALT snpPos motifPos
## <Rle> <IRanges> <Rle> | <DNAStringSet> <DNAStringSet> <integer> <integer>
## rs1006140 chr19 38778913-38778923 - | A G 38778915 9
## rs1006140 chr19 38778913-38778923 - | A G 38778915 9
## rs1006140 chr19 38778915-38778925 - | A G 38778915 11
## rs1006140 chr19 38778912-38778924 - | A G 38778915 10
## rs1006140 chr19 38778910-38778920 - | A G 38778915 6
## rs1006140 chr19 38778905-38778926 - | A G 38778915 12
## rs1006140 chr19 38778907-38778925 - | A G 38778915 11
## rs1006140 chr19 38778910-38778920 - | A G 38778915 6
## rs1006140 chr19 38778907-38778923 - | A G 38778915 9
## rs1006140 chr19 38778905-38778924 - | A G 38778915 10
## rs1006140 chr19 38778912-38778928 - | A G 38778915 14
## geneSymbol dataSource providerName providerId seqMatch
## <character> <character> <character> <character> <character>
## rs1006140 EGR1 HOCOMOCO EGR1_f2 EGR1_HUMAN ccAcccaccct
## rs1006140 EGR2 HOCOMOCO EGR2_si EGR2_HUMAN ccAcccaccct
## rs1006140 EGR4 HOCOMOCO EGR4_f1 EGR4_HUMAN Acccaccctcc
## rs1006140 EPAS1 HOCOMOCO EPAS1_si EPAS1_HUMAN cccAcccaccctc
## rs1006140 KLF1 HOCOMOCO KLF1_f1 KLF1_HUMAN tccccAcccac
## rs1006140 RREB1 HOCOMOCO RREB1_si RREB1_HUMAN cccactccccAcccaccctcct
## rs1006140 SP1 HOCOMOCO SP1_f2 SP1_HUMAN cactccccAcccaccctcc
## rs1006140 SP1 HOCOMOCO SP1_f1 SP1_HUMAN tccccAcccac
## rs1006140 SP2 HOCOMOCO SP2_si SP2_HUMAN cactccccAcccaccct
## rs1006140 SP3 HOCOMOCO SP3_f1 SP3_HUMAN cccactccccAcccaccctc
## rs1006140 ZBTB4 HOCOMOCO ZBTB4_si ZBTB4_HUMAN cccAcccaccctcctgt
## pctRef pctAlt scoreRef scoreAlt
## <numeric> <numeric> <numeric> <numeric>
## rs1006140 0.872576362913107 0.919856011083053 9.92842685969322 10.466388405248
## rs1006140 0.960333809386458 0.996981439092749 10.0425988795408 10.4258379589411
## rs1006140 0.904466067177726 0.915572992725868 7.42410931066732 7.5152780480774
## rs1006140 0.918315516234323 0.79870400254378 7.09686011612808 6.17248699389117
## rs1006140 0.962133797324208 0.808767865451189 8.99202486538536 7.5586792363881
## rs1006140 0.825601779803683 0.716284210979446 10.0684130278547 8.73525888375256
## rs1006140 0.900129993865763 0.927682614187614 10.092293503227 10.4012145845887
## rs1006140 0.907157686157437 0.940596439155853 9.55675262677591 9.90902422787688
## rs1006140 0.850346403843652 0.905003499448131 6.21074088981548 6.60994415222112
## rs1006140 0.917522805258346 0.948382954812404 8.96608373907422 9.26765083202962
## rs1006140 0.803054830846981 0.764173920551312 12.7442409323818 12.127212468764
## Refpvalue Altpvalue alleleRef alleleAlt
## <numeric> <numeric> <numeric> <numeric>
## rs1006140 0.000141620635986328 3.98159027099609e-05 0.155509026321165 0.748436131852517
## rs1006140 2.88486480712891e-05 2.14576721191406e-06 0.233333333333333 0.716666666666667
## rs1006140 6.91413879394531e-05 5.60283660888672e-05 0.100612612234726 0.530474476237945
## rs1006140 1.02818012237549e-05 0.000711917877197266 0.869570911818731 0
## rs1006140 3.24249267578125e-05 0.0010833740234375 1 0
## rs1006140 1.83549503276481e-05 0.000875384435971682 0.935589712108114 0
## rs1006140 2.50729608524125e-05 6.10251299804077e-06 0.148364998792684 0.647582503097021
## rs1006140 4.69684600830078e-05 1.23977661132812e-05 0.106401451920105 0.644070946496345
## rs1006140 0.000467563804704696 3.93299269489944e-05 0 0.619047619047619
## rs1006140 1.95318170881364e-05 2.19532830669777e-06 0.0873493975903614 0.627560240963857
## rs1006140 1.21764605864882e-05 5.19903842359781e-05 0.799044082390768 0.200955917609232
## effect
## <character>
## rs1006140 weak
## rs1006140 weak
## rs1006140 weak
## rs1006140 strong
## rs1006140 strong
## rs1006140 strong
## rs1006140 weak
## rs1006140 weak
## rs1006140 weak
## rs1006140 weak
## rs1006140 weak
## -------
## seqinfo: 93 sequences (1 circular) from hg19 genomeAnd here we see that for each SNP we have at least one allele achieving a p-value below 1e-4 threshold that we required. The seqMatch column shows what the reference genome sequence is at that location, with the variant position appearing in an uppercase letter. pctRef and pctAlt display the the score for the motif in the sequence as a percentage of the best score that motif could achieve on an ideal sequence. In other words \((scoreVariant-minscorePWM)/(maxscorePWM-minscorePWM)\). We can also see the absolute scores for our method in scoreRef and scoreAlt and thier respective p-values.
3.2.1 Parallel Evaluation (an aside)
Important to note, is that motifbreakR uses the BiocParallel parallel back-end, and one may modify what type of parallel evaluation it uses (or if it runs in parallel at all). Here we can see the versions available on the machine this vignette was compiled on.
BiocParallel::registered()## $MulticoreParam
## class: MulticoreParam
## bpisup: FALSE; bpnworkers: 4; bptasks: 0; bpjobname: BPJOB
## bplog: FALSE; bpthreshold: INFO; bpstopOnError: TRUE
## bpRNGseed: ; bptimeout: 2592000; bpprogressbar: FALSE
## bpexportglobals: TRUE
## bplogdir: NA
## bpresultdir: NA
## cluster type: FORK
##
## $SnowParam
## class: SnowParam
## bpisup: FALSE; bpnworkers: 4; bptasks: 0; bpjobname: BPJOB
## bplog: FALSE; bpthreshold: INFO; bpstopOnError: TRUE
## bpRNGseed: ; bptimeout: 2592000; bpprogressbar: FALSE
## bpexportglobals: TRUE
## bplogdir: NA
## bpresultdir: NA
## cluster type: SOCK
##
## $SerialParam
## class: SerialParam
## bpisup: FALSE; bpnworkers: 1; bptasks: 0; bpjobname: BPJOB
## bplog: FALSE; bpthreshold: INFO; bpstopOnError: TRUE
## bpRNGseed: ; bptimeout: 2592000; bpprogressbar: FALSE
## bpexportglobals: TRUE
## bplogdir: NA
## bpresultdir: NABiocParallel::bpparam()## class: MulticoreParam
## bpisup: FALSE; bpnworkers: 4; bptasks: 0; bpjobname: BPJOB
## bplog: FALSE; bpthreshold: INFO; bpstopOnError: TRUE
## bpRNGseed: ; bptimeout: 2592000; bpprogressbar: FALSE
## bpexportglobals: TRUE
## bplogdir: NA
## bpresultdir: NA
## cluster type: FORKBy default motifbreakR uses bpparam() as an argument to BPPARAM and will use all available cores on the machine on which it is running. However on Windows machines this reverts to using a serial evaluation model, so if you wish to run in parallel on a Windows machine consider using a different parameter shown in BiocParallel::registered() such as SnowParam passing BPPARAM = SnowParam().
3.3 Step 3 | Visualize
Now that we have our results, we can visualize them with the function plotMB. Lets take a look at rs1006140.
plotMB(results = results, rsid = "rs1006140", effect = "strong")
4 Methods
motifBreakR works with position probability matrices (PPM). PPM are derived as the fractional occurrence of nucleotides A,C,G, and T at each position of a position frequency matrix (PFM). PFM are simply the tally of each nucleotide at each position across a set of aligned sequences. With a PPM, one can generate probabilities based on the genome, or more practically, create any number of position specific scoring matrices (PSSM) based on the principle that the PPM contains information about the likelihood of observing a particular nucleotide at a particular position of a true transcription factor binding site. What follows is a discussion of the three different algorithms that may be employed in calls to the motifBreakR function via the method argument.
Suppose we have a frequency matrix \(M\) of width \(n\) (i.e. a PPM as described above). Furthermore, we have a sequence \(s\) also of length \(n\), such that \(s_{i} \in \{ A,T,C,G \}, i = 1,\ldots n\). Each column of \(M\) contains the frequencies of each letter in each position.
Commonly in the literature sequences are scored as the sum of log probabilities:
4.1 Sum of Log Probabilities
\[F( s,M ) = \sum_{i = 1}^{n}{\log( \frac{M_{s_{i},i}}{b_{s_{i}}} )}\]
where \(b_{s_{i}}\) is the background frequency of letter \(s_{i}\) in the genome of interest. This method can be specified by the user as method='log'.
As an alternative to this method, we introduced a scoring method to directly weight the score by the importance of the position within the match sequence. This method of weighting is accessed by specifying method='default' (weighted sum). A general representation of this scoring method is given by:
4.2 Weighted Sum
\[F( s,M ) = p( s ) \cdot \omega_{M}\]
where \(p_{s}\) is the scoring vector derived from sequence \(s\) and matrix \(M\), and \(w_{M}\) is a weight vector derived from \(M\). First, we compute the scoring vector of position scores \(p\):
4.3 Calculating Scoring Vector \(p\)
\[p( s ) = ( M_{s_{i},i} ) \textrm{ where } \frac{i = 1,\ldots n}{s_{i} \in \{ A,C,G,T \}}\]
and second, for each \(M\) a constant vector of weights \(\omega_{M} = ( \omega_{1},\omega_{2},\ldots,\omega_{n} )\).
There are two methods for producing \(\omega_{M}\). The first, which we call weighted sum, is the difference in the probabilities for the two letters of the polymorphism (or variant), i.e. \(\Delta p_{s_{i}}\), or the difference of the maximum and minimum values for each column of \(M\):
4.4 Calculating \(\omega\) For Weighted Sum
\[\omega_{i} = \max \{ M_{i} \} - \min \{ M_{i} \}\textrm{ where }i = 1,\ldots n\]
The second variation of this theme is to weight by relative entropy. Thus the relative entropy weight for each column \(i\) of the matrix is given by:
4.5 Calculating \(\omega\) For Relative Entropy
\[\omega_{i} = \sum_{j \in \{ A,C,G,T \}}^{}{M_{j,i}\log_2( \frac{M_{j,i}}{b_{i}} )}\textrm{ where }i = 1,\ldots n\]
where \(b_{i}\) is again the background frequency of the letter \(i\).
Thus, there are 3 possible algorithms to apply via the method argument. The first is the standard summation of log probabilities (method='log'). The second and third are the weighted sum and information content methods (method='default' and method='ic') specified by equations for Weighted Sum and Relative Entropy, respectively. motifBreakR assumes a uniform background nucleotide distribution (\(b\)) in equations 4.1 and 4.5 unless otherwise specified by the user. Since we are primarily interested in the difference between alleles, background frequency is not a major factor, although it can change the results. Additionally, inclusion of background frequency introduces potential bias when collections of motifs are employed, since motifs are themselves unbalanced with respect to nucleotide composition. With these cautions in mind, users may override the uniform distribution if so desired. For all three methods, motifBreakR scores and reports the reference and alternate alleles of the sequence (\(F( s_{\textrm{REF}},M )\) and \(F( s_{\textrm{ALT}},M )\)), and provides the matrix scores \(p_{s_{\textrm{REF}}}\) and \(p_{s_{\textrm{ALT}}}\) of the SNP (or variant). The scores are scaled as a fraction of scoring range 0-1 of the motif matrix, \(M\). If either of \(F( s_{\textrm{REF}},M )\) and \(F( s_{\textrm{ALT}},M )\) is greater than a user-specified threshold (default value of 0.85) the SNP is reported. By default motifBreakR does not display neutral effects, (\(\Delta p_{i} < 0.4\)) but this behaviour can be overridden.
4.6 Calculate P-values for PWM match
Additionally, now, with the use of TFMPvalue, we may filter by p-value of the match. This is unfortunately a two step process. First, by invoking filterp=TRUE and setting a threshold at a desired p-value e.g 1e-4, we perform a rough filter on the results by rounding all values in the PWM to two decimal place, and calculating a scoring threshold based upon that. The second step is to use the function calculatePvalue() on a selection of results which will change the Refpvalue and Altpvalue columns in the output from NA to the p-value calculated by TFMsc2pv. This can be (although not always) a very memory and time intensive process if the algorithm doesn’t converge rapidly.
5 Session Info
sessionInfo()## R version 3.6.0 (2019-04-26)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 18.04.2 LTS
##
## Matrix products: default
## BLAS: /home/biocbuild/bbs-3.9-bioc/R/lib/libRblas.so
## LAPACK: /home/biocbuild/bbs-3.9-bioc/R/lib/libRlapack.so
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C LC_TIME=en_US.UTF-8
## [4] LC_COLLATE=C LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C LC_ADDRESS=C
## [10] LC_TELEPHONE=C LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] stats4 parallel grid stats graphics grDevices utils datasets methods
## [10] base
##
## other attached packages:
## [1] BSgenome.Drerio.UCSC.danRer7_1.4.0 BSgenome.Hsapiens.UCSC.hg19_1.4.0
## [3] SNPlocs.Hsapiens.dbSNP142.GRCh37_0.99.5 BSgenome_1.52.0
## [5] rtracklayer_1.44.0 GenomicRanges_1.36.0
## [7] GenomeInfoDb_1.20.0 motifbreakR_1.14.0
## [9] MotifDb_1.26.0 Biostrings_2.52.0
## [11] XVector_0.24.0 IRanges_2.18.0
## [13] S4Vectors_0.22.0 BiocGenerics_0.30.0
## [15] BiocStyle_2.12.0
##
## loaded via a namespace (and not attached):
## [1] ProtGenerics_1.16.0 bitops_1.0-6 matrixStats_0.54.0
## [4] bit64_0.9-7 RColorBrewer_1.1-2 progress_1.2.0
## [7] httr_1.4.0 backports_1.1.4 tools_3.6.0
## [10] R6_2.4.0 rGADEM_2.32.0 rpart_4.1-15
## [13] Hmisc_4.2-0 seqLogo_1.50.0 DBI_1.0.0
## [16] splitstackshape_1.4.8 Gviz_1.28.0 lazyeval_0.2.2
## [19] colorspace_1.4-1 nnet_7.3-12 ade4_1.7-13
## [22] motifStack_1.28.0 gridExtra_2.3 tidyselect_0.2.5
## [25] prettyunits_1.0.2 grImport2_0.1-4 curl_3.3
## [28] bit_1.1-14 compiler_3.6.0 Biobase_2.44.0
## [31] htmlTable_1.13.1 grImport_0.9-1.1 DelayedArray_0.10.0
## [34] checkmate_1.9.1 scales_1.0.0 stringr_1.4.0
## [37] digest_0.6.18 Rsamtools_2.0.0 foreign_0.8-71
## [40] rmarkdown_1.12 dichromat_2.0-0 base64enc_0.1-3
## [43] jpeg_0.1-8 pkgconfig_2.0.2 htmltools_0.3.6
## [46] highr_0.8 ensembldb_2.8.0 htmlwidgets_1.3
## [49] rlang_0.3.4 rstudioapi_0.10 RSQLite_2.1.1
## [52] BiocParallel_1.18.0 acepack_1.4.1 dplyr_0.8.0.1
## [55] VariantAnnotation_1.30.0 RCurl_1.95-4.12 magrittr_1.5
## [58] GenomeInfoDbData_1.2.1 Formula_1.2-3 Matrix_1.2-17
## [61] Rcpp_1.0.1 munsell_0.5.0 stringi_1.4.3
## [64] yaml_2.2.0 MASS_7.3-51.4 SummarizedExperiment_1.14.0
## [67] zlibbioc_1.30.0 plyr_1.8.4 blob_1.1.1
## [70] crayon_1.3.4 lattice_0.20-38 splines_3.6.0
## [73] GenomicFeatures_1.36.0 hms_0.4.2 knitr_1.22
## [76] pillar_1.3.1 MotIV_1.40.0 codetools_0.2-16
## [79] biomaRt_2.40.0 TFMPvalue_0.0.8 XML_3.98-1.19
## [82] glue_1.3.1 evaluate_0.13 biovizBase_1.32.0
## [85] latticeExtra_0.6-28 data.table_1.12.2 BiocManager_1.30.4
## [88] png_0.1-7 gtable_0.3.0 purrr_0.3.2
## [91] assertthat_0.2.1 ggplot2_3.1.1 xfun_0.6
## [94] AnnotationFilter_1.8.0 survival_2.44-1.1 tibble_2.1.1
## [97] GenomicAlignments_1.20.0 AnnotationDbi_1.46.0 memoise_1.1.0
## [100] cluster_2.0.9Website: encode-motifs Paper: Systematic discovery and characterization of regulatory motifs in ENCODE TF binding experiments↩
Website: Factorbook Paper: Sequence features and chromatin structure around the genomic regions bound by 119 human transcription factors↩
Website: HOCOMOCO Paper: HOCOMOCO: a comprehensive collection of human transcription factor binding sites models↩
Website: Homer Paper: http://www.sciencedirect.com/science/article/pii/S1097276510003667↩