Have DNA motifs 6-12bp long, trying to get conservation scores

Question

I have about 200 short nucleotide motifs (6-12 bp in length) from the human genome, and I'm trying to see how conserved they are across vertebrates.

I was thinking that I'd need to make a bed file for each motif that lists all of its occurrences in the human genome. From there, I could map the beds to a bigwig files of PhastCons scores (essentially doing the reverse of what the PhastCons software was designed to do). Does that sound like the best approach?

I'm getting stuck at the step of going from motifs to bed files. I've tried using BLAST to find all occurrences of motifs, but their short length is causing issues.
I've tried messing with the e-value threshold, word size, and filter parameters, but I still don't get any hits.

Is there a work-around for this issue, or should I just rethink my entire approach?

Answer 1

In case you have only ACGT in your motifs

The short motifs make it sound as if you are in the business of looking for a kmer counter. You can choose to use existing software or build your own.

1. Using existing software might be your easiest path. An older post from 2014 will probably give you a first idea what's out there: http://homolog.us/blogs/blog/2014/04/07/kmer-counting-a-2014-recap/ . Note that a couple of algorithms mentioned there got successors, so it is worthwhile digging a bit around. The small kmer size will make most of them usable for your needs.
2. As the maximum size of your kmers is comparatively small (12 nt need 24 bits, i.e., max 16.7 million entries in your kmer table), you should be able to easily roll your own kmer counting in about any language you like and on about any of nowadays computer. The pseudocode section on the Wikipedia entry for kmers will give you first pointers for that. Might be a bit more work, but maybe more flexible depending on your needs.

In case you have IUPAC bases (N, W, etc.) in your motifs

I don't know any any pre-existing software doing what you need. I could imagine that the short motifs make using regular expressions doable for this kind of search, but I may be wrong. Testing this should be easy in a simple script as all major programming languages have modules or libraries for REs. Even if it should take a couple of hours to run on your data set, that would be good enough for a one-off calculation.

Answer 2

To scan motifs in a genome (or database) I would use FIMO which will give you the exact locations of these motifs in your genome.

Once you have the locations, you can use a phastCons bigiwig from UCSC to calculate the basewise conservation scores. However, please remember that phastCons scores are smoothed across windows and it might not be the best metric if you are trying to compare the conservation levels at your motif matching sites as compared to the sequences flanking them.

I wrote a package a while back to do this, including doing de-novo motif discovery. However, it might be an overkill for your use case.

Answer 3

Where transcription factor binding sites are concerned, in our papers we have used per-base phyloP conservation data instead of smoothed phastCons scores.

We use BEDOPS bedmap to map scores over multiple, padded binding sites (BED-formatted) for a given motif model. The phyloP scores are WIG files obtained from UCSC goldenpath and converted to BED via wig2bed.

The resulting matrix of binding site intervals and their per-base scores can be turned into a ranked heatmap or aggregated to determine average per-base conservation for a motif model.

While TF binding sites innately show high information content and therefore high conservation, matrices can be further sorted by score-maps of ChIP-seq or DNaseI-seq tag density.

This can help filter for interesting low-noise, high-signal patterns and is useful when the matrix is windowed and the window includes low-information regions with a few high-information residues outside of the main motif — as is the case with CTCF, for instance.