pydna_helpers.py

from Bio import SeqIO
from Bio.Seq import Seq
from Bio.SeqFeature import FeatureLocation, SeqFeature
from Bio.SeqRecord import SeqRecord
from pydna.assembly2 import homologous_recombination_integration
from pydna.dseqrecord import Dseqrecord
from pydna.genbank import Genbank


def fetch_genome(accession, email):
    "Fetches genome from GenBank"
    gb = Genbank(email)
    genome = gb.nucleotide(accession)  # Returns Dseqrecord (pydna)
    return genome


def get_plus_strand_cds_stops(genome):
    "Defines non-disruptive insertion positions by identifying CDS sequences and stop codons (for +strand CDS features) - potential for upgrade with (-) strands too"
    stops = []
    for seq in genome.features:
        if seq.type == "CDS" and int(getattr(seq.location, "strand", 0) or 0) == 1:
            stops.append(int(seq.location.end))
    return stops


def extract_homology_arms(genome, insert_pos, left_len, right_len):
    "Extracts left and right homology arms around insertion position"
    if insert_pos < left_len or insert_pos + right_len > len(genome.seq):
        return None
    left_h_arm = str(genome.seq[insert_pos - left_len : insert_pos])
    right_h_arm = str(genome.seq[insert_pos : insert_pos + right_len])
    return left_h_arm, right_h_arm


def count_exact_matches(haystack, needle):
    """ "Counts overlapping occurrences using a 'needle in haystack' approach:
    haystack = genome; needle = left/right_h_arm
    """

    count = 0
    start = 0
    while True:
        pos = haystack.find(needle, start)
        if pos == -1:
            break
        count += 1
        start = pos + 1
    return count


def score_homology_arm_pydna(
    genome_dsr, arm_str, use_subsets: bool, subset_k: int, left_len: int, right_len: int
):
    #### Scores a homology arm using pydna's dsDNA functions
    total = 0

    # Create Dseqrecord from arm string for pydna operations
    arm_dsr = Dseqrecord(arm_str, circular=False)
    genome_str = str(genome_dsr.seq)

    # Count full-length matches
    # Forward strand
    total += count_exact_matches(genome_str, arm_str)

    # Reverse complementary sequence counts
    # A recombination event could use the plus strand or the minus strand as template
    # So we must check both directions to know if the homology arm is unique
    arm_rc_str = str(arm_dsr.rc().seq)
    total += count_exact_matches(genome_str, arm_rc_str)

    # But if there are many shorter perfect matches scattered across the genome,
    # there’s a chance of unintended recombination or micro-homology-mediated events.

    # Count subset matches if "True"
    if use_subsets and 0 < subset_k <= len(arm_str):
        for seq in (arm_str, arm_rc_str):
            for i in range(0, len(seq) - subset_k + 1):
                total += (subset_k / min(left_len, right_len)) * count_exact_matches(
                    str(genome_dsr.seq), seq[i : i + subset_k]
                )
                # multiply by a coefficient because finding a subset match is better than finding an exact one

    return total


def score_candidate(
    genome_dsr,
    left_h,
    right_h,
    count_subsets: bool,
    subset_k: int,
    left_len: int,
    right_len: int,
):
    ### Score a candidate by summing scores of both homology arms using pydna
    return score_homology_arm_pydna(
        genome_dsr=genome_dsr,
        arm_str=left_h,
        use_subsets=count_subsets,
        subset_k=subset_k,
        left_len=left_len,
        right_len=right_len,
    ) + score_homology_arm_pydna(
        genome_dsr=genome_dsr,
        arm_str=right_h,
        use_subsets=count_subsets,
        subset_k=subset_k,
        left_len=left_len,
        right_len=right_len,
    )


