blasr read length filter?

[ewilbanks]

Hey folks,

I'm working with aligning some pacbio reads to a reference and have a question about blasr that I can't seem to figure out from the documentation. What I want is to only map the corrected reads which are longer than 3kb, to try to avoid spurious mismapping errors in repeat/transposase regions.

The flag -minReadLength doesn't behave as I would expect it to though. When I peak at the results by running the following:

head -5000 blasr.sam | awk '{if ($1 ~/@.*/) next; else if (length($10) < 3000) print length($10);}'

I see that loads of alignments that are shorter than 3k which leaves me puzzled about a few things.

• Anyone know what this flag actually does?
• Am I missing something about what that tenth sequence field in the sam file means?
• The default behavior is to not clip the alignments, so why are do I see alignments shorter than the shortest read I'm aligning? Shouldn't it only be reporting full length alignments?

Thanks for your thoughts!
Lizzy

[rhall]

The -minReadLength sets the minimum length of the reads getting aligned, not the minimum alignment length. Blasr finds best local alignment between a read and a reference, it does not align the entire read to the reference.
To align, and only report alignments above a certain length you can use the pbalign wrapper to blasr, https://github.com/PacificBiosciences/pbalign also part of SMRT Analysis. Or, for aligning low error corrected reads Mummer may be a better option.
Richard.