<p>I'm trying to find all the reads (by name) from a BAM file that align to various regions in a bed file. Right now I can do this with <code>bedtools</code> using <code>intersectBed</code>:</p>

<pre><code>intersectBed -abam reads.bam -wo -f 1 -b regions.bed -bed
</code></pre>

<p>From this one can parse all the read ids that land in every interval in <code>regions.bed</code>, but it's not very compact. Is there a way to get <code>bedtools</code> to natively transform this into a more compact format, e.g.</p>

<pre><code>chr1 x y .... read_id1,read_id2,read_id3
</code></pre>

<p>where <code>chr1 x y</code> is a given interval in <code>regions.bed</code> and the comma separated <code>read_id1,...</code> is the list of read ids from <code>reads.bam</code> that fall in that interval. In this compact format, the output BED file would have at most as many entries as there are regions in <code>regions.bed</code>, whereas with the <code>-wo</code> option it can be even larger than the number of reads in <code>reads.bam</code>. Thanks.</p>

The BEDOPS <a href="http://code.google.com/p/bedops/wiki/bedmap">`bedmap` tool</a> with the `--echo`, `--echo-map-id`, `--delim` and `--multidelim` options can also answer this very quickly. You can pipe in reads into `bedmap` by using the <a href="http://code.google.com/p/bedops/wiki/bam2bed">`bam2bed` script</a>:

    $ bam2bed.csh < reads.bam \
        | bedmap --echo --echo-map-id --delim '\t' --multidelim ',' sorted-regions.bed - \
        > answer.bed

The `--echo` option returns the `region` element, while `--echo-map-id` returns a list of `read` elements which overlap the region by one or more bases.

The `answer.bed` file will look like:

    [ region-A in BED format ] \t [ comma-delimited list of read IDs overlapping region-A ]
    [ region-B in BED format ] \t [ comma-delimited list of read IDs overlapping region-B ] 
    ...

As the answer is BED-formatted output, you can very easily pipe this to downstream processes that consume BED data.

Just be sure to sort regions, if not sorted, e.g.:

    $ sort-bed regions.bed > sorted-regions.bed

And, likewise, for reads:

    $ bam2bed.csh < reads.bam | sort-bed - > sorted-reads.bed

Or to put it together with the overall analysis:

    $ bam2bed.csh < reads.bam \
        | sort-bed - \
        | bedmap --echo --echo-map-id --delim '\t' --multidelim ',' sorted-regions.bed - \
        > answer.bed