CRISPR Guide RNA Design: Best Practices for Precision Gene Editing

class CRISPRSystem:
    """
    CRISPR-Cas9 gene editing system.

    Components:
    - Cas9: Enzyme that cuts DNA
    - Guide RNA: 20-nucleotide sequence that directs Cas9 to target
    - PAM sequence: NGG motif required next to target (Cas9 recognition)
    """

    def __init__(self, target_sequence, pam='NGG'):
        self.target = target_sequence  # 20 nt guide sequence
        self.pam = pam  # Protospacer Adjacent Motif
        self.cas9 = Cas9Enzyme()

    def cut_dna(self, genome):
        """
        Cut DNA at target location.

        Process:
        1. Guide RNA binds to complementary DNA sequence
        2. Cas9 checks for PAM sequence (NGG)
        3. If match: Cas9 cuts both DNA strands
        4. Cell repairs cut (can introduce edits)
        """
        # Find target in genome
        cut_site = genome.find(self.target + self.pam)

        if cut_site == -1:
            return None  # Target not found

        # Create double-strand break
        genome = genome[:cut_site] + '[CUT]' + genome[cut_site:]

        return genome

import re
from Bio import SeqIO
from Bio.Seq import Seq

class GuideRNADesigner:
    """Design optimal guide RNAs for CRISPR editing."""

    def __init__(self, gene_sequence, edit_location):
        """
        gene_sequence: Full gene sequence (DNA)
        edit_location: Base pair position to edit
        """
        self.sequence = gene_sequence
        self.edit_location = edit_location

    def find_all_pam_sites(self):
        """
        Find all possible PAM sites (NGG) near edit location.

        PAM can be NGG where N is any nucleotide:
        AGG, TGG, CGG, GGG
        """
        pam_pattern = r'[ATCG]GG'  # Regex for NGG

        # Search in 100bp window around edit location
        window_start = max(0, self.edit_location - 50)
        window_end = min(len(self.sequence), self.edit_location + 50)

        search_region = self.sequence[window_start:window_end]

        # Find all PAM sites
        pam_sites = []
        for match in re.finditer(pam_pattern, search_region):
            pam_position = window_start + match.start()
            pam_sites.append({
                'position': pam_position,
                'pam_sequence': match.group(),
                'distance_to_edit': abs(pam_position - self.edit_location)
            })

        return sorted(pam_sites, key=lambda x: x['distance_to_edit'])

    def design_guide_rna(self, pam_site):
        """
        Design 20nt guide RNA targeting upstream of PAM.

        Guide RNA: 20 nucleotides immediately 5' of PAM
        PAM: NGG motif (not included in guide RNA)
        """
        pam_position = pam_site['position']

        # Extract 20nt upstream of PAM
        guide_start = pam_position - 20
        guide_end = pam_position

        guide_rna = self.sequence[guide_start:guide_end]

        return {
            'guide_sequence': guide_rna,
            'pam': pam_site['pam_sequence'],
            'start': guide_start,
            'end': guide_end,
            'full_target': guide_rna + pam_site['pam_sequence']
        }

    def score_guide_rna(self, guide_rna):
        """
        Score guide RNA for on-target efficiency.

        Based on Doench et al. 2016 scoring algorithm.
        Factors:
        - GC content (40-60% optimal)
        - Position-specific nucleotide preferences
        - Avoiding poly-T sequences (terminates transcription)
        """
        sequence = guide_rna['guide_sequence']

        # 1. GC content (optimal: 40-60%)
        gc_count = sequence.count('G') + sequence.count('C')
        gc_content = gc_count / len(sequence)

        if 0.4 <= gc_content <= 0.6:
            gc_score = 1.0
        else:
            gc_score = 0.5

        # 2. Poly-T check (TTTT terminates pol III transcription)
        if 'TTTT' in sequence:
            poly_t_score = 0.0  # Fail
        else:
            poly_t_score = 1.0

        # 3. Position-specific preferences (simplified)
        position_score = self._calculate_position_score(sequence)

        # Combined score
        total_score = (gc_score + poly_t_score + position_score) / 3

        return {
            'total_score': total_score,
            'gc_content': gc_content,
            'has_poly_t': 'TTTT' in sequence,
            'predicted_efficiency': total_score * 100  # Percentage
        }

    def _calculate_position_score(self, sequence):
        """
        Position-specific nucleotide scoring.

        Certain positions prefer certain nucleotides for efficiency.
        """
        # Simplified scoring (real version uses ML model)
        score = 0

        # Position 1: prefer G
        if sequence[0] == 'G':
            score += 0.1

        # Position 20: avoid T
        if sequence[19] != 'T':
            score += 0.1

        # Middle positions: balanced
        middle = sequence[7:13]
        if 2 <= middle.count('G') + middle.count('C') <= 4:
            score += 0.1

        return min(score, 1.0)

