{"spec_id":"ma-differential-expression","library":"pygal","language":"python","code":"\"\"\" anyplot.ai\nma-differential-expression: MA Plot for Differential Expression\nLibrary: pygal 3.1.3 | Python 3.13.14\nQuality: 86/100 | Updated: 2026-06-21\n\"\"\"\n\nimport os\nimport sys\n\n\n# Script filename shadows the installed pygal package when run as `python pygal.py`;\n# dropping the script directory from sys.path lets the real package resolve.\nsys.path.pop(0)\n\nimport numpy as np\nimport pygal\nfrom pygal.style import Style\n\n\n# Theme-adaptive chrome tokens (Imprint style guide)\nTHEME = os.getenv(\"ANYPLOT_THEME\", \"light\")\nPAGE_BG = \"#FAF8F1\" if THEME == \"light\" else \"#1A1A17\"\nINK = \"#1A1A17\" if THEME == \"light\" else \"#F0EFE8\"\nINK_SOFT = \"#4A4A44\" if THEME == \"light\" else \"#B8B7B0\"\nINK_MUTED = \"#6B6A63\" if THEME == \"light\" else \"#A8A79F\"\n\n# Semantic exception: upregulated → Imprint green (#009E73, up/gain),\n# downregulated → Imprint matte red (#AE3030, loss/bad).\n# Non-significant uses INK_MUTED (theme-adaptive muted anchor).\nSERIES_COLORS = (\n    INK_MUTED,  # non-significant (muted background layer)\n    \"#009E73\",  # upregulated (Imprint green — semantic up/gain)\n    \"#AE3030\",  # downregulated (Imprint matte red — semantic loss/bad)\n    INK,  # M=0 reference line (neutral structural element)\n    INK_SOFT,  # +2-fold dashed threshold\n    INK_SOFT,  # -2-fold dashed threshold\n    \"#C475FD\",  # LOESS trend (Imprint lavender)\n    \"#4467A3\",  # top DE genes (Imprint blue)\n)\n\n# --- Data: Simulated RNA-seq differential expression results ---\nnp.random.seed(42)\nn_genes = 15000\n\n# Mean expression (A values) — log2 scale, typical RNA-seq range\nmean_expression = np.random.exponential(scale=3, size=n_genes) + 1\n\n# Log fold change (M values) — most genes near zero, some truly DE\nlog_fold_change = np.random.normal(0, 0.3, n_genes)\n\n# Add truly differentially expressed genes (~8% up, ~7% down)\nn_up = 1200\nn_down = 1050\nup_idx = np.random.choice(n_genes, n_up, replace=False)\nremaining = np.setdiff1d(np.arange(n_genes), up_idx)\ndown_idx = np.random.choice(remaining, n_down, replace=False)\n\nlog_fold_change[up_idx] = np.random.normal(2.5, 0.8, n_up)\nlog_fold_change[down_idx] = np.random.normal(-2.2, 0.7, n_down)\n\n# Simulate p-values (significant for DE genes, uniform noise otherwise)\np_values = np.ones(n_genes)\np_values[up_idx] = 10 ** (-np.random.uniform(2, 10, n_up))\np_values[down_idx] = 10 ** (-np.random.uniform(2, 10, n_down))\nnoise_idx = np.setdiff1d(np.arange(n_genes), np.concatenate([up_idx, down_idx]))\np_values[noise_idx] = np.random.uniform(0.01, 1.0, len(noise_idx))\n\nsignificant = p_values < 0.05\n\n# Notable gene names for the top most-significant hits\ngene_names = [f\"Gene{i}\" for i in range(n_genes)]\ntop_genes = [\"BRCA1\", \"TP53\", \"MYC\", \"EGFR\", \"KRAS\", \"PTEN\", \"CDK2\", \"RB1\", \"AKT1\", \"VEGFA\"]\ntop_idx = np.argsort(p_values)[: len(top_genes)]\nfor i, name in zip(top_idx, top_genes, strict=False):\n    gene_names[i] = name\n\n# LOESS-like smoothing curve (binned moving average)\nsort_order = np.argsort(mean_expression)\nsorted_a = mean_expression[sort_order]\nsorted_m = log_fold_change[sort_order]\n\nn_bins = 30\nbin_edges = np.percentile(sorted_a, np.linspace(0, 100, n_bins + 1))\nraw_x: list[float] = []\nraw_y: list[float] = []\nfor b in range(n_bins):\n    mask = (sorted_a >= bin_edges[b]) & (sorted_a < bin_edges[b + 1])\n    if mask.sum() > 20:\n        raw_x.append(float(np.median(sorted_a[mask])))\n        raw_y.append(float(np.mean(sorted_m[mask])))\n\nsmooth_y = np.array(raw_y)\nfor _ in range(4):\n    smoothed = np.copy(smooth_y)\n    for j in range(1, len(smoothed) - 1):\n        smoothed[j] = (smooth_y[j - 1] + smooth_y[j] + smooth_y[j + 1]) / 3\n    smooth_y = smoothed\nsmooth_x = raw_x\n\n# --- Subsample for SVG rendering performance ---\nnp.random.seed(42)\nsig_indices = np.where(significant)[0]\nnonsig_indices = np.where(~significant)[0]\nnonsig_sample = np.random.choice(nonsig_indices, min(1800, len(nonsig_indices)), replace=False)\n\nnonsig_points = []\nfor i in nonsig_sample:\n    nonsig_points.append(\n        {\n            \"value\": (round(float(mean_expression[i]), 2), round(float(log_fold_change[i]), 2)),\n            \"label\": f\"{gene_names[i]} | A={mean_expression[i]:.1f}, M={log_fold_change[i]:.2f}, p={p_values[i]:.2e}\",\n        }\n    )\n\ntop_idx_set = set(top_idx.tolist())\nsig_up_points = []\nsig_down_points = []\nfor i in sig_indices:\n    if i in top_idx_set:\n        continue\n    point = {\n        \"value\": (round(float(mean_expression[i]), 2), round(float(log_fold_change[i]), 2)),\n        \"label\": f\"{gene_names[i]} | A={mean_expression[i]:.1f}, M={log_fold_change[i]:.2f}, p={p_values[i]:.2e}\",\n    }\n    if log_fold_change[i] > 0:\n        sig_up_points.append(point)\n    else:\n        sig_down_points.append(point)\n\n# Top 10 most-significant DE genes — star marker in tooltip\nlabeled_points = []\nfor i in top_idx:\n    labeled_points.append(\n        {\n            \"value\": (round(float(mean_expression[i]), 2), round(float(log_fold_change[i]), 2)),\n            \"label\": f\"★ {gene_names[i]} | A={mean_expression[i]:.1f}, M={log_fold_change[i]:.2f}, p={p_values[i]:.2e}\",\n        }\n    )\n\nx_min = 0\nx_max = float(np.percentile(mean_expression, 99.5))\n\n# --- Style: Imprint palette + theme-adaptive chrome ---\ncustom_style = Style(\n    background=PAGE_BG,\n    plot_background=PAGE_BG,\n    foreground=INK,\n    foreground_strong=INK,\n    foreground_subtle=INK_SOFT,\n    colors=SERIES_COLORS,\n    title_font_size=66,\n    label_font_size=56,\n    major_label_font_size=44,\n    legend_font_size=44,\n    value_font_size=36,\n    stroke_width=2.5,\n    opacity=0.72,\n    opacity_hover=0.95,\n)\n\n# --- Chart (landscape 3200×1800 — canonical size) ---\nchart = pygal.XY(\n    width=3200,\n    height=1800,\n    style=custom_style,\n    title=\"ma-differential-expression · python · pygal · anyplot.ai\",\n    x_title=\"Mean Expression (A)\",\n    y_title=\"Log₂ Fold Change (M)\",\n    show_legend=True,\n    legend_at_bottom=True,\n    legend_at_bottom_columns=4,\n    legend_box_size=24,\n    dots_size=5,\n    stroke=False,\n    show_x_guides=True,\n    show_y_guides=True,\n    truncate_legend=-1,\n    print_values=False,\n    dynamic_print_values=True,\n    js=[],\n    x_label_rotation=0,\n    margin_bottom=110,\n)\n\n# Non-significant genes — muted background cloud\nchart.add(\"Not Significant\", nonsig_points, dots_size=3)\n\n# Upregulated significant — Imprint green; larger dots than downregulated for CVD redundancy\nchart.add(\"Upregulated (p<0.05)\", sig_up_points, dots_size=10)\n\n# Downregulated significant — Imprint matte red; smaller to distinguish from upregulated\nchart.add(\"Downregulated (p<0.05)\", sig_down_points, dots_size=5)\n\n# M = 0 reference line (no change)\nchart.add(\"M = 0\", [(x_min, 0), (x_max, 0)], stroke=True, show_dots=False, stroke_style={\"width\": 5})\n\n# M = +1 threshold (2-fold up)\nchart.add(\n    \"+2-fold\", [(x_min, 1), (x_max, 1)], stroke=True, show_dots=False, stroke_style={\"width\": 4, \"dasharray\": \"14, 8\"}\n)\n\n# M = -1 threshold (2-fold down)\nchart.add(\n    \"−2-fold\", [(x_min, -1), (x_max, -1)], stroke=True, show_dots=False, stroke_style={\"width\": 4, \"dasharray\": \"14, 8\"}\n)\n\n# LOESS smoothing curve — Imprint lavender; no dots so the curve stands out from scatter\nchart.add(\n    \"LOESS trend\",\n    [(round(x, 2), round(y, 3)) for x, y in zip(smooth_x, smooth_y, strict=False)],\n    stroke=True,\n    show_dots=False,\n    stroke_style={\"width\": 11},\n)\n\n# Top DE genes — Imprint blue, large prominent dots\nchart.add(\"Top DE genes\", labeled_points, dots_size=16, stroke=False)\n\n# --- Save ---\nwith open(f\"plot-{THEME}.html\", \"wb\") as f:\n    f.write(chart.render())\nchart.render_to_png(f\"plot-{THEME}.png\")\n"}