import matplotlib.pyplot as plt
import seaborn as sns
import numpy as np
import pandas as pd
from scipy import stats
# --- 1. Global Plot Settings ---
# Set general plot style for Nature-like aesthetics
sns.set_theme(
    context="paper",      # 'paper' for smaller elements, 'talk' for presentations
    style="white",        # No background grid by default
    palette="colorblind", # Excellent colorblind-friendly palette
    font_scale=1.2,       # Scale font sizes (adjust for desired text size in final output)
    rc={
        "font.family": "sans-serif",
        "font.sans-serif": ["Arial", "Helvetica", "DejaVu Sans"],
        "axes.edgecolor": "#333333", # Darker grey for axis lines
        "axes.labelcolor": "#333333",
        "xtick.color": "#333333",
        "ytick.color": "#333333",
        "figure.dpi": 300, # Default DPI for figures
        "savefig.dpi": 300 # Default DPI for saved figures
    }
)
# Define custom colors for conditions using the 'colorblind' palette for consistency
palette_custom = sns.color_palette("colorblind", 3)
colors = {
    'Vehicle': palette_custom[0],  # Blue
    'Drug X (Low)': palette_custom[1], # Orange
    'Drug X (High)': palette_custom[2] # Green
}
# --- 2. Simulate Data ---
np.random.seed(42) # for reproducibility
# Group sizes (number of replicates per condition)
n_samples = 15
# Data for Panel A (Cell Viability)
data_a = {
    'Condition': ['Vehicle'] * n_samples + ['Drug X (Low)'] * n_samples + ['Drug X (High)'] * n_samples,
    'Cell Viability (%)': np.concatenate([
        np.random.normal(100, 5, n_samples),   # Vehicle
        np.random.normal(115, 6, n_samples),   # Drug X (Low) - slight increase
        np.random.normal(130, 7, n_samples)    # Drug X (High) - significant increase
    ])
}
df_a = pd.DataFrame(data_a)
df_a['Cell Viability (%)'] = np.clip(df_a['Cell Viability (%)'], 80, 150) # Clip to reasonable range
# Data for Panel B (Gene A Expression)
data_b = {
    'Condition': ['Vehicle'] * n_samples + ['Drug X (Low)'] * n_samples + ['Drug X (High)'] * n_samples,
    'Gene A Expression (Arbitrary Units)': np.concatenate([
        np.random.normal(10, 1.5, n_samples),  # Vehicle
        np.random.normal(8, 1.2, n_samples),    # Drug X (Low) - slight decrease
        np.random.normal(5, 1.0, n_samples)     # Drug X (High) - significant decrease
    ])
}
df_b = pd.DataFrame(data_b)
df_b['Gene A Expression (Arbitrary Units)'] = np.clip(df_b['Gene A Expression (Arbitrary Units)'], 2, 12)
# Data for Panel C (Gene A vs. Gene B Correlation)
# Gene B is positively correlated with Gene A, but Drug X decreases Gene A
gene_a_values_c = df_b['Gene A Expression (Arbitrary Units)'].values
gene_b_values_c = 1.5 * gene_a_values_c + np.random.normal(5, 2, len(gene_a_values_c))
df_c = pd.DataFrame({
    'Gene A Expression (Arbitrary Units)': gene_a_values_c,
    'Gene B Expression (Arbitrary Units)': gene_b_values_c,
    'Condition': df_b['Condition'] # Keep conditions for potential coloring
})
# Data for Panel D (Heatmap of Differential Gene Expression)
gene_names = [f'Gene {i+1}' for i in range(15)]
conditions_heatmap = ['Vehicle', 'Drug X (Low)', 'Drug X (High)']
heatmap_data = pd.DataFrame(index=gene_names, columns=conditions_heatmap)
# Simulate different gene expression patterns
for i, gene in enumerate(gene_names):
    if i < 5: # Genes downregulated by Drug X
        heatmap_data.loc[gene, 'Vehicle'] = np.random.normal(10, 1)
        heatmap_data.loc[gene, 'Drug X (Low)'] = np.random.normal(8, 1)
        heatmap_data.loc[gene, 'Drug X (High)'] = np.random.normal(5, 0.8)
    elif i < 10: # Genes upregulated by Drug X
        heatmap_data.loc[gene, 'Vehicle'] = np.random.normal(5, 0.8)
        heatmap_data.loc[gene, 'Drug X (Low)'] = np.random.normal(7, 1)
        heatmap_data.loc[gene, 'Drug X (High)'] = np.random.normal(10, 1)
    else: # Genes largely unchanged
        heatmap_data.loc[gene, 'Vehicle'] = np.random.normal(7, 1)
        heatmap_data.loc[gene, 'Drug X (Low)'] = np.random.normal(7.5, 1)
        heatmap_data.loc[gene, 'Drug X (High)'] = np.random.normal(6.8, 0.9)
