Statistical Analysis

BioLang includes a comprehensive statistics library designed for biological data. This tutorial covers hypothesis testing, correlation, dimensionality reduction, and clustering -- all the tools you need to analyze experimental results.

Step 1 — Descriptive Statistics

# A small experiment is embedded so every Run button works in the browser.
let data = table({
  sample: ["S01", "S02", "S03", "S04", "S05", "S06",
           "S07", "S08", "S09", "S10", "S11", "S12"],
  group: ["control", "control", "control", "control", "control", "control",
          "treated", "treated", "treated", "treated", "treated", "treated"],
  treatment: ["A", "A", "A", "A", "B", "B", "B", "B", "C", "C", "C", "C"],
  expression: [10.2, 9.8, 10.5, 11.0, 12.4, 12.8, 13.1, 13.5, 15.0, 15.4, 15.8, 16.1],
  before_treatment: [8.2, 8.0, 8.5, 8.7, 9.1, 9.0, 9.4, 9.2, 9.8, 9.7, 10.0, 10.1],
  after_treatment:  [8.8, 8.6, 9.1, 9.0, 9.8, 9.7, 10.2, 10.0, 10.6, 10.5, 10.9, 11.0],
  gene_a_expression: [4.1, 4.4, 4.8, 5.0, 5.4, 5.7, 6.0, 6.3, 6.7, 7.0, 7.4, 7.7],
  gene_b_expression: [3.8, 4.2, 4.5, 4.9, 5.1, 5.5, 5.9, 6.1, 6.5, 6.9, 7.1, 7.6],
  gene_a: [4.1, 4.4, 4.8, 5.0, 5.4, 5.7, 6.0, 6.3, 6.7, 7.0, 7.4, 7.7],
  gene_b: [3.8, 4.2, 4.5, 4.9, 5.1, 5.5, 5.9, 6.1, 6.5, 6.9, 7.1, 7.6],
  gene_c: [7.2, 7.0, 6.8, 6.5, 6.2, 6.0, 5.8, 5.5, 5.2, 5.0, 4.7, 4.5],
  gene_d: [2.0, 2.5, 2.1, 2.8, 2.4, 3.0, 2.7, 3.2, 2.9, 3.5, 3.1, 3.8],
  gene_e: [5.0, 4.9, 5.1, 5.0, 5.2, 5.1, 5.3, 5.2, 5.4, 5.3, 5.5, 5.4],
  dosage: [0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 5.5],
})

# Extract a column as a list (table field access returns column)
let values = col(data, "expression")

println("=== Descriptive Statistics ===")
println(f"N:        {len(values)}")
println(f"Mean:     {round(mean(values), 3)}")
println(f"Median:   {round(median(values), 3)}")
println(f"Std Dev:  {round(stdev(values), 3)}")
println(f"Variance: {round(variance(values), 3)}")
println(f"Min:      {min(values)}")
println(f"Max:      {max(values)}")
println(f"Range:    {max(values) - min(values)}")

# Quantiles
let q25 = quantile(values, 0.25)
let q50 = quantile(values, 0.50)
let q75 = quantile(values, 0.75)
println(f"\n25th percentile: {round(q25, 3)}")
println(f"50th percentile: {round(q50, 3)}")
println(f"75th percentile: {round(q75, 3)}")
println(f"IQR:             {round(q75 - q25, 3)}")

# Summary for all numeric columns at once
let summary = describe(data)
println(summary)

Step 2 — Comparing Distributions

# Extract groups for comparison
data |> filter(|r| r.group == "control") |> into ctrl_rows
data |> filter(|r| r.group == "treated") |> into treat_rows
let control = col(ctrl_rows, "expression")
let treated = col(treat_rows, "expression")

# Kolmogorov-Smirnov test: do two samples come from the same distribution?
let ks = ks_test(control, treated)
println("=== KS Test (control vs treated) ===")
println(f"Statistic: {round(ks.statistic, 4)}")
println(f"p-value:   {round(ks.pvalue, 4)}")

if ks.pvalue < 0.05 {
  println("The samples show a detectable distribution difference.")
} else {
  println("No clear distribution difference was detected; this does not prove identical or normal data.")
}

# Compare basic distribution properties
println(f"\nControl — mean: {round(mean(control), 3)}, stdev: {round(stdev(control), 3)}")
println(f"Treated — mean: {round(mean(treated), 3)}, stdev: {round(stdev(treated), 3)}")

Step 3 — t-Tests and Non-Parametric Alternatives

State the design and method rather than asking software to diagnose them from the values. These are independent groups, so Welch's test is a useful default when a difference in means answers the question.