def create_candidate_dseqrecord(insert_seq, insert_pos, left, right, score, rank):
    ### Create pydna Dseqrecord for a candidate with full insert_w_homologies sequence
    insert_w_homologies_str = left + insert_seq + right

    # Create Dseqrecord
    candidate = Dseqrecord(insert_w_homologies_str, circular=False)
    candidate.id = "candidate_%03d_pos_%d" % (rank, insert_pos)
    candidate.name = "cand_%03d" % rank
    candidate.description = "Candidate %d at position %d, score=%d" % (
        rank,
        insert_pos,
        score,
    )
    left_len = len(left)
    ins_len = len(insert_seq)
    total_len = len(insert_w_homologies_str)

    # Add features using pydna format
    candidate.features = [
        SeqFeature(
            FeatureLocation(0, left_len),
            type="misc_feature",
            qualifiers={
                "note": [
                    "homology_left (pos %d-%d)" % (insert_pos - left_len, insert_pos)
                ]
            },
        ),
        SeqFeature(
            FeatureLocation(left_len, left_len + ins_len),
            type="misc_feature",
            qualifiers={"note": ["insert"]},
        ),
        SeqFeature(
            FeatureLocation(left_len + ins_len, total_len),
            type="misc_feature",
            qualifiers={
                "note": [
                    "homology_right (pos %d-%d)" % (insert_pos, insert_pos + len(right))
                ]
            },
        ),
        SeqFeature(
            FeatureLocation(0, 0),
            type="misc_feature",
            qualifiers={
                "note": [
                    "score=%d" % score,
                    "rank=%d" % rank,
                    "insertion_pos=%d" % insert_pos,
                ]
            },
        ),
    ]
    return candidate


def create_insert_with_homologies(insert_seq, left, right, rank, pos):
    """Create clean pydna Dseqrecord for insert with homology arms"""
    insert_w_homologies_str = left + insert_seq + right

    insert_dsr = Dseqrecord(insert_w_homologies_str, circular=False)
    insert_dsr.id = "insert_w_homologies_candidate_%03d" % rank
    insert_dsr.name = "insert_%03d" % rank
    insert_dsr.description = (
        "Insert with homologies for candidate %d at position %d" % (rank, pos)
    )

    insert_dsr.features = [
        SeqFeature(
            FeatureLocation(0, len(left)),
            type="misc_feature",
            qualifiers={"note": ["homology_left"]},
        ),
        SeqFeature(
            FeatureLocation(len(left), len(left) + len(insert_seq)),
            type="misc_feature",
            qualifiers={"note": ["insert"]},
        ),
        SeqFeature(
            FeatureLocation(len(left) + len(insert_seq), len(insert_w_homologies_str)),
            type="misc_feature",
            qualifiers={"note": ["homology_right"]},
        ),
    ]

    return insert_dsr


def integrate_features_pydna(genome, insert_pos: int, insert_len: int):
    """Create transformed genome with features shifted around insertion site using pydna"""
    # Create new Dseqrecord with same sequence (will be updated after integration)
    new_genome = Dseqrecord(str(genome.seq), circular=genome.circular)
    new_genome.id = genome.id
    new_genome.name = genome.name
    new_genome.description = genome.description
    new_genome.annotations = dict(genome.annotations)
    new_feats = []
    for feat in genome.features:

        s = int(feat.location.start)
        e = int(feat.location.end)
        strand = int(getattr(feat.location, "strand", 0) or 0)
        # direction of new features is unknown

        if e <= insert_pos:  # feature entirely upstream of insertion site
            loc = FeatureLocation(s, e, strand=strand)  # so we do nothing
        elif (
            s >= insert_pos
        ):  # feature is after the insertion site, so everything is moved
            loc = FeatureLocation(s + insert_len, e + insert_len, strand=strand)
        else:
            loc = FeatureLocation(s, e, strand=strand)
        new_feats.append(
            SeqFeature(
                location=loc,
                type=feat.type,
                id=feat.id,
                qualifiers=dict(feat.qualifiers),
            )
        )

    # Add insertion site markers, adding note in "qualifiers" for a genome positions of insertion
    new_feats.append(
        SeqFeature(
            FeatureLocation(insert_pos, insert_pos),
            type="misc_feature",
            qualifiers={"note": ["insertion_site_after_stop_codon"]},
        )
    )
    new_feats.append(
        SeqFeature(
            FeatureLocation(insert_pos, insert_pos + insert_len),
            type="misc_feature",
            qualifiers={"note": ["insert"]},
        )
    )
    new_genome.features = new_feats

    return new_genome


