feat: speed improvements for primer pair hit building from single primer hits #99
There are no files selected for viewing
| Original file line number | Diff line number | Diff line change | ||||||
|---|---|---|---|---|---|---|---|---|
|
|
@@ -90,6 +90,7 @@ | |||||||
| from prymer.api.oligo import Oligo | ||||||||
| from prymer.api.primer_pair import PrimerPair | ||||||||
| from prymer.api.span import Span | ||||||||
| from prymer.api.span import Strand | ||||||||
| from prymer.offtarget.bwa import BWA_EXECUTABLE_NAME | ||||||||
| from prymer.offtarget.bwa import BwaAlnInteractive | ||||||||
| from prymer.offtarget.bwa import BwaHit | ||||||||
|
|
@@ -126,6 +127,24 @@ class OffTargetResult: | |||||||
| right_primer_spans: list[Span] = field(default_factory=list) | ||||||||
|
|
||||||||
|
|
||||||||
| @dataclass(init=True) | ||||||||
| class PrimerPairBwaHitsBySideAndStrand: | ||||||||
| """A helper class for storing BWA hits for the left and right primers of a primer pair organized | ||||||||
| by left/right and hit strand. | ||||||||
|
|
||||||||
| Attributes: | ||||||||
| left_positive: A list of BwaHit objects on the positive strand for the left primer. | ||||||||
| left_negative: A list of BwaHit objects on the negative strand for the left primer. | ||||||||
| right_positive: A list of BwaHit objects on the positive strand for the right primer. | ||||||||
| right_negative: A list of BwaHit objects on the negative strand for the right primer. | ||||||||
| """ | ||||||||
|
|
||||||||
| left_positive: list[BwaHit] = field(default_factory=list) | ||||||||
| left_negative: list[BwaHit] = field(default_factory=list) | ||||||||
| right_positive: list[BwaHit] = field(default_factory=list) | ||||||||
| right_negative: list[BwaHit] = field(default_factory=list) | ||||||||
|
|
||||||||
|
|
||||||||
| class OffTargetDetector(AbstractContextManager): | ||||||||
| """A class for detecting off-target mappings of primers and primer pairs that uses a custom | ||||||||
| version of "bwa aln" named "bwa-aln-interactive". | ||||||||
|
|
@@ -334,26 +353,71 @@ def _build_off_target_result( | |||||||
| result: OffTargetResult | ||||||||
|
|
||||||||
| # Get the mappings for the left primer and right primer respectively | ||||||||
| left_bwa_result: BwaResult = hits_by_primer[primer_pair.left_primer.bases] | ||||||||
| right_bwa_result: BwaResult = hits_by_primer[primer_pair.right_primer.bases] | ||||||||
|
|
||||||||
| # If there are too many hits, this primer pair will not pass. Exit early. | ||||||||
| if ( | ||||||||
| left_bwa_result.hit_count > self._max_primer_hits | ||||||||
| or right_bwa_result.hit_count > self._max_primer_hits | ||||||||
| ): | ||||||||
| result = OffTargetResult(primer_pair=primer_pair, passes=False) | ||||||||
| return result | ||||||||
|
|
||||||||
| # Get the set of reference names with hits | ||||||||
| hits_by_refname: dict[str, PrimerPairBwaHitsBySideAndStrand] = { | ||||||||
| hit.refname: PrimerPairBwaHitsBySideAndStrand() | ||||||||
| for hit in left_bwa_result.hits + right_bwa_result.hits | ||||||||
|
There was a problem hiding this comment. These feels pretty inefficient .. it will concatenate two fairly large lists and (I think?) generate a I think something like this would work: for refname in set(hit.refname for hit in itertools.chain(left_bwa_result.hits, right_bwa_result.hits) |
||||||||
| } | ||||||||
|
|
||||||||
| # Split the hits for left and right by reference name and strand | ||||||||
| for hit in left_bwa_result.hits: | ||||||||
| if hit.negative: | ||||||||
| hits_by_refname[hit.refname].left_negative.append(hit) | ||||||||
| else: | ||||||||
| hits_by_refname[hit.refname].left_positive.append(hit) | ||||||||
|
|
||||||||
| for hit in right_bwa_result.hits: | ||||||||
| if hit.negative: | ||||||||
| hits_by_refname[hit.refname].right_negative.append(hit) | ||||||||
| else: | ||||||||
| hits_by_refname[hit.refname].right_positive.append(hit) | ||||||||
|
There was a problem hiding this comment. Now that I see this ... I think actually maybe you don't care any more about left vs. right here and beyond this point? Sorry, as I think I may have sent you down this path. Specifically if a single left primer happens to bind to the genome in two places, close together, in the right relative orientations, that could still cause problems. So I think it's really just one collection of + strand mappings and one collection of - strand mappings. |
||||||||
|
|
||||||||
| # Build amplicons from hits on the same reference with valid relative orientation | ||||||||
| amplicons: list[Span] = [] | ||||||||
| for hits in hits_by_refname.values(): | ||||||||
| amplicons.extend( | ||||||||
| self._to_amplicons( | ||||||||
| positive_hits=hits.left_positive, | ||||||||
| negative_hits=hits.right_negative, | ||||||||
| max_len=self._max_amplicon_size, | ||||||||
| strand=Strand.POSITIVE, | ||||||||
|
There was a problem hiding this comment. I'm not sure how to think about strand here. I think it's probably mostly irrelevant and could just always be set to positive? There was a problem hiding this comment. In the When a hit to a primer pair is in the same orientation as the primer pair definition, e.g. left positive and right negative, then it makes sense for the amplicon to also be on the positive strand. If the hit is in the opposite orientation to the primer pair definition, the amplicon should be on the negative strand. I think it makes sense to assign strandedness like this:
There was a problem hiding this comment. As discussed elsewhere, will leave questions of strand and L/L R/R hits to a separate PR. |
||||||||
| ) | ||||||||
| ) | ||||||||
| amplicons.extend( | ||||||||
| self._to_amplicons( | ||||||||
| positive_hits=hits.right_positive, | ||||||||
| negative_hits=hits.left_negative, | ||||||||
| max_len=self._max_amplicon_size, | ||||||||
| strand=Strand.NEGATIVE, | ||||||||
| ) | ||||||||
| ) | ||||||||
|
|
||||||||
| result = OffTargetResult( | ||||||||
| primer_pair=primer_pair, | ||||||||
| passes=self._min_primer_pair_hits <= len(amplicons) <= self._max_primer_pair_hits, | ||||||||
| spans=amplicons if self._keep_spans else [], | ||||||||
| left_primer_spans=( | ||||||||
| [self._hit_to_span(h) for h in left_bwa_result.hits] | ||||||||
| if self._keep_primer_spans | ||||||||
| else [] | ||||||||
| ), | ||||||||
| right_primer_spans=( | ||||||||
| [self._hit_to_span(h) for h in right_bwa_result.hits] | ||||||||
| if self._keep_primer_spans | ||||||||
| else [] | ||||||||
| ), | ||||||||
| ) | ||||||||
|
|
||||||||
| if self._cache_results: | ||||||||
| self._primer_pair_cache[primer_pair] = replace(result, cached=True) | ||||||||
|
|
@@ -410,19 +474,50 @@ def mappings_of(self, primers: list[PrimerType]) -> dict[str, BwaResult]: | |||||||
| return hits_by_primer | ||||||||
|
|
||||||||
| @staticmethod | ||||||||
| def _to_amplicons( | ||||||||
| positive_hits: list[BwaHit], negative_hits: list[BwaHit], max_len: int, strand: Strand | ||||||||
| ) -> list[Span]: | ||||||||
| """Takes a set of hits for one or more left primers and right primers and constructs | ||||||||
| amplicon mappings anywhere a left primer hit and a right primer hit align in F/R | ||||||||
| orientation up to `maxLen` apart on the same reference. Primers may not overlap. | ||||||||
|
|
||||||||
| Args: | ||||||||
| positive_hits: List of hits on the positive strand for one of the primers in the pair. | ||||||||
| negative_hits: List of hits on the negative strand for the other primer in the pair. | ||||||||
| max_len: Maximum length of amplicons to consider. | ||||||||
| strand: The strand of the amplicon to generate. Set to Strand.POSITIVE if | ||||||||
| `positive_hits` are for the left primer and `negative_hits` are for the right | ||||||||
| primer. Set to Strand.NEGATIVE if `positive_hits` are for the right primer and | ||||||||
| `negative_hits` are for the left primer. | ||||||||
|
|
||||||||
| Raises: | ||||||||
| ValueError: If any of the positive hits are not on the positive strand, or any of the | ||||||||
| negative hits are not on the negative strand. If hits are present on more than one | ||||||||
| reference. | ||||||||
| """ | ||||||||
| if any(h.negative for h in positive_hits): | ||||||||
| raise ValueError("Positive hits must be on the positive strand.") | ||||||||
| if any(not h.negative for h in negative_hits): | ||||||||
| raise ValueError("Negative hits must be on the negative strand.") | ||||||||
|
|
||||||||
| refnames: set[str] = {h.refname for h in positive_hits + negative_hits} | ||||||||
|
There was a problem hiding this comment.
Suggested change
|
||||||||
| if len(refnames) > 1: | ||||||||
| raise ValueError("Hits are present on more than one reference.") | ||||||||
|
|
||||||||
|
|
||||||||
| amplicons: list[Span] = [] | ||||||||
| for positive_hit, negative_hit in itertools.product(positive_hits, negative_hits): | ||||||||
|
There was a problem hiding this comment. I think this is good enough if it meets your performance needs @ameynert. As noted elsewhere, if you had many many hits to one reference, you could speed this up by sorting the positive and negative strand hits, looping over list indices, and avoiding whole swathes of pairs. There was a problem hiding this comment. I've tested this separately and the speed improvement by sorting the pos/neg strand hits is significant when there are a large number of hits returned, e.g. in a scenario where the max mismatches parameters are relaxed, one of the primer sequences is getting 100+ hits, and there are multiple highly similar reference contigs in the indexed genome that's being searched. |
||||||||
| if ( | ||||||||
| negative_hit.start > positive_hit.end | ||||||||
| and (negative_hit.end - positive_hit.start + 1) <= max_len | ||||||||
|
There was a problem hiding this comment. Should we allow for situations where the two primers' mappings overlap on the genome? E.g. maybe the check should be: negative_hit.start > positive_hit.start and negative_hit.end > positive_hit.endThere was a problem hiding this comment. That seems reasonable. Do you think we'd ever see a situation where a primer is a palindrome and hits the same target both forward and reverse? There was a problem hiding this comment. This change causes a number of tests to fail. Since it's a breaking change, I'd like to address it in a separate PR. |
||||||||
| ): | ||||||||
| amplicons.append( | ||||||||
| Span( | ||||||||
| refname=positive_hit.refname, | ||||||||
| start=positive_hit.start, | ||||||||
| end=negative_hit.end, | ||||||||
| strand=strand, | ||||||||
| ) | ||||||||
| ) | ||||||||
|
|
||||||||
| return amplicons | ||||||||
|
|
||||||||
|
|
||||||||
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Could we just return here? Is there any value to caching the result in this case, where it's so little work to compute the result?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Caching is required by the unit test at
prymer/tests/offtarget/test_offtarget.py
Line 99 in 2656adf
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Ah, ok. Perhaps I'll "fix" this in a separate PR.