Hi all,
I'm currently working on a PacBio amplicon sequencing project. The libraries were constructed using M13-barcoded primers, and our goal is to detect indels or large deletions within specific regions of the amplicons.
After demultiplexing, I performed additional trimming to ensure that all reads in the final dataset were full-length amplicons. Specifically, I filtered for reads containing both the first and last 30 bp of the target sequence (allowing up to 5 mismatches), so the expectation is that:
The reads may contain indels or large internal deletions, but both the 5′ and 3′ ends should match the reference and therefore align properly.
However, when I aligned these reads to the reference using minimap2 and visualized the results in IGV, I noticed that many reads were only partially aligned, including in the negative control library, where no deletions are expected. (See figure.)
Here is another example where I am expecting some of the reads to have indels or large deletions:
minimap2 -x map-pb -a --secondary=no --end-bonus=100 --score-N=0 "$ref_file" "$fastq_file" | samtools view -Sb - | samtools sort -o "$output_bam" | samtools index "$output_bam"
To validate this, I also tried aligning the same reads using pbmm2 with the CCS preset:
pbmm2 align /path/to/MYC_PCR.fa yll504.bam /path/to/yll504.aligned.bam --preset CCS --sort --min-length 50 --best-n 1 --log-level DEBUG --num-threads 4
Unfortunately, the result was similar—many reads still showed only partial alignment, even though they should cover the full amplicon ends.
My questions: Could this be caused by an issue during library construction (e.g., degradation, incomplete extension)? Could the alignment settings be suboptimal for detecting reads with large internal deletions while preserving full-end alignment? Are there better aligners or parameters suited for this type of amplicon with variable internal structure? Has anyone successfully aligned PacBio amplicons like this with consistent end-to-end mapping?
Any suggestions or insights would be very helpful—thank you in advance!
•
written
by
louisflower1999 •
1 answer
The soft-clipped bases will tell you which of your three hypotheses is right, and it's a two-minute check rather than a guess. Pull them out and look at what they actually are:
samtools view your.bam | awk '$6 ~ /^[0-9]+S/ {print $1, $6}' | head
then extract a few of those clipped stretches and BLAST them.
If they turn out to be M13 tail plus barcode, then your filtering step selected for full-length reads but never actually trimmed the ends, so there is nothing for the aligner to place and IGV shows the read stopping short of the amplicon. If instead they are sequence from somewhere else in your own amplicon, you have PCR chimeras from template switching, which is common at 4.4 kb and would appear in the negative control exactly as you describe. If they are neither, it's a library problem.
Two parameter notes. If these are HiFi reads then -x map-pb is the wrong preset, that one is for noisy CLR - use map-hifi. And for the large internal deletions you're actually hunting, raise the chaining bandwidth (something like -r 2k,10k), because a deletion longer than the default gap breaks the chain and gets reported as a clipped partial alignment instead of as a deletion.
Log in to answer this question.
-
PacBio HiFi alignment: am I doing this right?
written by matvergara19 •
-
Understanding Re-use of Query Sequences in Long Read Alignments
written by Charles-Alexandre Roy •
-
Processing and calling SVs from PacBio data
written by adarsh •
-
Split large deletion D in CIGAR string and filtering reads
written by marco.barr •
-
different result using minimap2 and pbmm2
written by pingu77 •
-
Amplicon variant calling from PacBio CCS reads
written by vaslanzadeh •
-
Samtools merge Illumina and PB bam file empty
written by talbots •
-
Aligner for large indels (up to 200bp out of 500bp reads)
written by Vitis
-
why does qualimap2 shows different results for number of mismatches for bowtie2 and BWAMEM?
written by maria2019
-
How to shorten runtimes of freebayes on very deep long reads with a smattering of largish indels
written by rleach
Only including this as an alternate option to try (if your reads are HiFi) : https://github.com/PacificBiosciences/pbAA
How long are these amplicons? What is the size of the deletions that you are expecting. What is the length of the two ends that should match (looks like 3'-ends go to the end of the amplicon, if the shots show the full length) so it must be the 5'end that may be getting soft-clipped. IGV will not display soft-clipped reads unless an option is turned on so you may want to try that.
It is certainly possible that something unexpected happened during the experimental phase and that is reflected here.
GenoMax has a good point about trying
pbaa. The tool will output consensus sequences of the different alleles. You can then use MSA to align the consensus sequences and the reference sequence. I suspect that the mapping looks funky because there's an insertion that is causing the mappings to terminate (via clipping).Hi GenoMax,
Thanks so much for your suggestion! Those amplicons are around 4.4kb. And there should be at least 1kb from the left and 800 bp from the right that should match the reference.