One way to do this is to make a pipeline that:

1. Converts the annotations to BED via convert2bed -i gff
2. Uses awk to get one base elements, if the previous stop and current start position values are equal
3. Uses BEDOPS bedops --range to expand (for example) a 10-base window around the resulting one-base elements (adjusting this window to whatever you need)
4. Maps the ranged elements to find overlapping reads via bedmap and convert2bed -i bam

$ convert2bed -i gff < annotations.gff \
    | awk ' \
        BEGIN { \
            previous_chr = "chr0"; \
            previous_stop = -1; \
        } \
        { \
            current_chr = $1; \
            current_start = $2; \
            current_stop = $3; \
            if ((current_start == previous_stop) && (current_chr == previous_chr)) { \
                print current_chr"\t"current_start"\t"(current_start + 1); \
            } \
           previous_chr = current_chr; \
           previous_stop = current_stop; \
        } \
    }' \
    | bedops --everything --range 10 - \
    | bedmap --echo --echo-map --skip-unmapped - <(convert2bed -i bam < reads.bam) \
    > answer.bed

Often the way to solve these problems lies in building one ore more new dataset that when used with bedtools intersect produces what you need.

For example here create file that contains a marker in the form of a 1bp long interval for its start of each feature. Then in another file create markers for the 1bp interval of the end of each feature. Now intersect the two files and you'll end up with an file with intervals that are common. Use this file to intersect with the reads.