I'm most familiar with R/Bioconductor, and the general workflow I would use would be as follows (using the ShortRead, Rsamtools, rtracklayer and IRanges packages)

Firstly, run import() on a bed file of your SureSelect regions - this allows you to use them to find overlaps with your actual data later.

Use scanBam() to read in a BAM file of aligned reads, convert to a GRanges (or GRangesList for multiple samples) using as() and then call coverage() on each GRanges object - this gives you "times coverage" per basepair of the genome in a relatively compact run-length encoded format.

Edit: alternatively once you have GRanges objects you could just call countOverlaps() to count how many reads are in each of your SureSelect regions if you don't need basepair resolution information).

You can then look for regions with less than say 20x coverage using slice(lower=0, upper=20) and then find overlaps of these regions of low coverage with your SureSelect targets using findOverlaps().

Alternatively you could subset the genome wide coverage objects to the SureSelect target regions using Views() and then perform per-region coverage summarisations using viewMins(), viewsMaxs(), viewMeans() etc.

All up this shouldn't take more than a few hours start to finish for a first run through - Good luck!