The other way to do this is to use hgWiggle program from UCSC source code. You can download the program from here, and you will need to prepare a three line ".hg.conf" file in your home directory.
.hg.conf should contain the exact following lines:
You should download the related .wib file from UCSC. For human hg18 assembly they are located at ftp://hgdownload.cse.ucsc.edu/gbdb/hg18/ . Then do "chmod 600 .hg.conf " on the file. Now we are ready to use hgWiggle. See this wiki page for all the other details.
As an example, you can get the statistics for chr16 for phastCons scores as shown below.
hgWiggle -db=hg18 -chr=chr16 -doStats phastCons44wayPlacental
you can also give a bed file with genomic coordinates using -bedFile option, then you will get phastCons summary for the regions in the bed file.
It is also pretty straight forward to call the hgWiggle from R using system() function, so you can include this functionality in your R code. Although I intended to do that here when I started the post, I thought that it is trivial to use the system() command and call any command line tool, so I won't be doing that for now. I might add it later on when I have more time.