<p>Another option would be to adjust the quality of the primer bases to be below the threshold that your variant caller uses.  In my case I have "masked" the primer bases by dropping the quality to 0.  This prevents the variants from being called based that are synthetic but if you have an overlapping amplicon it will call bases from those reads.</p>

CC BY-SA 4.0
2019-07-23T08:22:56.100Z
That sounds promising. How would I go about that given a bedpe file of primer locations and a bam file? Is there a tool that can take those two files as input and set the MAPQ to 0? Or do you mean changing the quality score in the original fastq file?
terdon
CC BY-SA 4.0
2019-07-24T14:28:06.060Z
The tools we use are custom built but in general we detect the primer for each amplicon and then edit the synthetic (primer) baseQs not the MapQ as 0 for those reads.  You could do it either in the FASTQs prior to alignment or on the BAMs before variant calling.  If you don't have overlapping amplicons this would probably be overkill and filtering the result VCF or redefine the regions of interested to the variant caller might be easier
Bioathlete
CC BY-SA 4.0
2019-07-23T03:13:50.403Z
2019-07-23T03:13:50.403Z
Bioathlete
<blockquote>
  <p>I know that the primer sequences are part of the genome and a variant
  there is no less valid than a variant anywhere else,</p>
</blockquote>

<p>You won't see variants under PCR primers.  You get the primer sequence itself in the vast majority of reads.</p>

<p>I wouldn't clip the bam file, I'd clip either the fastq, or clip within bcl2fastq while making the fastqs.</p>

CC BY-SA 4.0
2019-07-22T19:56:45.643Z
Thanks, but I do actually get variants in the primer. That's the problem :( And I really don't want to clip the fastq. Both because the primers can help the alignment process and because that's harder to do because I would require fasta sequences of the primers and not regions which is what I have and what assay manufacturers provide. Finally, I never have the bcl files, I will always start from FASTQ (people send me their fastq files).
terdon
CC BY-SA 4.0
2019-07-22T20:13:31.520Z
Variants under the primers can't be real.  Just look at how PCR works.  You will only see a variant under the forward primer if amplification happens from the reverse primer sitting on the orignal template. Once you get amplification happening off of templates made by extending off of the forward primer, the variant is gone.  The orignal template will be outnumbered a billion to 1.
swbarnes2
CC BY-SA 4.0
2019-07-23T05:50:53.000Z
@swbarnes2 primers may be inaccurately synthetized and lead to false positives, theoretically, but I would be more worried about false negatives. A true variant under the primer will be lost, for which you should not call a 0/0 but a ./.
Wouter De Coster
CC BY-SA 4.0
2019-07-23T08:09:21.453Z
@swbarnes2 that's a very valid point. But my main issue is that I have a client who doesn't want to see any variants in the primers, ever. And I am calling variants in th primers, for whatever reason. It's possible that the primer itself has a mismatch, I guess, and that's what I'm seeing. But I need to completely hide these regions from the variant caller.
terdon
CC BY-SA 4.0
2019-07-23T15:47:15.983Z
Wouldn't it be easier to filter the vcf?
swbarnes2
CC BY-SA 4.0
2019-07-24T07:33:02.230Z
@swbarnes2 indeed it would and that's probably what I'll end up doing. I was hoping there'd be a tool to do this at the bam level to avoid calling them in the first place, but since it looks like there isn't, I'll probably either filter the VCF or extract the target regions from my bedpe file and direct the caller to only call there.
terdon
CC BY-SA 4.0
2019-07-24T16:02:26.147Z
Trimming at the bam stage would altering at multiple fields in the bam; the sequence itself, the quality string, the CIGAR string, the sequence length, plus optional tags at the end, it's a massive headache.
swbarnes2
CC BY-SA 4.0
2019-07-22T19:45:36.713Z
2019-07-22T19:45:36.713Z
swbarnes2
<p>I am working on a sample which was sequenced using an amplicon kit and I would like to ensure that no variants are called in the PCR primer sequences. I know that the primer sequences are part of the genome and a variant there is no less valid than a variant anywhere else, but for reasons I cannot control, I need to ensure that no variants can be called in the primer regions. </p>