import "statistics" as stat

# Independent groups; unequal variances are allowed.
let t_result = stat.compare_groups(control, treated, {method: "welch"})
println("=== Welch two-sample t-test ===")
println(stat.explain_task(t_result))
println(f"mean difference: {round(t_result.analysis.mean_diff, 4)}")
println(f"p-value:         {round(t_result.analysis.p_value, 6)}")
println(f"df:              {round(t_result.analysis.df, 1)}")

# Paired t-test (for before/after measurements)
let before = col(data, "before_treatment")
let after  = col(data, "after_treatment")
let paired = stat.paired_change(before, after, {method: "paired_t"})
println(f"\nPaired t-test p-value: {round(paired.analysis.p_value, 6)}")

# Non-parametric: Wilcoxon rank-sum test (two independent samples)
let w = stat.compare_groups(control, treated, {method: "mann_whitney"})
println("\n=== Wilcoxon Rank-Sum Test ===")
println(f"U-statistic: {round(w.analysis.u_statistic, 1)}")
println(f"p-value:     {round(w.analysis.p_value, 6)}")

Mann-Whitney is a rank/distribution comparison; it is not automatically a test of medians. Pairing, independence, and the experimental unit come from the study design, not the histogram.

Step 4 — ANOVA

# One-way ANOVA for multiple groups
# group_by returns a Map of group_name -> Table
let groups = group_by(data, "treatment")

# Extract expression columns from each group into a list of lists
let group_names = colnames(data) |> filter(|c| c == "treatment")
let treatment_names = col(data, "treatment") |> unique()

let group_lists = treatment_names |> map(|name| {
  let subset = data |> filter(|r| r.treatment == name)
  col(subset, "expression")
})

let anova_result = anova(group_lists)
println("=== One-Way ANOVA ===")
println(f"F-statistic: {round(anova_result.f_statistic, 3)}")
println(f"p-value:     {round(anova_result.p_value, 6)}")
println(f"df between:  {anova_result.df_between}")
println(f"df within:   {anova_result.df_within}")

# If significant, do pairwise t-tests with p-value adjustment
if anova_result.p_value < 0.05 {
  println("\n=== Pairwise Comparisons ===")
  let pvals = []
  let labels = []
  for i in 0..len(treatment_names) {
    for j in (i+1)..len(treatment_names) {
      let a = group_lists[i]
      let b = group_lists[j]
      let t = ttest(a, b)
      pvals = pvals + [t.p_value]
      labels = labels + [f"{treatment_names[i]} vs {treatment_names[j]}"]
    }
  }
  let adjusted = p_adjust(pvals, "bh")
  for k in 0..len(labels) {
    let sig = if adjusted[k] < 0.05 { "*" } else { "ns" }
    println(f"  {labels[k]}: p={round(adjusted[k], 4)} {sig}")
  }
}

Step 5 — Correlation

# Pearson correlation (returns a Float)
let x = col(data, "gene_a_expression")
let y = col(data, "gene_b_expression")

let r = cor(x, y)
println("=== Pearson Correlation ===")
println(f"r: {round(r, 4)}")

# Spearman rank correlation (returns a Record)
let rho = spearman(x, y)
println(f"\nSpearman rho: {round(rho.coefficient, 4)}")
println(f"Spearman p:   {round(rho.pvalue, 4)}")

# Kendall rank correlation (returns a Record)
let tau = kendall(x, y)
println(f"Kendall tau:  {round(tau.coefficient, 4)}")
println(f"Kendall p:    {round(tau.pvalue, 4)}")

# Compute correlations for multiple gene pairs
let gene_cols = ["gene_a", "gene_b", "gene_c", "gene_d", "gene_e"]
println("\n=== Pairwise Correlations ===")
for i in 0..len(gene_cols) {
  for j in (i+1)..len(gene_cols) {
    let a = col(data, gene_cols[i])
    let b = col(data, gene_cols[j])
    let r = cor(a, b)
    println(f"  {gene_cols[i]} ~ {gene_cols[j]}: r={round(r, 3)}")
  }
}

Step 6 — Linear Regression

# Simple linear regression: lm(x_list, y_list)
let dosage = col(data, "dosage")
let expression = col(data, "expression")
let model = lm(dosage, expression)

println("=== Linear Regression ===")
println(f"Intercept:  {round(model.intercept, 3)}")
println(f"Slope:      {round(model.slope, 3)}")
println(f"R-squared:  {round(model.r_squared, 4)}")
println(f"p-value:    {round(model.p_value, 6)}")

