import scanpy as sc

# Load annotated data
adata = sc.read_h5ad('results/basic/annotated_data.h5ad')

# Verify cell type annotations exist
print("Cell type columns:", [col for col in adata.obs.columns 
                            if col in ['leiden', 'cell_type', 'cluster']])

pip install heartmap[communication]

from heartmap import Config

# Load configuration
config = Config.default()

# Communication-specific settings
config.analysis.use_liana = True  # Use LIANA database
config.update_paths('./communication_analysis')
config.create_directories()

from heartmap.pipelines import AdvancedCommunicationPipeline

# Initialize pipeline
pipeline = AdvancedCommunicationPipeline(config)

# Run on annotated data
results = pipeline.run(
    data_path='results/basic/annotated_data.h5ad',
    output_dir='results/communication'
)

print("Communication analysis completed!")

import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns

# Extract communication scores
comm_scores = results['results']['communication_scores']

print(f"Found {len(comm_scores)} communication interactions")
print("\nTop 10 strongest interactions:")
top_interactions = comm_scores.nlargest(10, 'communication_score')
for _, row in top_interactions.iterrows():
    if 'ligand' in row:
        print(f"{row['source']} → {row['target']}: "
              f"{row['ligand']}-{row['receptor']} "
              f"(score: {row['communication_score']:.3f})")
    else:
        print(f"{row['source']} → {row['target']}: "
              f"{row['communication_score']:.3f}")

# Extract hub scores
hub_scores = results['results']['hub_scores']
adata = results['adata']

# Calculate mean hub score per cell type
cluster_col = 'leiden'  # or 'cell_type'
hub_by_type = adata.obs.groupby(cluster_col)['hub_score'].agg(['mean', 'std'])
hub_by_type = hub_by_type.sort_values('mean', ascending=False)

print("\nCommunication Hub Scores:")
for cell_type, scores in hub_by_type.head(10).iterrows():
    print(f"Cluster {cell_type}: {scores['mean']:.4f} ± {scores['std']:.4f}")

import scanpy as sc

# Load pre-generated heatmap
from pathlib import Path
from matplotlib.image import imread

fig_path = Path('results/communication/figures/communication_heatmap.png')
if fig_path.exists():
    img = imread(fig_path)
    plt.figure(figsize=(12, 10))
    plt.imshow(img)
    plt.axis('off')
    plt.title('Cell-Cell Communication Heatmap')
    plt.show()

# Visualize hub scores on UMAP
sc.pl.umap(adata, color='hub_score', 
           cmap='viridis',
           title='Communication Hub Scores',
           save='_hub_scores.png')

# Filter for specific ligand-receptor pairs
if 'ligand' in comm_scores.columns:
    # Example: VEGF signaling
    vegf_comm = comm_scores[
        comm_scores['ligand'].str.contains('VEGF', na=False) |
        comm_scores['receptor'].str.contains('FLT|KDR', na=False)
    ]
    
    print("\nVEGF Signaling:")
    for _, row in vegf_comm.iterrows():
        print(f"  {row['source']} → {row['target']}: "
              f"{row['ligand']}-{row['receptor']}")
    
    # Create pivot table for heatmap
    if not vegf_comm.empty:
        pivot = vegf_comm.pivot_table(
            values='communication_score',
            index='source',
            columns='target',
            aggfunc='sum',
            fill_value=0
        )
        
        plt.figure(figsize=(10, 8))
        sns.heatmap(pivot, annot=True, fmt='.2f', cmap='Reds')
        plt.title('VEGF Signaling Network')
        plt.tight_layout()
        plt.savefig('results/communication/vegf_network.png', dpi=300)

# Save communication scores
comm_scores.to_csv(
    'results/communication/communication_scores.csv',
    index=False
)

# Save hub scores
hub_by_type.to_csv(
    'results/communication/hub_scores_by_type.csv'
)

print("Communication data exported!")