Pattern Matching for Transcription Factor Binding Sites

[This article was first published on rforcancer, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

I’m trying to search for binding sites for the transcription factor MAF (i.e. TFBS for MAF) in the promoter regions of various genes. I initially started out looking at a precomputed database of binding sites MAPPER. However the TFBS models that have been used from TRANSFAC or JASPAR are quite unlike the human binding sites  shown in the recent literature (the sites in the recent literature are palindromic sequences called MAREs).

Unfortunately the palindromic sites shown in the literature do not have alignment or position-weight-matrices (PWM) and they are quite variable. It would appear that these binding sites have not been derived by an unbiased screen i.e. ChIP with an antibody against a genome. Rather they have mutated a known site in a single gene promoter and tested which ones still work. Jeez!

So I cannot create my own alignments or position weight matrices. Instead I will try some cruder pattern matching and test a few mismatch and indel options (possibly with IUPAC coding later). Here is scripting for a first effort:

## I am firstly downloading upstream sequences using BSGenome and the latest human sequence build ##
## this creates subject a so-called DNAStringSet of upstream gene sequences i.e. promoters##
## the query string and its reverse complement - for now: will experiment later - ##
CMARE_complement <- complement(CMARE)
## extract the accession name with a bit of messing about ##
nm=strsplit(nm, '_up')
nm=unlist(lapply(nm, function(x)x[1]))
## this is matching forward and complement CMARE DNAStrings to the promoters DNAStringSet ##
mindexF <- vmatchPattern(CMARE, subject, max.mismatch=1)
mindexR <- vmatchPattern(CMARE_complement, subject, max.mismatch=1)
## I'm using 1 mismatch here but will probably adapt later ##
## create indices whF + whR of the matches
count_indexF <- countIndex(mindexF)
count_indexR <- countIndex(mindexR)
total_matches=  sum(count_indexF)+ sum(count_indexR)  # Total number of matches.
whF <- which(count_indexF > 0)
whR <- which(count_indexR > 0)
## next we extract the data for each and put it in a sort of list view ##
for(i in 1: length(whF))  
result_listF[[i]]= Views(subject[[whF[i]]], mindexF[[whF[i]]])
for(i in 1: length(whR))  
result_listR[[i]]= Views(subject[[whR[i]]], mindexR[[whR[i]]])
## and add the original accesion names to the correct list index
result_list=c(result_listF, result_listR)

Now… I don’t like this script much as it seems a faff to do everything twice for forward and reverse strands… but I was in a hurry an didn’t have time to really understand the object classes and combine them properly. I also never like using for loops but in this case I couldn’t see any way to make it more functional stylee..

Anyway that gives a result something like this:

> result_list[1:2]
  Views on a 1000-letter DNAString subject
    start end width
[1]   716 729    14 [TCCTGACGTCAGCA]
  Views on a 1000-letter DNAString subject
    start end width
[1]   929 942    14 [TGCTGACGTCAGCT]

Next I need to use biomart and Ensembl to try and get gene symbols, chromosomes, start positions, strand etc.. for the promoters. There is an R interface package biomaRt for this service. Once again however Ensembl/biomart is having a nervous breakdown– so I’ll probably continue this tomorrow.

To leave a comment for the author, please follow the link and comment on their blog: rforcancer. offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Never miss an update!
Subscribe to R-bloggers to receive
e-mails with the latest R posts.
(You will not see this message again.)

Click here to close (This popup will not appear again)