PCA_complex_hulls.py

import json
import pandas as pd
import copy
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler
import matplotlib.pyplot as plt
from matplotlib.patches import Circle
from scipy.spatial import ConvexHull
import numpy as np
from sklearn.cluster import KMeans
from sklearn.metrics import silhouette_samples, silhouette_score
import matplotlib.cm as cm


colors = {'Cell_control_rat': '#da3500',
          'Cell_control_A': '#e66101',
          'Cell_control_B': '#f28d32',
          'Cell_control_D': '#fdb863',
          'Cell_control_dSPN' : '#e66101',
          'Cell_control_iSPN' : '#fdb863',
          'Cell_control_SPN': "#da3500",
          "Cell_control": '#da3500',
          'Cell_HT': '#1a9850',
          'Cell_MCh':'#7f46c2',
          'Cell_ACh':'#7f46c2',
          'Cell_SPN_ACh': '#7f46c2',
          'Cell_NA': '#2166ac',
          'Cell_5HT':'#6ce058',
          'Cell_DA': '#0c89f0',
          'Cell_Ox': '#0c89f0',
          'Cell_DA_rat': '#0c89f0',
          'Cell_human_DA': '#9bc1e0',
          'Cell_d1_DA': '#0c89f0',
          'Cell_d2_DA': '#23e0ff',
          'Cell_avg_human': '#ce0900'}

colors_tot = {'Cell_ctx_control_rat': '#da3500',
          'Cell_ctx_control_A': '#e66101',
          'Cell_ctx_control_B': '#f28d32',
          'Cell_ctx_control_D': '#fdb863',
          'Cell_control_dSPN' : '#e66101',
          'Cell_control_iSPN' : '#fdb863',
          'Cell_control_SPN': "#da3500",
          'Cell_control': '#da3500',
          'Cell_DG_control': '#da3500',
          'Cell_TPS_control': '#da3500',
          'Cell_TRN_control': '#da3500', 
          'Cell_GC_control': '#da3500',    
          'Cell_ctx_HT': '#1a9850',
          'Cell_DG_HT': '#1a9850',    
          'Cell_ctx_MCh':'#7f46c2',
          'Cell_ACh':'#7f46c2',
          'Cell_TPS_ACh':'#7f46c2',    
          'Cell_DG_ACh':'#7f46c2',
          'Cell_TRN_ACh':'#7f46c2',
          'Cell_GC_ACh':'#7f46c2',     
          'Cell_SPN_ACh': '#7f46c2',
          'Cell_ctx_NA': '#2166ac',
          'Cell_TPS_NA': '#2166ac',   
          'Cell_DG_NA': '#2166ac',    
          'Cell_ctx_5HT':'#6ce058',
          'Cell_GC_5HT':'#6ce058',    
          'Cell_DA': '#0c89f0',
          'Cell_TPS_Ox': '#0c89f0',
          'Cell_ctx_DA_rat': '#0c89f0',
          'Cell_human_DA': '#9bc1e0',
          'Cell_d1_DA': '#0c89f0',
          'Cell_d2_DA': '#23e0ff',
          'Cell_avg_human': '#ce0900'}

colors_region = {'Cell_ctx_control_rat': '#d7191c',
          'Cell_ctx_control_A': '#d7191c',
          'Cell_ctx_control_B': '#d7191c',
          'Cell_ctx_control_D': '#d7191c',
          'Cell_control_dSPN' : '#fdae61',
          'Cell_control_iSPN' : '#fdae61',
          'Cell_control_SPN': "#fdae61",
          'Cell_control': '#da3500',
          'Cell_DG_control': '#ffffbf',
          'Cell_TPS_control': '#abd9e9',
          'Cell_TRN_control': '#abd9e9', 
          'Cell_GC_control': '#2c7bb6',    
          'Cell_ctx_HT': '#d7191c',
          'Cell_DG_HT': '#ffffbf',    
          'Cell_ctx_MCh':'#d7191c',
          'Cell_ACh':'#7f46c2',
          'Cell_TPS_ACh':'#abd9e9',    
          'Cell_DG_ACh':'#ffffbf',
          'Cell_TRN_ACh':'#abd9e9',
          'Cell_GC_ACh':'#2c7bb6',     
          'Cell_SPN_ACh': '#fdae61',
          'Cell_ctx_NA': '#d7191c',
          'Cell_TPS_NA': '#abd9e9',   
          'Cell_DG_NA': '#ffffbf',    
          'Cell_ctx_5HT':'#d7191c',
          'Cell_GC_5HT':'#2c7bb6',    
          'Cell_DA': '#0c89f0',
          'Cell_TPS_Ox': '#abd9e9',
          'Cell_ctx_DA_rat': '#d7191c',
          'Cell_human_DA': '#9bc1e0',
          'Cell_d1_DA': '#fdae61',
          'Cell_d2_DA': '#fdae61',
          'Cell_avg_human': '#ce0900'}

# Extract_model_features function to include the model ID
def extract_model_features(cell_data, cell_name, n_models):
    # Extract model features from each dictionary and assign the IDs
    models = [entry['model'] for entry in cell_data][:n_models]
    df = pd.DataFrame(models)
    df['cell_type'] = cell_name  
    df['model_id'] = [f"{i}" for i in range(len(df))]
    return df


