Example 3: Clustering Analysis

import pandas as pd
import numpy as np
from sklearn.preprocessing import StandardScaler
from sklearn.cluster import KMeans
import matplotlib.pyplot as plt
import seaborn as sns

# Set random seed for reproducibility
np.random.seed(42)

# Read the TSV file
df = pd.read_csv('Filtered_Diff_Exp.tsv', sep='\t')

# Let's look at the first few rows and basic information about our dataset
print("First few rows of the dataset:")
print(df.head())
print("\nDataset information:")
print(df.info())

# Calculate mean expression for wild type (WT) and hy5 mutant
df['WT_mean'] = df[['Ath_WT_R1_expression', 'Ath_WT_R2_expression']].mean(axis=1)
df['hy5_mean'] = df[['Ath_hy5_R1_expression', 'Ath_hy5_R2_expression']].mean(axis=1)

# Create a new dataframe with just the mean values
expression_data = df[['feature_ids', 'WT_mean', 'hy5_mean']]

print("First few rows with mean expression values:")
print(expression_data.head())

# Extract features for clustering (excluding gene IDs)
X = expression_data[['WT_mean', 'hy5_mean']].values

# Scale the data
scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)

print("Shape of scaled data:", X_scaled.shape)
print("\nFirst few rows of scaled data:")
print(X_scaled[:5])

# Calculate inertia for different numbers of clusters
inertias = []
K = range(1, 11)

for k in K:
    kmeans = KMeans(n_clusters=k, random_state=42)
    kmeans.fit(X_scaled)
    inertias.append(kmeans.inertia_)

# Plot elbow curve
plt.figure(figsize=(10, 6))
plt.plot(K, inertias, 'bx-')
plt.xlabel('k')
plt.ylabel('Inertia')
plt.title('Elbow Method For Optimal k')
plt.grid(True)
plt.show()

# Create a copy of the data first
kmeans = KMeans(n_clusters=4, random_state=42)
clusters = kmeans.fit_predict(X_scaled)         
expression_data = expression_data.copy()
expression_data['Cluster'] = clusters
print("Number of genes in each cluster:")
print(expression_data['Cluster'].value_counts())

# Create scatter plot
plt.figure(figsize=(12, 8))
scatter = plt.scatter(expression_data['WT_mean'], 
                     expression_data['hy5_mean'],
                     c=expression_data['Cluster'], 
                     cmap='viridis')
plt.xlabel('Wild Type Expression (mean)')
plt.ylabel('hy5 Mutant Expression (mean)')
plt.title('Gene Expression Clusters')
plt.colorbar(scatter, label='Cluster')

# Add a diagonal line for reference
min_val = min(expression_data['WT_mean'].min(), expression_data['hy5_mean'].min())
max_val = max(expression_data['WT_mean'].max(), expression_data['hy5_mean'].max())
plt.plot([min_val, max_val], [min_val, max_val], 'r--', alpha=0.5)

plt.grid(True, alpha=0.3)
plt.show()

def calculate_cluster_patterns(df, fold_change_threshold=1.5):
    """
    Analyze gene expression patterns for each cluster, providing comprehensive statistics
    about how genes respond to the absence of HY5.
    
    Args:
        df: DataFrame containing gene expression data with WT_mean and hy5_mean columns
        fold_change_threshold: Minimum fold change to consider a gene significantly regulated
        
    Returns:
        Three DataFrames containing basic statistics, expression ratios, and regulation patterns
    """
    # Calculate basic statistics for each cluster
    # We separate WT and hy5 calculations for clarity
    wt_stats = df.groupby('Cluster')['WT_mean'].agg(['mean', 'std', 'count']).round(2)
    hy5_stats = df.groupby('Cluster')['hy5_mean'].agg(['mean', 'std']).round(2)
    
    # Combine statistics into a multi-level column DataFrame
    cluster_stats = pd.concat({
        'WT_mean': wt_stats,
        'hy5_mean': hy5_stats
    }, axis=1)
    
    # Calculate expression ratios with careful handling of low expression values
    def calculate_stable_ratio(group):
        """
        Calculate expression ratios while avoiding division by very small numbers.
        Only considers genes with meaningful WT expression (>0.1) to prevent inflated ratios.
        """
        meaningful_expression = group['WT_mean'] > 0.1
        if meaningful_expression.any():
            ratios = group[meaningful_expression]['hy5_mean'] / group[meaningful_expression]['WT_mean']
            return ratios.mean()
        return 0  # Return 0 if no genes meet the criteria
    
    # Calculate ratios using our stable method
    ratios = df.groupby('Cluster').apply(
        calculate_stable_ratio,
        include_groups=False
    ).round(2)
    
    # Calculate the percentage of genes showing significant regulation
    patterns = pd.DataFrame({
        'Upregulated in hy5': df.groupby('Cluster').apply(
            lambda x: (x['hy5_mean'] > x['WT_mean'] * fold_change_threshold).mean() * 100,
            include_groups=False
        ),
        'Downregulated in hy5': df.groupby('Cluster').apply(
            lambda x: (x['hy5_mean'] < x['WT_mean'] / fold_change_threshold).mean() * 100,
            include_groups=False
        )
    }).round(1)
    
    return cluster_stats, ratios, patterns

# Calculate all statistics
cluster_stats, cluster_ratios, regulation_patterns = calculate_cluster_patterns(expression_data)

# Print results with informative headers
print("Cluster Statistics:")
print("This shows the mean expression, standard deviation, and number of genes in each cluster")
print(cluster_stats)
print("\nAverage hy5/WT expression ratio for each cluster:")
print("(Ratios calculated only for genes with meaningful WT expression)")
print(cluster_ratios)
print("\nRegulation patterns (percentage of genes):")
print("(Shows what fraction of genes are significantly up- or down-regulated)")
print(regulation_patterns)