class OffTargetPredictor:
    """
    Predict off-target binding sites for guide RNA.

    ⚠️ CRITICAL: Off-target effects can edit unintended genes!
    """

    def __init__(self, genome_fasta):
        """
        genome_fasta: Path to reference genome FASTA file
        """
        self.genome = self._load_genome(genome_fasta)

    def _load_genome(self, fasta_path):
        """Load genome from FASTA file."""
        genome_sequence = ""
        for record in SeqIO.parse(fasta_path, "fasta"):
            genome_sequence += str(record.seq)
        return genome_sequence

    def find_off_targets(self, guide_rna, max_mismatches=3):
        """
        Find potential off-target sites in genome.

        CRISPR tolerates 1-4 mismatches and still cuts!

        max_mismatches: Maximum allowed mismatches (default: 3)
        """
        guide_seq = guide_rna['guide_sequence']
        pam = guide_rna['pam']

        potential_off_targets = []

        # Search entire genome for similar sequences
        # (In production: use Bowtie2 or BLAST for speed)
        for i in range(len(self.genome) - 23):
            # Extract 20nt + PAM
            candidate = self.genome[i:i+20]
            candidate_pam = self.genome[i+20:i+23]

            # Check PAM matches
            if not self._pam_matches(candidate_pam, pam):
                continue

            # Count mismatches
            mismatches = self._count_mismatches(guide_seq, candidate)

            if mismatches <= max_mismatches:
                potential_off_targets.append({
                    'position': i,
                    'sequence': candidate,
                    'mismatches': mismatches,
                    'pam': candidate_pam,
                    'cutting_probability': self._calculate_cutting_prob(mismatches)
                })

        return sorted(potential_off_targets, key=lambda x: x['mismatches'])

    def _count_mismatches(self, seq1, seq2):
        """Count number of mismatched nucleotides."""
        return sum(a != b for a, b in zip(seq1, seq2))

    def _pam_matches(self, candidate_pam, target_pam):
        """Check if PAM sequence is compatible."""
        # NGG pattern allows any first nucleotide
        return candidate_pam[1:] == 'GG'

    def _calculate_cutting_prob(self, mismatches):
        """
        Estimate probability of off-target cutting.

        0 mismatches: ~100% cutting
        1 mismatch: ~40% cutting
        2 mismatches: ~10% cutting
        3 mismatches: ~2% cutting
        4+ mismatches: <1% cutting
        """
        probabilities = {
            0: 1.0,
            1: 0.4,
            2: 0.1,
            3: 0.02,
            4: 0.005
        }
        return probabilities.get(mismatches, 0.001)

def design_optimal_guide_rna(gene_sequence, edit_position, genome_fasta):
    """
    Complete guide RNA design workflow with off-target screening.

    Returns: Best guide RNA with safety assessment
    """

    # Step 1: Find all possible guides near edit position
    designer = GuideRNADesigner(gene_sequence, edit_position)
    pam_sites = designer.find_all_pam_sites()

    if not pam_sites:
        raise ValueError("No PAM sites found near edit position")

    # Step 2: Design guide RNAs for all PAM sites
    guide_candidates = []
    for pam_site in pam_sites[:10]:  # Top 10 closest PAMs
        guide = designer.design_guide_rna(pam_site)
        score = designer.score_guide_rna(guide)

        guide['on_target_score'] = score['total_score']
        guide['predicted_efficiency'] = score['predicted_efficiency']

        guide_candidates.append(guide)

    # Step 3: Screen for off-targets
    predictor = OffTargetPredictor(genome_fasta)

    for guide in guide_candidates:
        off_targets = predictor.find_off_targets(guide, max_mismatches=3)
        guide['off_targets'] = off_targets
        guide['off_target_count'] = len(off_targets)

        # Calculate specificity score (penalize off-targets)
        guide['specificity_score'] = 1.0 / (1.0 + len(off_targets))

    # Step 4: Rank by combined score
    for guide in guide_candidates:
        guide['combined_score'] = (
            guide['on_target_score'] * 0.6 +
            guide['specificity_score'] * 0.4
        )

    # Sort by combined score
    ranked_guides = sorted(
        guide_candidates,
        key=lambda x: x['combined_score'],
        reverse=True
    )

    best_guide = ranked_guides[0]

    # Step 5: Safety assessment
    safety_report = assess_safety(best_guide)

    return best_guide, safety_report

