Not with R, but you could use BEDOPS tools:

• Convert the per-chromosome phyloP WIG files on UCSC's servers to sorted BED using wig2bed.

• Prepare your genomic regions into a sorted BED file containing per-base elements.

• Run bedmap on the phyloP data and your regions to map per-base phyloP scores to your per-base regions.

Retrieving phyloP data and converting it to sorted BED

We can make per-chromosome phyloP BED files and ultimately use these with bedmap or other BEDOPS tools. As an example, we can use the vertebrate hg19 phyloP data located at http://hgdownload-test.cse.ucsc.edu/goldenPath/hg19/phyloP46way/vertebrate/

$ rsync -avz --progress rsync://hgdownload.cse.ucsc.edu/goldenPath/hg19/phyloP46way/vertebrate .

...

$ for fn in `ls vertebrate/*.gz`; \
    do gunzip -c $fn \
        | wig2bed - \
        > ${fn%%.*}.bed; \
    done

Preparing your regions-of-interest

Let's assume your regions-of-interest are already in a four-column (or more) BED file called roi.contiguous.unsorted.bed. To make sure it works with BEDOPS correctly, we sort it and create a file called roi.contiguous.bed:

$ sort-bed roi.contiguous.unsorted.bed > roi.contiguous.bed

Next, we want to split the ranges into per-base elements. We assume there is an ID field in the fourth column. We also assume that the ID for a given region-of-interest is uniquely named. When splitting our regions, we can use this ID as a unique prefix, in order to specifically and uniquely identify with which range the per-base element associates:

$ awk ' \
    { \
         regionChromosome = $1; \
         regionStart = $2; \
         regionStop = $3; \
         regionID = $4; \
         baseIdx = 0; \
         for (baseStart = regionStart; baseStart < regionStop; baseStart++) { \
             baseStop = baseStart + 1; \
             print regionChromosome"\t"baseStart"\t"baseStop"\t"regionID"-"baseIdx; \
             baseIdx++; \
         } \
    }' roi.contiguous.bed > roi.perBase.bed

Run head or less on the roi.perBase.bed file so that you see what this script does. When we map phyloP scores, we will be able to map them base-by-base.

Note: If your BED file does not have an ID field, or region IDs are not unique, it is trivial to use awk to add or replace the ID field with a unique, per-row identifier.

Performing bedmap operations

Finally, we're ready to do the map operation and see which elements associate:

$ for chrLabel in `seq 1 22` X Y; \ 
    do \
        mapFn=/path/to/phyloP/chr${chrLabel}.bed; \
        bedmap --chrom chr${chrLabel} --echo --echo-map-score roi.perBase.bed $mapFn \
            > perBase.chr${chrLabel}.answer.bed; \
    done

The file perBase.chr${chrLabel}.answer.bed (where chr${chrLabel} is replaced with the chromosome name from the per-chromosome phyloP files) contains per-base-split regions-of-interest along with the phyloP scores that associate with each base.

To explain the bedmap statement, the --echo operator prints out the region-of-interest (a per-base element from your original regions-of-interest) and --echo-map-score prints out the phyloP score that associates with that base, for each chromosome.

In the case where no phyloP score maps to a base, a NAN is printed.

You may or may not want to use sed or other tools to replace NAN with a zero, particularly if the presence of a zero is a biologically meaningful score, while the absence of a score is also biologically meaningful, but in a different way.

Some more work to do

A couple more things you could do that might make this easier, if you have to automate this process:

• Concatenate the per-chromosome phyloP BED files into one master BED file.

• Merge the split elements in the perBase.chr${chrLabel}.answer.bed result back into contiguous ranges, using the ID prefix to group them. In a column adjacent to the newly re-merged contiguous range, print all the scores into a row vector that is delimited with some useful character.

Doing step 1 would remove the need to use a for loop on the bedmap statement. You then just have one bedmap operation against the master phyloP file, instead of looping through each per-chromosome phyloP file.

Step 2 would let you more easily process contiguous regions. For example, you might be working with a window around a set of transcription start sites and you want to know what the mean conservation is across each window.

Hope this gets you started!