#!/usr/bin/env python3
"""C4-adaptive residue scan: palm PEPC vs maize C4-PEPC1 (P04711).
Verifies absence of Ser774 across all palm PEPC isoforms.
MBE manuscript §4 / Supplementary Analysis — generated 2026-07-07."""

from Bio import Align
import re

# Load sequences
with open("all_palms_pepc_v2_dedup.faa") as f:
    raw = f.read()

# Parse FASTA
seqs = {}
current_id = None
current_seq = []
for line in raw.split('\n'):
    if line.startswith('>'):
        if current_id:
            seqs[current_id] = ''.join(current_seq)
        current_id = line[1:].strip()
        current_seq = []
    else:
        current_seq.append(line.strip())
if current_id:
    seqs[current_id] = ''.join(current_seq)

# Identify palm sequences
palm_sps = ['Cocos','Elaeis','Phoenix','Nypa','Calamus','Metroxylon',
            'Areca','Euterpe','Chamaedorea','Borassus']
palm_seqs = {k:v for k,v in seqs.items() 
             if any(sp in k for sp in palm_sps)}

# Filter >= 200 aa
palm_long = {k:v for k,v in palm_seqs.items() if len(v) >= 200}

print(f"Palm PEPC sequences: {len(palm_seqs)}")
print(f"Palm >= 200 aa: {len(palm_long)}")

# C4 reference: P04711 (ZmPEPC1, 970 aa)
ref_id = "sp|P04711|CAPP1_MAIZE"
ref_seq = seqs.get(ref_id)
if not ref_seq:
    raise ValueError("P04711 not found in FASTA")

# Align each palm sequence to reference
aligner = Align.PairwiseAligner()
aligner.substitution_matrix = Align.substitution_matrices.load("BLOSUM62")
aligner.open_gap_score = -10
aligner.extend_gap_score = -0.5

ser774_position = 774  # 1-indexed in P04711

results = []
for seq_id, palm_seq in palm_long.items():
    alignments = aligner.align(ref_seq, palm_seq)
    best = alignments[0]
    aligned_ref, aligned_palm = best.alignments[0]
    
    # Find Ser774 in aligned reference
    ref_pos = 0
    ser774_found = False
    for i, aa in enumerate(aligned_ref):
        if aa != '-':
            ref_pos += 1
        if ref_pos == ser774_position:
            palm_aa = aligned_palm[i]
            ser774_found = (palm_aa == 'S')
            break
    
    results.append((seq_id, ser774_found, palm_aa if ser774_found else None))

# Report
ser774_count = sum(1 for _, found, _ in results if found)
print(f"\nSer774 scan: {ser774_count}/{len(results)} carry Ser774")
print(f"Conclusion: {'No' if ser774_count == 0 else 'SOME'} palm PEPC carries C4 diagnostic Ser774")
for seq_id, found, aa in results:
    if found:
        print(f"  FOUND: {seq_id}")

# Also check N-terminal SIDAQLR motif
sidaqlr = "SIDAQLR"
for seq_id, seq in palm_long.items():
    if sidaqlr in seq[:50]:
        print(f"  SIDAQLR found in N-term of {seq_id}")