def assess_safety(guide_rna):
    """
    Assess safety of guide RNA.

    ⚠️ WARNINGS to check:
    - Off-target sites in critical genes
    - High cutting probability off-targets
    - Potential horizontal gene transfer
    """

    warnings = []

    # Check 1: Off-target count
    if guide_rna['off_target_count'] > 5:
        warnings.append({
            'severity': 'HIGH',
            'message': f"{guide_rna['off_target_count']} potential off-targets detected"
        })

    # Check 2: High-probability off-targets
    high_prob_off_targets = [
        ot for ot in guide_rna['off_targets']
        if ot['cutting_probability'] > 0.1
    ]
    if high_prob_off_targets:
        warnings.append({
            'severity': 'CRITICAL',
            'message': f"{len(high_prob_off_targets)} off-targets with >10% cutting probability"
        })

    # Check 3: PAM density (for gene drive risk)
    if guide_rna['guide_sequence'].count('GG') > 3:
        warnings.append({
            'severity': 'MEDIUM',
            'message': "High GG content - potential gene drive substrate"
        })

    return {
        'approved': len([w for w in warnings if w['severity'] == 'CRITICAL']) == 0,
        'warnings': warnings,
        'recommendation': 'APPROVE' if not warnings else 'REVIEW_REQUIRED'
    }

# Example usage
gene_seq = "ATCGATCGATCG..."  # Your target gene
edit_pos = 500  # Base pair to edit
genome = "path/to/human_genome.fasta"

best_guide, safety = design_optimal_guide_rna(gene_seq, edit_pos, genome)

print(f"Best guide RNA: {best_guide['guide_sequence']}")
print(f"On-target efficiency: {best_guide['predicted_efficiency']:.1f}%")
print(f"Off-targets found: {best_guide['off_target_count']}")
print(f"Safety assessment: {safety['recommendation']}")

class BaseEditor:
    """
    Base editors: Precise single-nucleotide changes without double-strand breaks.

    Types:
    - CBE (Cytosine Base Editor): C → T conversions
    - ABE (Adenine Base Editor): A → G conversions
    """

    def __init__(self, editor_type='CBE'):
        self.editor_type = editor_type

    def design_base_edit(self, sequence, target_position, desired_edit):
        """
        Design guide RNA for base editing.

        sequence: Gene sequence
        target_position: Position of base to edit
        desired_edit: e.g., 'C->T' or 'A->G'
        """

        # Base editors have editing window (typically positions 4-8 of guide)
        editing_window = (4, 8)

        # Design guide RNA so target falls in editing window
        # Target should be at position 4-8 of the guide (counting from 5' end)

        # Calculate where guide should start
        guide_start = target_position - editing_window[0]
        guide_end = guide_start + 20

        guide_rna = sequence[guide_start:guide_end]

        # Find PAM (must be 3' of guide)
        pam_start = guide_end
        pam = sequence[pam_start:pam_start+3]

        if pam[1:] != 'GG':
            raise ValueError(f"No valid PAM found. Got: {pam}")

        return {
            'guide_sequence': guide_rna,
            'pam': pam,
            'edit_position_in_guide': target_position - guide_start,
            'edit_type': desired_edit,
            'editor_type': self.editor_type
        }

# Example: Change C to T at position 450
editor = BaseEditor(editor_type='CBE')
base_edit_guide = editor.design_base_edit(
    sequence=gene_seq,
    target_position=450,
    desired_edit='C->T'
)

print(f"Base editor guide: {base_edit_guide['guide_sequence']}")
print(f"Edit at position {base_edit_guide['edit_position_in_guide']} in guide")

class CRISPRDelivery:
    """Methods to deliver CRISPR into cells."""

    @staticmethod
    def plasmid_delivery(guide_rna_sequence):
        """
        Plasmid-based delivery (research standard).

        Plasmid contains:
        - Cas9 gene
        - Guide RNA gene
        - Selection marker
        """
        plasmid_sequence = f"""
        // Plasmid: pCRISPR-Cas9-{guide_rna_sequence[:10]}
        Origin of replication: pUC ori
        Cas9 gene: Human codon-optimized
        U6 promoter: {guide_rna_sequence}  // Guide RNA expression
        Selection: Ampicillin resistance
        """
        return plasmid_sequence

    @staticmethod
    def viral_vector_delivery(guide_rna):
        """
        AAV (Adeno-Associated Virus) delivery (clinical use).

        ⚠️ WARNING: Viral delivery is permanent and can integrate into genome
        """
        return {
            'vector_type': 'AAV',
            'cargo': 'Cas9 + guide RNA',
            'tropism': 'Liver/muscle/brain (depending on serotype)',
            'integration_risk': 'Low but non-zero',
            'immune_response': 'Possible anti-AAV antibodies'
        }

    @staticmethod
    def rnp_delivery(guide_rna):
        """
        Ribonucleoprotein (RNP) delivery (safest for clinical).

        Cas9 protein + guide RNA delivered directly.
        No DNA integration, transient editing.
        """
        return {
            'components': 'Cas9 protein + synthetic guide RNA',
            'half_life': '24 hours (degrades naturally)',
            'integration_risk': 'Zero (no DNA template)',
            'delivery_method': 'Electroporation or lipid nanoparticles'
        }

