# Python code to generate mock predicted sequences and mutation probability matrices
import numpy as np
import pandas as pd

def generate_mock_protein_data(length=100, num_aa_types=20, wild_type_sequence=None):
    """
    Generates mock protein data: wild-type sequence, predicted sequence,
    and a mutation probability matrix.
    """
    amino_acids = 'ACDEFGHIKLMNPQRSTVWY'
    if wild_type_sequence is None:
        wild_type_sequence = ''.join(np.random.choice(list(amino_acids), length))

    # Simulate a predicted sequence with some mutations
    predicted_sequence_list = list(wild_type_sequence)
    num_mutations = np.random.randint(5, 15) # Random number of mutations
    mutation_indices = np.random.choice(length, num_mutations, replace=False)
    for idx in mutation_indices:
        original_aa = predicted_sequence_list[idx]
        possible_mutations = [aa for aa in amino_acids if aa != original_aa]
        predicted_sequence_list[idx] = np.random.choice(possible_mutations)
    predicted_sequence = ''.join(predicted_sequence_list)

    # Simulate mutation probability matrix (position x AA type)
    # Each row sums to 1 (probability of changing TO a specific AA at that position)
    mutation_probabilities = np.random.rand(length, num_aa_types)
    mutation_probabilities = mutation_probabilities / mutation_probabilities.sum(axis=1, keepdims=True)
    
    # Create a DataFrame for mutation probabilities for easier handling
    mutation_df = pd.DataFrame(
        mutation_probabilities,
        columns=list(amino_acids),
        index=[f'Pos_{i+1}' for i in range(length)]
    )
    
    print(f"Wild-type Sequence: {wild_type_sequence}")
    print(f"Predicted Sequence: {predicted_sequence}")
    print(f"Mutation Probabilities (first 5 positions):
{mutation_df.head()}")
    
    return wild_type_sequence, predicted_sequence, mutation_df

# --- EXECUTION --- 
# Define a custom wild-type sequence for reproducibility or use None for random
custom_wt = "MVLSPADKTNVKAAWGKVGAHAGEYGAEALERMFLSFPTTKTYFPHFDLSHGSAQVKGHGKKVADALTNAVAHVDDMPNALSALSDLLHAHKLRVDPVNFKLLSHCLLVTLAAHLPAEFTPAVHASLDKFLASVSTVLTSKYR"

wild_type, predicted_seq, mutation_probs_df = generate_mock_protein_data(
    length=len(custom_wt),
    wild_type_sequence=custom_wt
)

# Save data for subsequent steps (optional, can pass variables directly)
# np.save('mutation_probabilities.npy', mutation_probs_df.values)
# with open('wild_type_sequence.txt', 'w') as f: f.write(wild_type)
# with open('predicted_sequence.txt', 'w') as f: f.write(predicted_seq)

# Python code to visualize differences between wild-type and predicted sequences
import matplotlib.pyplot as plt
import numpy as np

def visualize_sequence_differences(wild_type_seq, predicted_seq):
    """
    Visualizes differences between two protein sequences.
    Highlights differing amino acids.
    """
    if len(wild_type_seq) != len(predicted_seq):
        raise ValueError("Sequences must be of the same length for comparison.")

    length = len(wild_type_seq)
    differences = []
    for i in range(length):
        if wild_type_seq[i] != predicted_seq[i]:
            differences.append(i) # Store index of difference

    # Create a simple plot to highlight differences
    fig, ax = plt.subplots(figsize=(max(15, length / 5), 3))

    # Plot wild-type sequence
    ax.text(0, 1.5, 'Wild-type:', va='center', ha='right', fontsize=12)
    for i, aa in enumerate(wild_type_seq):
        ax.text(i, 1.5, aa, va='center', ha='center', fontsize=10, color='black')

    # Plot predicted sequence
    ax.text(0, 0.5, 'Predicted:', va='center', ha='right', fontsize=12)
    for i, aa in enumerate(predicted_seq):
        color = 'red' if i in differences else 'black'
        ax.text(i, 0.5, aa, va='center', ha='center', fontsize=10, color=color, fontweight='bold' if i in differences else 'normal')
        if i in differences:
            # Add a small marker or box around the differing residue
            ax.plot(i, 0.5, 'o', color='red', markersize=8, alpha=0.3)

    ax.set_xlim(-1, length)
    ax.set_ylim(0, 2)
    ax.axis('off') # Hide axes
    ax.set_title('Predicted Sequence vs. Wild-type: Mutations Highlighted', fontsize=14)
    plt.tight_layout()
    plt.show()
    
    print(f"Found {len(differences)} differences at positions (0-indexed): {differences}")

# --- EXECUTION --- 
# Use the data generated in the previous step
# wild_type, predicted_seq = 'MVLSPADKT...', 'MVLSKADKT...'
# For demonstration, ensure these variables are available from the previous run

visualize_sequence_differences(wild_type, predicted_seq)

# Python code to generate a heatmap of mutation probabilities
import matplotlib.pyplot as plt
import seaborn as sns
import pandas as pd

def visualize_mutation_heatmap(mutation_df, title="Mutation Probability Heatmap", top_n_aa=20):
    """
    Generates a heatmap of mutation probabilities across sequence positions
    and amino acid types.
    """
    amino_acids = 'ACDEFGHIKLMNPQRSTVWY'
    # Ensure amino acids are in the correct order for consistent plotting
    plot_df = mutation_df.loc[:, [aa for aa in amino_acids if aa in mutation_df.columns]]

    plt.figure(figsize=(max(15, len(mutation_df)/5), 10)) # Adjust figure size dynamically
    sns.heatmap(
        plot_df.T, # Transpose to have AA types on y-axis, positions on x-axis
        cmap='viridis', # Choose a suitable colormap (e.g., 'viridis', 'plasma', 'rocket')
        linewidths=.5, 
        linecolor='lightgray',
        cbar_kws={'label': 'Probability'},
        yticklabels=True, # Show AA labels on y-axis
        xticklabels=10 # Show every 10th position label for readability
    )
    plt.title(title, fontsize=16)
    plt.xlabel('Protein Position', fontsize=12)
    plt.ylabel('Amino Acid Type', fontsize=12)
    plt.yticks(rotation=0) # Keep AA labels horizontal
    plt.tight_layout()
    plt.show()

# --- EXECUTION --- 
# Use the mutation_probs_df generated in the first step
# mutation_probs_df = pd.DataFrame(...)

visualize_mutation_heatmap(mutation_probs_df, title="Predicted Amino Acid Probability Across Sequence Positions")

# Python code to highlight key mutations based on probability thresholds
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np

def visualize_key_mutations(wild_type_seq, mutation_df, threshold=0.3, num_top_mutations=10):
    """
    Identifies and visualizes key mutation sites based on a probability threshold
    and highlights the most likely substitution.
    """
    amino_acids = list('ACDEFGHIKLMNPQRSTVWY')
    
    # Calculate entropy or variance to find positions with high uncertainty/potential for change
    # For simplicity, let's identify positions where the max probability for *any* AA (other than WT) is high.
    
    # Find the wild-type AA at each position from the wild_type_seq
    wt_aas_at_pos = [wild_type_seq[i] for i in range(len(wild_type_seq))]
    
    # For each position, find the highest probability for a *non-wild-type* amino acid
    max_mutation_probs = []
    most_likely_mutations = []
    for i, wt_aa in enumerate(wt_aas_at_pos):
        row_probs = mutation_df.iloc[i].copy() # Get probabilities for current position
        if wt_aa in row_probs.index: # If WT AA is in our AA list, exclude its prob for 'mutation' analysis
            row_probs[wt_aa] = 0 # Set WT AA probability to 0 to find highest *mutation* prob
        
        max_prob_for_pos = row_probs.max() # Max probability among non-WT AAs
        most_likely_aa = row_probs.idxmax() if max_prob_for_pos > 0 else wt_aa # If no mutation prob > 0, default to WT
        
        max_mutation_probs.append(max_prob_for_pos)
        most_likely_mutations.append(most_likely_aa)

    # Filter positions based on threshold
    high_impact_positions = []
    for i, prob in enumerate(max_mutation_probs):
        if prob >= threshold:
            high_impact_positions.append({
                'position': i + 1, # 1-indexed
                'wild_type_aa': wt_aas_at_pos[i],
                'predicted_aa': most_likely_mutations[i],
                'probability': prob
            })
            
    # Sort by probability descending and take top N
    high_impact_positions_sorted = sorted(high_impact_positions, key=lambda x: x['probability'], reverse=True)[:num_top_mutations]

    if not high_impact_positions_sorted:
        print("No high-impact mutations found above the given threshold.")
        return

    # Plotting the top mutations
    positions = [d['position'] for d in high_impact_positions_sorted]
    probabilities = [d['probability'] for d in high_impact_positions_sorted]
    labels = [f"Pos {d['position']}: {d['wild_type_aa']}->{d['predicted_aa']}" for d in high_impact_positions_sorted]

    fig, ax = plt.subplots(figsize=(12, 6))
    bars = ax.bar(range(len(positions)), probabilities, color='skyblue')
    ax.set_xticks(range(len(positions)))
    ax.set_xticklabels(labels, rotation=45, ha='right')
    ax.set_ylabel('Probability of Most Likely Mutation')
    ax.set_title(f'Top {len(high_impact_positions_sorted)} High-Impact Mutation Sites (Probability > {threshold})')
    ax.set_ylim(0, 1)

    for bar in bars:
        yval = bar.get_height()
        ax.text(bar.get_x() + bar.get_width()/2, yval + 0.02, round(yval, 2), ha='center', va='bottom', fontsize=9)

    plt.tight_layout()
    plt.show()

    print("\nHigh-Impact Mutations Detected:")
    for h_imp in high_impact_positions_sorted:
        print(f" - Position {h_imp['position']} (Wild-type: {h_imp['wild_type_aa']}): "
              f"Most likely mutation to {h_imp['predicted_aa']} (Probability: {h_imp['probability']:.2f})")


# --- EXECUTION --- 
# Use wild_type and mutation_probs_df from previous steps
# wild_type = 'MVLSPADKT...'
# mutation_probs_df = pd.DataFrame(...)

visualize_key_mutations(wild_type, mutation_probs_df, threshold=0.15, num_top_mutations=15)