# Normalize heatmap data for visualization (e.g., Z-score scaling per gene)
heatmap_data_scaled = heatmap_data.apply(lambda x: (x - x.mean()) / x.std(), axis=1)
# --- 3. Create Multi-Panel Figure ---
fig, axes = plt.subplots(2, 2, figsize=(10, 9), constrained_layout=True) # 2 rows, 2 columns
# Adjust subplot spacing (optional, constrained_layout often handles well)
# plt.subplots_adjust(hspace=0.4, wspace=0.3)
# --- Panel A: Cell Viability (Violin Plot + Scatter) ---
ax = axes[0, 0]
sns.violinplot(
    data=df_a,
    x='Condition',
    y='Cell Viability (%)',
    ax=ax,
    inner='quartile', # Show quartiles inside violins
    linewidth=1,
    cut=0, # Trim violins to the data range
    palette=colors,
    saturation=0.7 # Reduce saturation slightly
)
sns.stripplot(
    data=df_a,
    x='Condition',
    y='Cell Viability (%)',
    ax=ax,
    jitter=0.2, # Add jitter to points to prevent overlap
    color='black',
    s=4, # Size of the individual points
    alpha=0.6
)
ax.set_title('A', loc='left', fontweight='bold')
ax.set_xlabel('') # Remove x-label as conditions are self-explanatory
ax.set_ylabel('Cell Viability (%)')
ax.set_ylim(70, 140) # Consistent y-axis limits
# Add statistical annotations (manual for demonstration, use `statannotations` for complex cases)
# Example: Drug X (High) vs Vehicle
y_max = df_a['Cell Viability (%)'].max()
ax.text(0, y_max + 5, '***', ha='center', va='bottom', fontsize=12) # Vehicle
ax.text(2, y_max + 5, '***', ha='center', va='bottom', fontsize=12) # Drug X (High)
ax.plot([0, 0, 2, 2], [y_max + 2, y_max + 4, y_max + 4, y_max + 2], lw=1, c='k')
ax.text(1, y_max + 8, '***', ha='center', va='bottom', fontsize=12) # Bracket for 0-2
ax.plot([0, 0, 1, 1], [y_max-5, y_max-2, y_max-2, y_max-5], lw=1, c='k') # Vehicle vs Drug X (Low)
ax.text(0.5, y_max-1, '*', ha='center', va='bottom', fontsize=12) # Bracket for 0-1
# --- Panel B: Gene A Expression (Violin Plot + Scatter) ---
ax = axes[0, 1]
sns.violinplot(
    data=df_b,
    x='Condition',
    y='Gene A Expression (Arbitrary Units)',
    ax=ax,
    inner='quartile',
    linewidth=1,
    cut=0,
    palette=colors,
    saturation=0.7
)
sns.stripplot(
    data=df_b,
    x='Condition',
    y='Gene A Expression (Arbitrary Units)',
    ax=ax,
    jitter=0.2,
    color='black',
    s=4,
    alpha=0.6
)
ax.set_title('B', loc='left', fontweight='bold')
ax.set_xlabel('')
ax.set_ylabel('Gene A Expression (AU)') # Use abbreviation
ax.set_ylim(2, 12)
# Add statistical annotations
y_min = df_b['Gene A Expression (Arbitrary Units)'].min()
ax.text(0, y_min - 2, '***', ha='center', va='top', fontsize=12) # Vehicle
ax.text(2, y_min - 2, '***', ha='center', va='top', fontsize=12) # Drug X (High)
ax.plot([0, 0, 2, 2], [y_min - 1, y_min - 1.5, y_min - 1.5, y_min - 1], lw=1, c='k')
ax.text(1, y_min - 2.5, '***', ha='center', va='top', fontsize=12)
# --- Panel C: Gene A vs. Gene B Correlation (Scatter Plot with Regression) ---
ax = axes[1, 0]
sns.scatterplot(
    data=df_c,
    x='Gene A Expression (Arbitrary Units)',
    y='Gene B Expression (Arbitrary Units)',
    hue='Condition', # Color points by condition
    palette=colors,
    s=50,
    alpha=0.7,
    ax=ax
)
# Add regression line for all data combined
slope, intercept, r_value, p_value, std_err = stats.linregress(
    df_c['Gene A Expression (Arbitrary Units)'],
    df_c['Gene B Expression (Arbitrary Units)']
)
sns.lineplot(
    x=df_c['Gene A Expression (Arbitrary Units)'],
    y=intercept + slope * df_c['Gene A Expression (Arbitrary Units)'],
    color='grey',
    linestyle='--',
    ax=ax
)
ax.text(
    0.05, 0.95,
    f'R={r_value:.2f}\nP={p_value:.3f}',
    transform=ax.transAxes,
    fontsize=10,
    verticalalignment='top',
    bbox=dict(boxstyle="round,pad=0.3", fc="white", ec="none", alpha=0.7)
)
ax.set_title('C', loc='left', fontweight='bold')
ax.set_xlabel('Gene A Expression (AU)')
ax.set_ylabel('Gene B Expression (AU)')
ax.legend(title='Condition', loc='upper left', bbox_to_anchor=(1, 1)) # Move legend outside
# Ensure legend does not overlap with text
ax.set_xlim(ax.get_xlim()[0], ax.get_xlim()[1] * 1.1) # Extend x-axis slightly if needed
ax.set_ylim(ax.get_ylim()[0], ax.get_ylim()[1] * 1.05) # Extend y-axis slightly if needed
# --- Panel D: Differential Gene Expression (Heatmap) ---
ax = axes[1, 1]
sns.heatmap(
    heatmap_data_scaled,
    cmap='RdBu_r', # Diverging colormap, red for high, blue for low
    center=0,      # Center the colormap at 0 (for scaled data)
    linewidths=0.5, # Lines between cells
    linecolor='white',
    cbar_kws={'label': 'Scaled Expression (Z-score)'}, # Color bar label
    yticklabels=True, # Show gene names
    xticklabels=True, # Show condition names
    ax=ax
)
ax.set_title('D', loc='left', fontweight='bold')
ax.set_xlabel('')
ax.tick_params(axis='y', rotation=0) # Ensure gene names are horizontal
ax.tick_params(axis='x', rotation=45, ha='right') # Rotate x-axis labels
ax.set_ylabel('Gene')
# --- 4. Final Touches & Save Figure ---
# Apply consistent formatting to all axes (remove top/right spines, inward ticks)
for ax in axes.flatten():
    ax.spines[['right', 'top']].set_visible(False)
    ax.tick_params(direction='in', length=4, width=1.5, colors='#333333') # Inward ticks
# Adjust layout to prevent overlaps
plt.tight_layout(rect=[0, 0, 0.98, 1]) # Adjust rect to make space for Panel C legend
# Save figure with high resolution
output_filename = 'nature_style_figure.png'
plt.savefig(output_filename, dpi=300, bbox_inches='tight')
plt.show()
print(f"Figure saved as '{output_filename}' with 300 DPI.")