def clinical_safety_checklist(guide_rna):
    """
    Safety checks for clinical CRISPR use.

    ⚠️ CRITICAL for human gene therapy
    """

    checks = {
        'off_targets_screened': len(guide_rna.get('off_targets', [])) > 0,
        'specificity_score': guide_rna.get('specificity_score', 0) > 0.9,
        'immune_response_predicted': False,  # Would check guide RNA immunogenicity
        'germline_editing_prevented': True,   # Only edit somatic cells
        'reversibility_planned': False,       # Most CRISPR edits are permanent
        'informed_consent_obtained': False,   # Patient understanding required
        'ethics_approval': False,             # IRB approval
        'regulatory_approval': False          # FDA/EMA approval
    }

    passed = sum(checks.values())
    total = len(checks)

    return {
        'checks_passed': passed,
        'checks_total': total,
        'approval_recommended': passed == total,
        'details': checks
    }

# Example
safety = clinical_safety_checklist(best_guide)
print(f"Safety checks: {safety['checks_passed']}/{safety['checks_total']} passed")

def assess_gene_drive_risk(guide_rna):
    """
    Assess if guide RNA could enable gene drive.

    Gene drives: Self-propagating genetic modifications
    - Guide RNA targets its own insertion site
    - Creates feedback loop
    - Spreads through population

    ⚠️ Can rewrite entire species' genome
    """

    # Check if guide targets genomic region that could harbor CRISPR cassette
    # (Simplified - real analysis more complex)

    risk_score = 0

    # High GC content facilitates insertion
    gc_content = (guide_rna['guide_sequence'].count('G') +
                  guide_rna['guide_sequence'].count('C')) / 20

    if gc_content > 0.6:
        risk_score += 1

    # PAM density affects drive efficiency
    if 'GG' in guide_rna['guide_sequence']:
        risk_score += 1

    # Check for homology to common genomic elements
    # (Would check against database of mobile genetic elements)

    return {
        'gene_drive_risk': 'HIGH' if risk_score >= 2 else 'LOW',
        'risk_score': risk_score,
        'recommendation': 'Additional safety testing required' if risk_score >= 2 else 'Approved'
    }

def predict_evolutionary_escape(guide_rna, generations=100):
    """
    Predict if target organism could evolve resistance.

    Relevant for:
    - Therapeutic CRISPR (cancer cells evolving resistance)
    - Agricultural CRISPR (pests evolving resistance)
    - Gene drives (populations evolving drive-resistance)
    """

    # Simulate evolution
    resistance_probability = 0

    for gen in range(generations):
        # Mutation could alter PAM or guide binding site
        mutation_rate = 1e-8  # Per base pair per generation

        # 23 bp target (20 guide + 3 PAM)
        escape_probability = 1 - (1 - mutation_rate) ** 23

        resistance_probability += escape_probability

    return {
        'generations': generations,
        'resistance_probability': resistance_probability,
        'expected_resistance_time': f"{int(1/resistance_probability)} generations"
    }

# Real-world guide RNA design tools
DESIGN_TOOLS = {
    'Benchling CRISPR': 'https://benchling.com',  # Commercial, web-based
    'CRISPOR': 'http://crispor.org',              # Free, academic
    'Cas-OFFinder': 'http://www.rgenome.net',     # Off-target prediction
    'IDT Custom gRNA': 'https://www.idtdna.com',  # Commercial synthesis
}

# Guide RNA synthesis (commercial services)
def order_guide_rna(sequence):
    """Order synthetic guide RNA from vendor."""
    return {
        'sequence': sequence,
        'modifications': '2\'-O-Methyl, 3\' phosphorothioate (increased stability)',
        'price': '$50-200',
        'turnaround': '3-5 days',
        'vendors': ['IDT', 'Synthego', 'GenScript']
    }