def write_gbk_pydna(dseqrecord, path):
    """Write a pydna Dseqrecord to a GenBank file.

    Parameters
    ----------
    dseqrecord : Dseqrecord
        The Dseqrecord to write.
    path : str
        The file path to write the GenBank file to.
    """
    # Convert Dseqrecord to BioPython SeqRecord for writing
    sr = SeqRecord(
        Seq(str(dseqrecord.seq)),
        id=dseqrecord.id,
        name=dseqrecord.name,
        description=dseqrecord.description,
        annotations=dseqrecord.annotations,
    )
    sr.features = dseqrecord.features
    sr.annotations["molecule_type"] = sr.annotations.get("molecule_type", "DNA")
    if not sr.id:
        sr.id = "record"
    if not sr.name:
        sr.name = sr.id[:16]
    SeqIO.write(sr, path, "genbank")


def build_and_score_candidates(
    genome,
    insert_seq,
    stops,
    left_len,
    right_len,
    count_subsets: bool = True,
    subset_k: int = 50,
):
    # Getting the list of all arms candidates with their scores
    # the list of stops is a result of function get_plus_strand_cds_stops()
    candidates = []  # a list for storing potential arms
    for pos in stops:
        arms = extract_homology_arms(
            genome=genome, insert_pos=pos, left_len=left_len, right_len=right_len
        )
        if not arms:
            continue
        left_h, right_h = arms
        sc = score_candidate(
            genome_dsr=genome,
            left_h=left_h,
            right_h=right_h,
            subset_k=subset_k,
            count_subsets=count_subsets,
            left_len=left_len,
            right_len=right_len,
        )
        insert_w_homologies_str = left_h + insert_seq + right_h
        insert_w_dsr = Dseqrecord(insert_w_homologies_str, circular=False)
        candidates.append(
            {
                "pos": pos,
                "left": left_h,
                "right": right_h,
                "score": sc,
                "insert_w_homologies": insert_w_dsr,
            }
        )
    return candidates


def selective_homology_integration(
    genome, ranked_candidates, limit_overlap=None, max_tries=50, progress_every=25
):
    ### Try to integrate a donor DNA with homology arms into a genome,
    ### testing candidate insertion sites one by one until we find a unique (single-product) homologous recombination.
    if not ranked_candidates:
        return None, []
    left_len = len(
        ranked_candidates[0]["left"]
    )  # ranked_candidates is a (dict), result of build_and_score_candidates() ranked
    right_len = len(ranked_candidates[0]["right"])
    if limit_overlap is None:
        limit_overlap = min(left_len, right_len)
    genome_linear = Dseqrecord(str(genome.seq), circular=False)  # making genome linear
    tried = []
    unique_hits = []
    for idx, cand in enumerate(ranked_candidates[:max_tries], 1):
        insert_w_dsr = cand["insert_w_homologies"]  # left_h + INSERT_SEQ + right_h
        try:
            products = homologous_recombination_integration(
                genome_linear, [insert_w_dsr], limit=limit_overlap
            )
        except Exception:
            products = []
        tried.append({**cand, "products": products})
        if (
            products and len(products) == 1
        ):  # exactly one recombination product means a unique, clean integration
            unique_hits.append({**cand, "product": products[0]})
            break
        if idx % progress_every == 0:
            print(
                "   Tried integration on %d / %d top sites..."
                % (idx, min(max_tries, len(ranked_candidates)))
            )
    best_result = unique_hits[0] if unique_hits else None
    return best_result, tried