I have a bed file and wanted to calculated average PhastCon score over the regions, is there a tool which could do so? I tried using cmotifs but it submits jobs and at times does not produce timely results.

You could use <a href="http://bedops.readthedocs.org/en/latest/content/reference/statistics/bedmap.html#score-operations">*bedmap --mean*</a> here, to map PhastCon signal onto your <a href="http://bedops.readthedocs.org/en/latest/content/reference/file-management/sorting/sort-bed.html">sorted</a> regions.

First, grab PhastCon data from UCSC and convert WIG to a compressed form of BED called <a href="http://bedops.readthedocs.org/en/latest/content/reference/file-management/compression/starch.html">Starch</a>, using a tool called <a href="http://bedops.readthedocs.org/en/latest/content/reference/file-management/conversion/wig2bed.html">*wig2starch*</a>:

```
$ rootPath="http://hgdownload.cse.ucsc.edu/goldenpath/hg19/phastCons100way/hg19.100way.phastCons"
$ for i in `seq 1 22` X Y; \

    do \
        echo "getting PhastCon data for chromosome chr${i}..."; \
        wigFn="chr${i}.phastCons100way.wigFix"; \
        url="${rootPath}/${wigFn}.gz"; \
        wget -qO- ${url} | gunzip -c - | wig2starch - > ${wigFn}.starch; \
    done
```

We could have used `wig2bed` instead of `wig2starch`, but note that the final set of PhastCons files could take between 10-20 GB compressed, and uncompressed BED files could require about 10x more disk space. It takes a little longer to make Starch files, but it will take less disk space. I'd recommend some kind of compression, if you plan to do repeated map calculations.

Second, map PhastCon signal files to your <a href="http://bedops.readthedocs.org/en/latest/content/reference/file-management/sorting/sort-bed.html">sorted</a> regions file by chromosome:

```
$ for i in `seq 1 22` X Y; \
    do \
        echo "mapping signal for chromosome chr${i}..."; \
        phastConFn="chr${i}.phastCons100way.wigFix.starch"; \
        bedmap --echo --mean --chrom chr${i} regions.bed ${phastConFn} > regions_with_avg_phastcons.${i}.bed;  \
    done
```

Third (optionally), union the per-chromosome results to a single file, using <a href="http://bedops.readthedocs.org/en/latest/content/reference/set-operations/bedops.html#everything-u-everything">*bedops --everything*</a>:

    $ bedops --everything regions_with_avg_phastcons.*.bed > regions_with_avg_phastcons.bed

Doing the *bedmap* step with `--chrom` allows parallelization with grid job schedulers or GNU Parallel, splitting the tasks by chromosome. This could reduce calculation time down to that required for the largest chromosome of regions.

One thing I would note is that the average value calculated is over scores for mapped elements. If you have gaps of scores along your mapped region, those gaps would not contribute to the calculation of the mean signal.

I do not think you can interpolate PhastCon conservation signal between a gap, but it has been a while since I have worked with this particular dataset - its authors could probably confirm or negate this.

It might generally be useful to use *bedmap* with the <a href="http://bedops.readthedocs.org/en/latest/content/reference/statistics/bedmap.html#score-operations">*--skip-unmapped*</a> option to filter out results for unmapped regions.