Answers.Body: <p>There are definitely better options than a simple mapped consensus, especially if you're using a low[er] accuracy / high-speed mapper like minimap2.</p>
  <p>Have a look at Clair3; that seems to be fairly popular at the moment for long reads:</p>
  <p><a href="https://github.com/HKU-BAL/Clair3" rel="nofollow noreferrer">https://github.com/HKU-BAL/Clair3</a></p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2024-11-20T03:18:48.727Z
Answers.LastActivityDate: 2024-11-20T03:18:48.727Z
Answers.UserDisplayName: gringer
Body: <p>I am trying to iterate over PacBio fltnc_mapped.bam reads for SNPs. I am getting 92% filtering success with T-base SNP iteration and 74% filtering success with C-base SNP at the same position. All the reads  have MAPQ 60. Is there a way of stringently calling the SNPs or this is as good as it gets? My script:</p>
  <pre class="lang-sh prettyprint-override"><code># Load module
  module load samtools
  
  # Extract reads at SNP position
  samtools view -h fltnc_mapped.bam chr1:138039992-138039992 &gt; temp_reads.sam
  
  # Separate the header and the reads for easier manipulation
  grep &quot;^@&quot; temp_reads2.sam &gt; header2.sam
  grep -v &quot;^@&quot; temp_reads2.sam &gt; reads2.sam
  
  # Filter reads based on the SNPs
  awk -v pos=138039992 -v ref_allele=T -v alt_allele=C '
    BEGIN { FS=&quot;\t&quot;; OFS=&quot;\t&quot; }
    /^@/ { print <span class="math-container">$0; next }  # Print header lines
    {
    # Calculate position in read corresponding to SNP
    for (i = 1; i &lt;= length($10); i++) {
      read_pos = $4 + i - 1
      if (read_pos == pos) {
      base = substr($10, i, 1)
      if (base == ref_allele) {
      print $0 &gt; "Ptprcb_reads2.sam"
      } else if (base == alt_allele) {
        print $</span>0 &gt; &quot;Ptprca_reads2.sam&quot;
      }
     }
    }
   }' temp_reads2.sam
  
  # Add the header back to the filtered reads SAM files
  cat header2.sam Ptprca_reads2.sam &gt; Ptprca_reads_with_header2.sam
  cat header2.sam Ptprcb_reads2.sam &gt; Ptprcb_reads_with_header2.sam
  
  # Convert SAM files to BAM files
  samtools view -bS Ptprca_reads_with_header2.sam &gt; Ptprca_reads2.bam
  samtools view -bS Ptprcb_reads_with_header2.sam &gt; Ptprcb_reads2.bam
  
  # Clean up intermediate files
  rm reads2.sam header2.sam Ptprca_reads2.sam Ptprcb_reads2.sam 
    Ptprca_reads_with_header2.sam Ptprcb_reads_with_header2.sam
  </code></pre>
  
Tags: <snp><pacbio>
Title: calling SNPs on PacBio fltnc_mapped.bam reads
ThreadText: # calling SNPs on PacBio fltnc_mapped.bam reads
  
  (Asked by: username_0 on 2024-11-19)
  
  I am trying to iterate over PacBio fltnc\_mapped.bam reads for SNPs. I am getting 92% filtering success with T-base SNP iteration and 74% filtering success with C-base SNP at the same position. All the reads have MAPQ 60. Is there a way of stringently calling the SNPs or this is as good as it gets? My script:
  
  ```sh
  # Load module
  module load samtools
  
  # Extract reads at SNP position
  samtools view -h fltnc_mapped.bam chr1:138039992-138039992 > temp_reads.sam
  
  # Separate the header and the reads for easier manipulation
  grep "^@" temp_reads2.sam > header2.sam
  grep -v "^@" temp_reads2.sam > reads2.sam
  
  # Filter reads based on the SNPs
  awk -v pos=138039992 -v ref_allele=T -v alt_allele=C '
    BEGIN { FS="\t"; OFS="\t" }
    /^@/ { print $0; next }  # Print header lines
    {
    # Calculate position in read corresponding to SNP
    for (i = 1; i <= length($10); i++) {
      read_pos = $4 + i - 1
      if (read_pos == pos) {
      base = substr($10, i, 1)
      if (base == ref_allele) {
      print $0 > "Ptprcb_reads2.sam"
      } else if (base == alt_allele) {
        print $0 > "Ptprca_reads2.sam"
      }
     }
    }
   }' temp_reads2.sam
  
  # Add the header back to the filtered reads SAM files
  cat header2.sam Ptprca_reads2.sam > Ptprca_reads_with_header2.sam
  cat header2.sam Ptprcb_reads2.sam > Ptprcb_reads_with_header2.sam
  
  # Convert SAM files to BAM files
  samtools view -bS Ptprca_reads_with_header2.sam > Ptprca_reads2.bam
  samtools view -bS Ptprcb_reads_with_header2.sam > Ptprcb_reads2.bam
  
  # Clean up intermediate files
  rm reads2.sam header2.sam Ptprca_reads2.sam Ptprcb_reads2.sam 
    Ptprca_reads_with_header2.sam Ptprcb_reads_with_header2.sam
  
  ```
  
  **Tags:** snp, pacbio
  
  
  ---
  
  ## 1 answers
  
  ### Answer by username_1 on 2024-11-20. Score: 0
  
  There are definitely better options than a simple mapped consensus, especially if you're using a low[er] accuracy / high-speed mapper like minimap2.
  
  Have a look at Clair3; that seems to be fairly popular at the moment for long reads:
  
  <https://github.com/HKU-BAL/Clair3>
  
  
  ---
  
  