Precisely clipping bam file to bed coordinates

Question

I have a bed file of coordinates and would like to clip anything outside these regions in my aligned bam files. However, trying programs such as bedtools intersect and samtools ampliconclip (on complementary bed file), they always let some reads slip through outside the ranges. I tried using hard-clip and both-ends options for ampliconclip aswell.

Example (top is clipped bam, middle is unclipped bam and bottom is bed file with regions I would like to keep):

Although most reads outside the designed amplicons are removed, there are still reads outside the regions in the Designed_pool1 file. The problem is that the bed file contains amplicons and outside the amplicons are primers where I do not want to call variants. I cannot get the exact primer sequences because I don't know their length and the manufacturer does not state them, so I cannot use programs that require primer sequences.

Does anyone have some experience with this and know why the programs behave this way?

EDIT:

The exact commands I have used are:

bedtools intersect -a $BamFile -b $BedFile > $BamOut

and

samtools ampliconclip --hard-clip --both-ends -b $ComplBedFile $BamFile > $BamOut

as well as without --hard-clip and --both-ends options. All methods let some sequences stay out of the regions in bed files.

I am not completely familiar with the bedpe format, but I am not using paired end data. I cannot get primer start and stop positions, only amplicon positions.

Answer 1

I am curious whether your pileup visualization tool is respecting the clipping. It is possible that there is a subtlety of samtools ampliconclip that your visualization is not understanding, with regard to how it treats e.g. supplementary alignments or something.

Perhaps more importantly, I notice that the reads that are not being properly clipped in your visualization tend to have one end anchored in each of your BED regions.

Looking at the samtools ampliconclip docs, it seems to be assuming that reads will only be clipped from the ends. They will not be clipped if the region outside the BED is in the middle, if I read it correctly. I think that you may need to include an additional step in which you either trim reads (so that they are not long enough to bridge your regions) or you somehow postprocess your alignments to explicitly clip in the middle.

If true, this seems like a significant enough oversight of samtools that it may be worth logging an issue with the samtools maintenance team.

Also, not clear if it's relevant but reading the samtools ampliconclip docs:

adjustments to the left most mapping position (POS) will mean that coordinate sorted files will need resorting. In such cases the sorting order in the header is set to unknown. Clipping of reads results in template length (TLEN) being incorrect. This can be corrected by samtools fixmates. Any MD and NM aux tags will also be incorrect*, which can be fixed by samtools calmd.

Answer 2

For single-end reads, generically:

$ bam2bed < reads.bam > reads.bed
$ bedmap --echo --fraction-ref 1 reads.bed roi.bed > answer.bed

The BED file answer.bed contains filtered reads which entirely overlap intervals in roi.bed. It contains fields required for backconversion to SAM/BAM.