# Visualize the fit
plot(data, {
  theme:  "publication",
  x:      "dosage",
  y:      "expression",
  type:   "scatter",
  title:  "Expression vs Dosage",
})

Step 7 — PCA Visualization

# Select an expression matrix from the embedded experiment.
let expr = select(data, "sample", "group", "gene_a", "gene_b", "gene_c", "gene_d", "gene_e")

# pca_plot performs PCA and renders a scatter plot of PC1 vs PC2
# It accepts a table and a config record
pca_plot(expr, {
  theme:     "publication",
  group_col: "group",
  labels:    true,
  title:     "PCA — RNA-seq Samples",
})

Step 8 — Clustering and Heatmaps

# Clustered heatmap with hierarchical clustering built in
let all_cols = colnames(expr)
let heatmap_cols = all_cols |> filter(|c| c != "sample" and c != "group")
let matrix = select(expr, ...heatmap_cols)

clustered_heatmap(matrix, {theme: "publication"})

Step 9 — Multiple Testing and FDR

# When running many tests, we need to correct for multiple comparisons

# Run a t-test per gene
let all_cols = colnames(expr)
let gene_cols = all_cols |> filter(|c| c != "sample" and c != "group")

let results = gene_cols |> map(|gene| {
  let ctrl_rows = expr |> filter(|r| r.group == "control")
  let treat_rows = expr |> filter(|r| r.group == "treated")
  let ctrl = col(ctrl_rows, gene)
  let treat = col(treat_rows, gene)
  let result = ttest(ctrl, treat)
  { gene: gene, pvalue: result.p_value, t_stat: result.t_statistic }
}) |> to_table()

# Collect raw p-values and apply correction methods
let raw_pvals = col(results, "pvalue")

let bh_adjusted   = p_adjust(raw_pvals, "bh")
let bonf_adjusted  = p_adjust(raw_pvals, "bonferroni")
let holm_adjusted  = p_adjust(raw_pvals, "holm")

# Build the corrected table with an explicit index. Table rows do not expose
# a hidden `_index` field.
let result_genes = col(results, "gene")
let result_pvalues = col(results, "pvalue")
let result_tstats = col(results, "t_stat")
let corrected = range(0, nrow(results)) |> map(|i| {
  { gene: result_genes[i], pvalue: result_pvalues[i], t_stat: result_tstats[i],
    bh_fdr: bh_adjusted[i], bonferroni: bonf_adjusted[i], holm: holm_adjusted[i] }
}) |> to_table()

# Compare methods
for method in ["bh_fdr", "bonferroni", "holm"] {
  let n_sig = corrected |> filter(|r| r[method] < 0.05) |> nrow()
  println(f"{method}: {n_sig} significant genes")
}

# Write results
corrected
  |> arrange("bh_fdr")
  |> write_csv("results/gene_tests.csv")

Step 10 — Complete Statistical Workflow

# stats_pipeline.bl — full statistical analysis

fn statistical_analysis(data_file, output_dir) {
  let data = read_csv(data_file)
  mkdir(output_dir)

  # 1. Descriptive statistics
  let desc = describe(data)
  desc |> write_csv(f"{output_dir}/descriptive_stats.csv")

  # 2. Summary per group
  let treatment_names = col(data, "treatment") |> unique()
  for name in treatment_names {
    let subset = data |> filter(|r| r.treatment == name)
    let vals = col(subset, "expression")
    println(f"{name}: mean={round(mean(vals), 3)}, stdev={round(stdev(vals), 3)}, n={len(vals)}")
  }

  # 3. ANOVA across groups
  let group_lists = treatment_names |> map(|name| {
    let subset = data |> filter(|r| r.treatment == name)
    col(subset, "expression")
  })
  let anova_result = anova(group_lists)
  println(f"\nANOVA: F={round(anova_result.f_statistic, 3)}, p={round(anova_result.p_value, 6)}")

  # 4. PCA visualization
  let all_cols = colnames(data)
  let numeric_cols = all_cols |> filter(|c| c != "sample" and c != "treatment" and c != "group")
  let pca_input = select(data, ...(numeric_cols + ["treatment"]))
  pca_plot(pca_input, {group_col: "treatment", theme: "publication"})
    |> save_svg(f"{output_dir}/pca.svg")

  # 5. Clustered heatmap
  clustered_heatmap(select(data, ...numeric_cols), {theme: "publication"})
    |> save_svg(f"{output_dir}/heatmap.svg")

  println(f"\nAnalysis complete. Results in {output_dir}/")
}

statistical_analysis("data/experiment.csv", "results")