def running_pca(data):
    # Separate the cell_type and model_id columns before standardizing the data
    cell_types = data['cell_type']
    model_ids = data['model_id']
    data_for_pca = data.drop(columns=['cell_type', 'model_id'])

    print("\nParameters variance:")
    print(np.var(data_for_pca, axis=0))
    
    # Standardize the features and perform PCA
    scaler = StandardScaler()
    standardized_data = scaler.fit_transform(data_for_pca)

    pca = PCA(n_components=2)
    pca_result = pca.fit_transform(standardized_data)
    
    # Create a DataFrame for PCA result
    pca_df = pd.DataFrame(pca_result, columns=['PC1', 'PC2'])
    pca_df['cell_type'] = cell_types
    pca_df['model_id'] = model_ids  # Add back the model IDs

    # Print out the top contributing features for PC1 and PC2
    feature_names = data_for_pca.columns
    loading_vectors = pca.components_.T  # shape (n_features, n_components)
    loading_df = pd.DataFrame(loading_vectors, index=feature_names, columns=['PC1', 'PC2'])
    
    print("\nTop features influencing PC1:")
    print(loading_df['PC1'].abs().sort_values(ascending=False).head(5))
    print("\nTop features influencing PC2:")
    print(loading_df['PC2'].abs().sort_values(ascending=False).head(5))
    # After fitting the PCA:
    explained_variance = pca.explained_variance_ratio_
    
    print(f"\nExplained variance by PC1: {explained_variance[0]:.2%}")
    print(f"Explained variance by PC2: {explained_variance[1]:.2%}")

    required_columns = {'PC1', 'PC2', 'cell_type'}
    if not required_columns.issubset(pca_df.columns):
        raise ValueError(f"pca_df must contain the following columns: {required_columns}")
    
    pca_points = pca_df[['PC1', 'PC2']].values  
    labels = pca_df['cell_type'].values         
    
    # Calculate the silhouette score
    try:
        silhouette = silhouette_score(pca_points, labels)
        print(f"\nSilhouette Score: {silhouette:.4f}")
    except Exception as e:
        print(f"Silhouette score could not be computed: {e}")
    
    return pca, pca_df


# Function to draw a convex hull around the points
def draw_convex_hull(points, color, radius = 0.75):
    if len(points) > 2:  # ConvexHull needs at least 3 points
        hull = ConvexHull(points)
        for simplex in hull.simplices:
            plt.plot(points[simplex, 0], points[simplex, 1], color=color, lw=2, alpha=0.2)
        plt.fill(points[hull.vertices, 0], points[hull.vertices, 1], color=color, alpha=0.1)
    elif len(points) == 1:  # Single point: draw a circle
        point = points[0]
        circle = Circle(point, radius, color=color, alpha=0.15)
        plt.gca().add_patch(circle)  # Add the circle to the current plot

def plotting_pca_results(pca_df, title):
    # Plot the PCA result with different colors for each cell type and model_id as labels
    fig = plt.figure(figsize=(6, 6))
    # Plot each cell type with labels for model_id
    for cell_type in pca_df['cell_type'].unique():
        subset = pca_df[pca_df['cell_type'] == cell_type]
        points = subset[['PC1', 'PC2']].values
        plt.scatter(points[:, 0], points[:, 1], color=colors[cell_type], label=cell_type)
        draw_convex_hull(points, colors[cell_type])
        
        # Add model_id labels to each point
        for i, model_id in enumerate(subset['model_id']):
            plt.text(points[i, 0], points[i, 1], model_id, fontsize=8, ha='right')
    
    # Finalize the plot
    plt.title(title)
    plt.xlabel('Principal Component 1')
    plt.ylabel('Principal Component 2')
    plt.grid(False)
    plt.legend()
    return fig

def plotting_pca_results_tot(pca_df, title):
    # Plot the PCA result with different colors for each cell type and model_id as labels
    fig = plt.figure(figsize=(15, 10))
    # Plot each cell type with labels for model_id
    for cell_type in pca_df['cell_type'].unique():
        subset = pca_df[pca_df['cell_type'] == cell_type]
        points = subset[['PC1', 'PC2']].values
        plt.scatter(points[:, 0], points[:, 1], color=colors_tot[cell_type], label=cell_type)
        draw_convex_hull(points, colors_tot[cell_type])
        
        # Add model_id labels to each point
        for i, model_id in enumerate(subset['model_id']):
            plt.text(points[i, 0], points[i, 1], model_id, fontsize=8, ha='right')
    
    # Finalize the plot
    plt.title(title)
    plt.xlabel('Principal Component 1')
    plt.ylabel('Principal Component 2')
    plt.grid(False)
    plt.legend()
    return fig

def plotting_pca_results_region(pca_df, title):
    # Plot the PCA result with different colors for each cell type and model_id as labels
    fig = plt.figure(figsize=(15, 10))
    # Plot each cell type with labels for model_id
    for cell_type in pca_df['cell_type'].unique():
        subset = pca_df[pca_df['cell_type'] == cell_type]
        points = subset[['PC1', 'PC2']].values
        plt.scatter(points[:, 0], points[:, 1], color=colors_region[cell_type], label=cell_type)
        draw_convex_hull(points, colors_region[cell_type])
        
        # Add model_id labels to each point
        for i, model_id in enumerate(subset['model_id']):
            plt.text(points[i, 0], points[i, 1], model_id, fontsize=8, ha='right')
    
    # Finalize the plot
    plt.title(title)
    plt.xlabel('Principal Component 1')
    plt.ylabel('Principal Component 2')
    plt.grid(False)
    plt.legend()
    return fig