1. Download dbSNP and convert from VCF to BED with BEDOPS vcf2bed.

   $ vcf2bed < (gunzip -c dbSNP.vcf.gz) > dbSNP.bed
   

2. Download phyloP44 or other alignments and convert from WIG to BED with BEDOPS wig2bed.

   $ wig2bed < (gunzip -c phyloP.wig.gz) > phyloP.bed
   

   File names will depend on what you download from NCBI and UCSC Goldenpath.

3. Map signal (score) to SNPs with BEDOPS bedmap:

   $ bedmap --echo --echo-map-score --skip-unmapped --delim '\t' dbSNP.bed phyloP.bed > answer.bed
   

4. Read answer.bed into R with read.table() and find the population of signals/scores for all mapped variants. You should be able to get the minimum-scoring variants from that distribution. Or use sort -nr on the signal column and head to get the minimum score, and awk to filter answer.bed for variants with that score.

In terms of published data, A. Quinlans group performed a nice analysis and released the data properly for others to use. Kudos.

Data are in BED and Bigwig format so hopefully useful. You can perhaps take the inverse of the constrained regions and restrict to exons to find the regions you want (use bedtools or similar?).

https://github.com/quinlan-lab/ccrhtml