<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>

CC BY-SA 4.0
2024-11-20T03:18:48.727Z
2024-11-20T03:18:48.727Z
gringer
<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>

<snp><pacbio>
calling SNPs on PacBio fltnc_mapped.bam reads
# 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:
# 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>


---

