Gibson primer design - allow for a sequence not to be amplified

src/opencloning/endpoints/primer_design.py

@@ -4,6 +4,8 @@
 from Bio.Restriction import RestrictionBatch
 from Bio.SeqUtils import gc_fraction
 
+from opencloning.temp_functions import pydna_primer_to_primer_model
+
 from ..dna_functions import get_invalid_enzyme_names
 from opencloning_linkml.datamodel import Primer as PrimerModel
 from opencloning.pydantic_models import PrimerDesignQuery
@@ -24,10 +26,11 @@
 from ..get_router import get_router
 from ..ebic.primer_design import ebic_primers
 from pydantic import BaseModel
+import copy
 
 router = get_router()
 
-PrimerDesignResponse = create_model('PrimerDesignResponse', primers=(list[PrimerModel], ...))
+PrimerDesignResponse = create_model('PrimerDesignResponse', primers=(list[PrimerModel | None], ...))
 
 
 def validate_spacers(spacers: list[str] | None, nb_templates: int, circular: bool):
@@ -101,7 +104,13 @@ async def primer_design_homologous_recombination(
 
 @router.post('/primer_design/gibson_assembly', response_model=PrimerDesignResponse)
 async def primer_design_gibson_assembly(
-    pcr_templates: list[PrimerDesignQuery],
+    pcr_templates: list[PrimerDesignQuery] = Body(
+        ...,
+        description='''The templates to design primers for. If the location is not
+        provided for a pcr_template, primers won't be designed for that part (the part will be included
+        in the Gibson as is).''',
+        min_length=1,
+    ),
     settings: PrimerDesignSettings = Body(description='Primer design settings.', default_factory=PrimerDesignSettings),
     spacers: list[str] | None = Body(
         None,
@@ -116,19 +125,26 @@ async def primer_design_gibson_assembly(
     ),
     circular: bool = Query(False, description='Whether the assembly is circular.'),
 ):
-    """Design primers for Gibson assembly"""
+    """Design primers for Gibson assembly. If the location is not provided for a pcr_template, primers won't be designed for that part."""
+
     # Validate the spacers
     validate_spacers(spacers, len(pcr_templates), circular)
     templates = list()
+    amplify_templates = list()
     for query in pcr_templates:
         dseqr = read_dsrecord_from_json(query.sequence)
-        location = query.location.to_biopython_location()
-        template = location.extract(dseqr)
-        if not query.forward_orientation:
-            template = template.reverse_complement()
-        # For naming the primers
-        template.name = dseqr.name
-        template.id = dseqr.id
+        if query.location is not None:
+            location = query.location.to_biopython_location()
+            template = location.extract(dseqr)
+            if not query.forward_orientation:
+                template = template.reverse_complement()
+            # For naming the primers
+            template.name = dseqr.name
+            template.id = dseqr.id
+            amplify_templates.append(True)
+        else:
+            template = copy.deepcopy(dseqr)
+            amplify_templates.append(False)
         templates.append(template)
     try:
         primers = gibson_assembly_primers(
@@ -139,11 +155,12 @@ async def primer_design_gibson_assembly(
             circular,
             spacers,
             tm_func=lambda x: primer3_calc_tm(x, settings),
+            amplify_templates=amplify_templates,
         )
     except ValueError as e:
         raise HTTPException(400, *e.args)
 
-    return {'primers': primers}
+    return {'primers': [pydna_primer_to_primer_model(primer) if primer is not None else None for primer in primers]}
 
 
 @router.post('/primer_design/simple_pair', response_model=PrimerDesignResponse)
@@ -193,7 +210,7 @@ async def primer_design_simple_pair(
     # This is to my knowledge the only way to get the enzymes
     rb = RestrictionBatch()
     try:
-        fwd, rvs = simple_pair_primers(
+        primers = simple_pair_primers(
             template,
             minimal_hybridization_length,
             target_tm,
@@ -208,7 +225,7 @@ async def primer_design_simple_pair(
     except ValueError as e:
         raise HTTPException(400, *e.args)
 
-    return {'primers': [fwd, rvs]}
+    return {'primers': [pydna_primer_to_primer_model(primer) for primer in primers]}
 
 
 @router.post('/primer_design/ebic', response_model=PrimerDesignResponse)

src/opencloning/primer_design.py

@@ -1,4 +1,5 @@
 from pydna.dseqrecord import Dseqrecord
+from pydna.primer import Primer
 from pydna.design import primer_design, assembly_fragments
 from Bio.SeqFeature import SimpleLocation
 from pydna.utils import locations_overlap, shift_location, location_boundaries
@@ -8,7 +9,6 @@
 from Bio.Data.IUPACData import ambiguous_dna_values as _ambiguous_dna_values
 from typing import Callable
 from .primer3_functions import primer3_calc_tm, PrimerDesignSettings
-from opencloning_linkml.datamodel import Primer as PrimerModel
 
 ambiguous_dna_values = _ambiguous_dna_values.copy()
 # Remove acgt
@@ -92,51 +92,92 @@ def gibson_assembly_primers(
     spacers: list[str] | None = None,
     tm_func: Callable[[str], float] = default_tm_func,
     estimate_function: Callable[[str], float] | None = None,
-) -> list[PrimerModel]:
-
-    initial_amplicons = [
-        primer_design(
-            template,
-            limit=minimal_hybridization_length,
-            target_tm=target_tm,
-            tm_func=tm_func,
-            estimate_function=estimate_function,
-        )
-        for template in templates
-    ]
-
-    for i, amplicon in enumerate(initial_amplicons):
-        if None in amplicon.primers():
+    amplify_templates: list[bool] | None = None,
+) -> list[Primer]:
+
+    if amplify_templates is None:
+        amplify_templates = [True] * len(templates)
+    if len(amplify_templates) != len(templates):
+        raise ValueError('The number of amplify_templates must be the same as the number of templates.')
+    for prev, next in zip(amplify_templates[:-1], amplify_templates[1:]):
+        if not prev and not next:
+            raise ValueError('Two consecutive templates with amplify_templates=False are not allowed.')
+    if circular and (not amplify_templates[-1] and not amplify_templates[0]):
+        raise ValueError('Two consecutive templates with amplify_templates=False are not allowed.')
+
+    if len(templates) == 1 and amplify_templates[0] is False:
+        raise ValueError('amplify_templates cannot be False for a single template.')
+
+    if spacers is not None and not circular:
+        if spacers[0] != '' and not amplify_templates[0]:
+            raise ValueError(
+                'The first spacer must be empty if the first template is not amplified in linear assembly.'
+            )
+        if spacers[-1] != '' and not amplify_templates[-1]:
+            raise ValueError('The last spacer must be empty if the last template is not amplified in linear assembly.')
+
+    # For the function assembly_fragments, maxlink is the maximum length of a Dseqrecord to be considered a spacer.
+    # It's important to check that the amplify_templates=False parts are not longer than maxlink, otherwise, they
+    # would be considered as spacers. This is perhaps not ideal, as there could be a case, for now we just do it
+    # like this.
+
+    maxlink = minimal_hybridization_length * 2
+    if spacers is not None:
+        maxlink = max(len(spacer) for spacer in spacers)
+
+    inputs: list[Amplicon | Dseqrecord] = list()
+    for i, template in enumerate(templates):
+        if amplify_templates[i]:
+            inputs.append(
+                primer_design(
+                    template,
+                    limit=minimal_hybridization_length,
+                    target_tm=target_tm,
+                    tm_func=tm_func,
+                    estimate_function=estimate_function,
+                )
+            )
+        else:
+            if len(template) < maxlink:
+                raise ValueError(
+                    f'Template {template.name} ({len(template)} bps) is shorter than the longest spacer or 2x the minimal hybridization length.'
+                )
+            inputs.append(template)
+
+    for i, amplicon in enumerate(inputs):
+        if amplify_templates[i] is True and None in amplicon.primers():
             raise ValueError(f'Primers could not be designed for template {templates[i].name}, try changing settings.')
 
-    maxlink = 40
     if spacers is not None:
-        maxlink = max(len(spacer) for spacer in spacers)
         spacers = [Dseqrecord(spacer) for spacer in spacers]
-        initial_amplicons_with_spacers = []
+        inputs_withspacers = []
         # For linear assemblies, the first spacer is the first thing
         if not circular:
-            initial_amplicons_with_spacers.append(spacers.pop(0))
-        for amplicon in initial_amplicons:
-            initial_amplicons_with_spacers.append(amplicon)
-            initial_amplicons_with_spacers.append(spacers.pop(0))
-        initial_amplicons = initial_amplicons_with_spacers
+            inputs_withspacers.append(spacers.pop(0))
+        for part in inputs:
+            inputs_withspacers.append(part)
+            inputs_withspacers.append(spacers.pop(0))
+        inputs = inputs_withspacers
         # Maxlink is used to define what is a spacer or what is a template (see docs)
-
-    assembly_amplicons: list[Amplicon] = assembly_fragments(
-        initial_amplicons, overlap=homology_length, circular=circular, maxlink=maxlink
+    inputs = [i for i in inputs if len(i) > 0]
+    assembly_output: list[Amplicon] = assembly_fragments(
+        inputs, overlap=homology_length, circular=circular, maxlink=maxlink
     )
 
-    all_primers = sum((list(amplicon.primers()) for amplicon in assembly_amplicons), [])
+    all_primers = list()
+    for i, part in enumerate(assembly_output):
+        all_primers.extend(list(part.primers() if amplify_templates[i] is True else [None, None]))
 
     for i in range(0, len(all_primers), 2):
         fwd, rvs = all_primers[i : i + 2]
+        if fwd is None or rvs is None:
+            continue
         template = templates[i // 2]
         template_name = template.name if template.name != 'name' else f'seq_{template.id}'
         fwd.name = f'{template_name}_fwd'
         rvs.name = f'{template_name}_rvs'
 
-    return [PrimerModel(id=0, name=primer.name, sequence=str(primer.seq)) for primer in all_primers]
+    return all_primers
 
 
 def sanitize_enzyme_site(site: str) -> str:
@@ -158,7 +199,7 @@ def simple_pair_primers(
     right_enzyme_inverted: bool = False,
     tm_func: Callable[[str], float] = default_tm_func,
     estimate_function: Callable[[str], float] | None = None,
-) -> tuple[PrimerModel, PrimerModel]:
+) -> tuple[Primer, Primer]:
     """
     Design primers to amplify a DNA fragment, if left_enzyme or right_enzyme are set, the primers will be designed
     to include the restriction enzyme sites.
@@ -202,10 +243,7 @@ def simple_pair_primers(
     fwd_primer_name = f'{template_name}_{left_enzyme}_fwd' if left_enzyme is not None else f'{template_name}_fwd'
     rvs_primer_name = f'{template_name}_{right_enzyme}_rvs' if right_enzyme is not None else f'{template_name}_rvs'
 
-    return (
-        PrimerModel(id=0, name=fwd_primer_name, sequence=str(fwd_primer_seq)),
-        PrimerModel(id=0, name=rvs_primer_name, sequence=str(rvs_primer_seq)),
-    )
+    return (Primer(fwd_primer_seq, name=fwd_primer_name), Primer(rvs_primer_seq, name=rvs_primer_name))
 
 
 # def gateway_attB_primers(

src/opencloning/pydantic_models.py

@@ -67,9 +67,11 @@ def all_children_source_ids(self, source_id: int, source_children: list | None =
 class PrimerDesignQuery(BaseModel):
     model_config = {'arbitrary_types_allowed': True}
     sequence: _TextFileSequence
-    location: SequenceLocationStr
+    location: SequenceLocationStr | None
     forward_orientation: bool = True
 
     @field_validator('location', mode='before')
     def parse_location(cls, v):
+        if v is None:
+            return None
         return SequenceLocationStr.field_validator(v)

src/opencloning/temp_functions.py

@@ -48,3 +48,7 @@ def restriction_sequence_cut_to_cutsite_tuple(
 
 def primer_model_to_pydna_primer(primer_model: PrimerModel) -> PydnaPrimer:
     return PydnaPrimer(primer_model.sequence, id=str(primer_model.id), name=primer_model.name)
+
+
+def pydna_primer_to_primer_model(pydna_primer: PydnaPrimer) -> PrimerModel:
+    return PrimerModel(id=0, name=pydna_primer.name, sequence=str(pydna_primer.seq))