calling SNPs on PacBio fltnc_mapped.bam reads

Question

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

Answer 1

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