<p>For single-end reads, generically:</p>
<pre><code>$ bam2bed &lt; reads.bam &gt; reads.bed
$ bedmap --echo --fraction-ref 1 reads.bed roi.bed &gt; answer.bed
</code></pre>
<p>The BED file <code>answer.bed</code> contains filtered reads which entirely overlap intervals in <code>roi.bed</code>. It contains fields required for backconversion to SAM/BAM.</p>

CC BY-SA 4.0
2021-07-06T00:00:38.343Z
2021-07-06T00:00:38.343Z
Alex Reynolds
<p>You can convert a SAM stream (mandatorily with a header) into FASTA (multiFASTA file padded with dashes) using <a href="https://github.com/virus-evolution/gofasta" rel="nofollow noreferrer">gofasta</a>:</p>
<p><code>samtools view -h some.bam | gofasta sam toMultiAlign --start 1171 --end 1177 &gt; region.aln</code></p>
<p>Beware <code>gofasta sam toMultiAlign</code> discards <code>INS</code>ertions from the sample relative to the reference, so the FASTA sequences may not be exactly sub-strings of the original data. You can detect DELetions in the sample reads when <code>-</code> are present relative to reference.</p>
<p>That is not an issue with <code>gofasta sam toPairAlign</code> but is terribly slow and it dumps zillions of single-entry FASTA files into your current working directory. Probably would be worth to add some flags to <code>samtools</code> to omit unmapped reads from the input to feed <code>gofasta</code> with less input.</p>
<p><code>samtools view -h some.bam | gofasta sam toPairAlign --reference someref.fna --start 1171 --end 1177 &gt; region.aln</code></p>

CC BY-SA 4.0
2024-11-05T15:32:11.333Z
2024-11-05T16:07:02.870Z
2024-11-05T16:07:02.870Z
Martin Mokrejs
<p>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.</p>
<p>Perhaps more importantly, <strong>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.</strong></p>
<p>Looking at the samtools ampliconclip docs, it seems to be assuming that <strong>reads will only be clipped from the ends</strong>. 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.</p>
<p>If true, this seems like a significant enough oversight of samtools that it may be worth <a href="https://github.com/samtools/samtools/issues" rel="nofollow noreferrer">logging an issue with the samtools maintenance team</a>.</p>
<p>Also, not clear if it's relevant but reading the samtools ampliconclip docs:</p>
<blockquote>
<p>adjustments to the left most mapping position (POS) will mean that <strong>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.</strong> This can be corrected by samtools fixmates. Any <em>MD and NM aux tags will also be incorrect</em>*, which can be fixed by samtools calmd.</p>
</blockquote>

CC BY-SA 4.0
2021-07-04T17:57:40.373Z
2021-07-04T17:57:40.373Z
Maximilian Press
<p>I wrote <a href="https://jvarkit.readthedocs.io/en/latest/PcrClipReads/" rel="nofollow noreferrer">https://jvarkit.readthedocs.io/en/latest/PcrClipReads/</a></p>
<pre><code>java -jar dist/jvarkit.jar pcrclipreads -B test.bed  samtools-0.1.19/examples/ex1.bam  
</code></pre>

CC BY-SA 4.0
2024-11-05T16:29:41.350Z
2024-11-05T16:29:41.350Z
Pierre
<p>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.</p>
<p>Example (top is clipped bam, middle is unclipped bam and bottom is bed file with regions I would like to keep):</p>
<p><a href="https://i.sstatic.net/VUC23.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/VUC23.png" alt="enter image description here" /></a></p>
<p>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.</p>
<p>Does anyone have some experience with this and know why the programs behave this way?</p>
<p>EDIT:</p>
<p>The exact commands I have used are:</p>
<pre><code>bedtools intersect -a $BamFile -b $BedFile &gt; $BamOut
</code></pre>
<p>and</p>
<pre><code>samtools ampliconclip --hard-clip --both-ends -b $ComplBedFile $BamFile &gt; $BamOut
</code></pre>
<p>as well as without --hard-clip and --both-ends options. All methods let some sequences stay out of the regions in bed files.</p>
<p>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.</p>

CC BY-SA 4.0
2021-06-18T00:51:39.707Z
Maybe you need to sort the bam file before trying to clip.  Did you try that already?
user3479780
CC BY-SA 4.0
2021-06-18T13:13:24.013Z
Yes, the bam files are sorted
KasperH
<base-clipping>
Precisely clipping bam file to bed coordinates
# Precisely clipping bam file to bed coordinates

(Asked by: username_0 on 2021-06-10)

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.

**Tags:** base-clipping
> **Comments:**
> **username_1**: Maybe you need to sort the bam file before trying to clip.  Did you try that already?  
> **username_0**: Yes, the bam files are sorted  

---

## 4 answers

### Answer by username_2 on 2021-07-06. Score: 1

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.


---

### Answer by username_3 on 2024-11-05. Score: 1

You can convert a SAM stream (mandatorily with a header) into FASTA (multiFASTA file padded with dashes) using [gofasta](https://github.com/virus-evolution/gofasta):

`samtools view -h some.bam | gofasta sam toMultiAlign --start 1171 --end 1177 > region.aln`

Beware `gofasta sam toMultiAlign` discards `INS`ertions from the sample relative to the reference, so the FASTA sequences may not be exactly sub-strings of the original data. You can detect DELetions in the sample reads when `-` are present relative to reference.

That is not an issue with `gofasta sam toPairAlign` but is terribly slow and it dumps zillions of single-entry FASTA files into your current working directory. Probably would be worth to add some flags to `samtools` to omit unmapped reads from the input to feed `gofasta` with less input.

`samtools view -h some.bam | gofasta sam toPairAlign --reference someref.fna --start 1171 --end 1177 > region.aln`


---

### Answer by username_4 on 2021-07-04. Score: 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](https://github.com/samtools/samtools/issues).

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 by username_5 on 2024-11-05. Score: 0

I wrote https://jvarkit.readthedocs.io/en/latest/PcrClipReads/
java -jar dist/jvarkit.jar pcrclipreads -B test.bed  samtools-0.1.19/examples/ex1.bam


---

