Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2019-07-23T08:22:56.100Z
Answers.Text: 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?
Answers.UserDisplayName: terdon
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2019-07-24T14:28:06.060Z
Answers.Text: 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
Answers.UserDisplayName: Bioathlete
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2019-07-23T03:13:50.403Z
Answers.LastActivityDate: 2019-07-23T03:13:50.403Z
Answers.UserDisplayName: Bioathlete
Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2019-07-22T19:56:45.643Z
Answers.Text: 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).
Answers.UserDisplayName: terdon
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2019-07-22T20:13:31.520Z
Answers.Text: 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.
Answers.UserDisplayName: swbarnes2
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2019-07-23T05:50:53.000Z
Answers.Text: @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 ./.
Answers.UserDisplayName: Wouter De Coster
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2019-07-23T08:09:21.453Z
Answers.Text: @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.
Answers.UserDisplayName: terdon
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2019-07-23T15:47:15.983Z
Answers.Text: Wouldn't it be easier to filter the vcf?
Answers.UserDisplayName: swbarnes2
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2019-07-24T07:33:02.230Z
Answers.Text: @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.
Answers.UserDisplayName: terdon
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2019-07-24T16:02:26.147Z
Answers.Text: 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.
Answers.UserDisplayName: swbarnes2
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2019-07-22T19:45:36.713Z
Answers.LastActivityDate: 2019-07-22T19:45:36.713Z
Answers.UserDisplayName: swbarnes2
Body: <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>
  
Tags: <sam><primer><base-clipping>
Title: Hard-clip primers in a bam file
ThreadText: # 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**: @username_2 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**: @username_2 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**: @username_2 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.  
  
  ---
  
  