I have a bam file and a bed file (with a thousand regions).
I would like, for each region in the bed, to get all the reads that mapped to this region (relatively straightforward) but I also want to get, if any, other regions mapped by these same reads.
I could not find something useful in this regard in the sam documentation.
bedtools bamtobed also seems promising but I think it output only the unique alignments (from the documentation: "There is not yet a standard method for reporting multiple alignments using BAM.").
I mapped Illumina reads with bwa, so I have access to the XA tags. But it will take a long time to parse the bam to get the XA tag for each region of the bed.
Does anyone tried to do something similar ?