<p>My first thought was to use a soft-clipping tool like <a href="https://github.com/tommyau/bamclipper" rel="nofollow noreferrer">BAMClipper</a>. However, the variant caller I am using (<a href="https://www.sentieon.com/products/" rel="nofollow noreferrer">Sentieon's implementation of GATK</a>) preforms a local realignment when calling and if it finds that soft-clipped regions align perfectly against the genome, it will remove the soft clips. This makes perfect sense for cases where the softclips were added by the aligner, but sadly breaks my workflow since I need the softclipped primers to remain clipped.</p>

<p>I would therefore like to hard clip, to completely remove the primer sequences from the bam file before variant calling. I have found one tool that seems to be able to do this, <a href="http://fulcrumgenomics.github.io/fgbio/tools/latest/TrimPrimers.html" rel="nofollow noreferrer">TrimPrimers</a> by Fulcrum Genomics, but this fails on my input bam file with <a href="https://pastebin.com/QRYK2kJY" rel="nofollow noreferrer">an error</a> I can't really make much of. </p>

<p>Is there any other tool that can take a list of regions (e.g. a bedpe file) and a BAM and hard-clip the regions from the BAM?</p>

<sam><primer><base-clipping>
Hard-clip primers in a bam file
# Hard-clip primers in a bam file

(Asked by: username_0 on 2019-07-22)

I am working on a sample which was sequenced using an amplicon kit and I would like to ensure that no variants are called in the PCR primer sequences. I know that the primer sequences are part of the genome and a variant there is no less valid than a variant anywhere else, but for reasons I cannot control, I need to ensure that no variants can be called in the primer regions.

My first thought was to use a soft-clipping tool like [BAMClipper](https://github.com/tommyau/bamclipper). However, the variant caller I am using ([Sentieon's implementation of GATK](https://www.sentieon.com/products/)) preforms a local realignment when calling and if it finds that soft-clipped regions align perfectly against the genome, it will remove the soft clips. This makes perfect sense for cases where the softclips were added by the aligner, but sadly breaks my workflow since I need the softclipped primers to remain clipped.

I would therefore like to hard clip, to completely remove the primer sequences from the bam file before variant calling. I have found one tool that seems to be able to do this, [TrimPrimers](http://fulcrumgenomics.github.io/fgbio/tools/latest/TrimPrimers.html) by Fulcrum Genomics, but this fails on my input bam file with [an error](https://pastebin.com/QRYK2kJY) I can't really make much of.

Is there any other tool that can take a list of regions (e.g. a bedpe file) and a BAM and hard-clip the regions from the BAM?

**Tags:** sam, primer, base-clipping

---

## 2 answers

### Answer by username_1 on 2019-07-23. Score: 2

Another option would be to adjust the quality of the primer bases to be below the threshold that your variant caller uses. In my case I have "masked" the primer bases by dropping the quality to 0. This prevents the variants from being called based that are synthetic but if you have an overlapping amplicon it will call bases from those reads.

> **Comments:**
> **username_0**: That sounds promising. How would I go about that given a bedpe file of primer locations and a bam file? Is there a tool that can take those two files as input and set the MAPQ to 0? Or do you mean changing the quality score in the original fastq file?  
> **username_1**: The tools we use are custom built but in general we detect the primer for each amplicon and then edit the synthetic (primer) baseQs not the MapQ as 0 for those reads.  You could do it either in the FASTQs prior to alignment or on the BAMs before variant calling.  If you don't have overlapping amplicons this would probably be overkill and filtering the result VCF or redefine the regions of interested to the variant caller might be easier  

---

### Answer by username_2 on 2019-07-22. Score: 1

> I know that the primer sequences are part of the genome and a variant
> there is no less valid than a variant anywhere else,

You won't see variants under PCR primers. You get the primer sequence itself in the vast majority of reads.

I wouldn't clip the bam file, I'd clip either the fastq, or clip within bcl2fastq while making the fastqs.

> **Comments:**
> **username_0**: Thanks, but I do actually get variants in the primer. That's the problem :( And I really don't want to clip the fastq. Both because the primers can help the alignment process and because that's harder to do because I would require fasta sequences of the primers and not regions which is what I have and what assay manufacturers provide. Finally, I never have the bcl files, I will always start from FASTQ (people send me their fastq files).  
> **username_2**: Variants under the primers can't be real.  Just look at how PCR works.  You will only see a variant under the forward primer if amplification happens from the reverse primer sitting on the orignal template. Once you get amplification happening off of templates made by extending off of the forward primer, the variant is gone.  The orignal template will be outnumbered a billion to 1.  
> **username_3**: @swbarnes2 primers may be inaccurately synthetized and lead to false positives, theoretically, but I would be more worried about false negatives. A true variant under the primer will be lost, for which you should not call a 0/0 but a ./.  
> **username_0**: @swbarnes2 that's a very valid point. But my main issue is that I have a client who doesn't want to see any variants in the primers, ever. And I am calling variants in th primers, for whatever reason. It's possible that the primer itself has a mismatch, I guess, and that's what I'm seeing. But I need to completely hide these regions from the variant caller.  
> **username_2**: Wouldn't it be easier to filter the vcf?  
> **username_0**: @swbarnes2 indeed it would and that's probably what I'll end up doing. I was hoping there'd be a tool to do this at the bam level to avoid calling them in the first place, but since it looks like there isn't, I'll probably either filter the VCF or extract the target regions from my bedpe file and direct the caller to only call there.  
> **username_2**: Trimming at the bam stage would altering at multiple fields in the bam; the sequence itself, the quality string, the CIGAR string, the sequence length, plus optional tags at the end, it's a massive headache.  

---

