Bio Applied Numerical Methods For Bioinformatics

Numerical Methods for Bioinformatics

When to Use

• Reconstructing missing or unmeasured time points in a qPCR/microarray time series
• Estimating a rate of change (e.g. dC/dt in a pharmacokinetic profile) from noisy discrete measurements
• Computing area-under-curve for dose-response curves, ROC curves, or PK exposure (AUC)
• Fitting nonlinear kinetics models (Michaelis-Menten, Hill equation) to enzyme/receptor-binding data
• Choosing a sampling/interpolation strategy that avoids Runge-phenomenon oscillation artifacts

Version Compatibility

NumPy >= 1.24, SciPy >= 1.11, scikit-learn >= 1.3 (for roc_curve), Python >= 3.10.

Prerequisites

pip install numpy scipy scikit-learn matplotlib

Basic calculus (derivatives, integrals) and familiarity with scipy.optimize.curve_fit are helpful. See bio-applied-enzyme-kinetics for a deeper dive into kinetics model fitting and bio-applied-statistics-for-bioinformatics for goodness-of-fit tests.

Interpolation

Goal: Reconstruct a value at an unmeasured x (e.g. a missing qPCR time point, ZT hour) from a small set of measured (x, y) pairs without introducing artificial oscillation.

Approach: Newton's divided differences give the same polynomial as Lagrange but let you add points incrementally with one extra term. High-degree polynomials on a uniform grid suffer the Runge phenomenon (oscillation near interval endpoints); Chebyshev nodes minimize that error, and scipy.interpolate.CubicSpline (piecewise cubic, continuous 1st/2nd derivatives) is the robust default for smooth biological signals such as circadian expression.

import numpy as np
from scipy.interpolate import CubicSpline


def newton_divided_differences(x, y):
    """Build the divided-difference table; return Newton coefficients c0..c_{n-1}."""
    x, y = np.asarray(x, float), np.asarray(y, float)
    n = len(x)
    table = np.zeros((n, n))
    table[:, 0] = y
    for j in range(1, n):
        for i in range(n - j):
            table[i, j] = (table[i + 1, j - 1] - table[i, j - 1]) / (x[i + j] - x[i])
    return table[0, :]


def newton_eval(coeffs, x_nodes, x_eval):
    """Evaluate the Newton interpolation polynomial via a Horner-like recurrence."""
    x_eval = np.asarray(x_eval, float)
    result = np.full_like(x_eval, coeffs[-1])
    for i in range(len(coeffs) - 2, -1, -1):
        result = result * (x_eval - x_nodes[i]) + coeffs[i]
    return result


def chebyshev_nodes(n, a=-1.0, b=1.0):
    """Return n Chebyshev nodes on [a, b] to avoid Runge-phenomenon oscillation."""
    k = np.arange(1, n + 1)
    return 0.5 * (b + a) + 0.5 * (b - a) * np.cos(np.pi * (2 * k - 1) / (2 * n))


# Gene expression: interpolate missing qPCR time points
t_measured = np.array([0, 2, 4, 8, 12, 24], dtype=float)  # hours
expr = np.array([1.0, 2.3, 4.1, 6.8, 5.2, 3.1])           # log2 fold change
coeffs = newton_divided_differences(t_measured, expr)
missing_t = np.array([6.0, 18.0])
expr_pred = newton_eval(coeffs, t_measured, missing_t)
print(f"Predicted: t=6h {expr_pred[0]:.2f}, t=18h {expr_pred[1]:.2f}")

# Circadian Per1 expression: cubic spline with periodic boundary condition
t_sampled = np.array([0, 3, 6, 9, 12, 15, 18, 21, 24], dtype=float)  # ZT hours
per1_expr = np.array([1.0, 3.2, 7.1, 9.4, 6.8, 2.9, 1.1, 0.8, 1.0])
cs = CubicSpline(t_sampled, per1_expr, bc_type="periodic")
t_fine = np.linspace(0, 24, 500)
peak_t = t_fine[np.argmax(cs(t_fine))]
print(f"Estimated peak at ZT {peak_t:.1f}h")

Numerical Differentiation and Integration

Goal: Estimate an instantaneous rate (e.g. drug clearance dC/dt) from noisy discrete measurements, and compute a stable area-under-curve (dose-response, ROC, PK exposure).

Approach: Differentiation amplifies noise (total error ~ M2*h + eps/h, minimized near h* ~ sqrt(eps/M2)) — use central differences (O(h^2)) and a moderate step, never h -> 0 on noisy real data. Integration is stable (errors average out): trapezoidal is O(h^2), Simpson's rule is O(h^4) but requires an even number of intervals.

import numpy as np


def central_diff(f, x, h=1e-3):
    """Central difference derivative estimate, O(h^2) accurate."""
    return (f(x + h) - f(x - h)) / (2 * h)


def trapezoid_rule(x, y):
    """Composite trapezoidal rule on a possibly non-uniform grid."""
    x, y = np.asarray(x, float), np.asarray(y, float)
    return np.sum(0.5 * (y[1:] + y[:-1]) * np.diff(x))


def simpsons_rule(x, y):
    """Composite Simpson's rule; requires a uniform grid with an even number of intervals."""
    x, y = np.asarray(x, float), np.asarray(y, float)
    n = len(x) - 1
    if n % 2 != 0:
        raise ValueError("Simpson's rule needs an even number of intervals")
    h = (x[-1] - x[0]) / n
    return h / 3 * (y[0] + y[-1] + 4 * np.sum(y[1:-1:2]) + 2 * np.sum(y[2:-2:2]))


# Drug clearance: dC/dt at interior points via central difference on measured points
t_pk = np.array([0, 1, 2, 4, 6, 8, 12, 24], dtype=float)
C_pk = np.array([100, 72, 52, 27, 14, 7.4, 2.0, 0.1])
for i in range(1, len(t_pk) - 1):
    dCdt = (C_pk[i + 1] - C_pk[i - 1]) / (t_pk[i + 1] - t_pk[i - 1])
    print(f"t={t_pk[i]:4.0f}h: dC/dt = {dCdt:6.2f} ng/mL/h")

# Dose-response AUC (% inhibition vs log10[drug])
dose_log = np.array([-9, -8, -7, -6, -5, -4, -3], dtype=float)
response = np.array([2.1, 3.4, 12.8, 48.3, 87.1, 97.2, 99.0])
auc = trapezoid_rule(dose_log, response)
print(f"AUC (dose-response): {auc:.2f} %*log-units")

Nonlinear Kinetics Curve Fitting

Goal: Fit an enzyme-kinetics or receptor-binding model to noisy data and pick the better model with an information criterion, not just R^2.

Approach: scipy.optimize.curve_fit (Levenberg-Marquardt) needs a sane initial guess (p0) — a bad guess converges to a local minimum or fails. Compare nested models (Michaelis-Menten vs Hill) with AIC; a lower AIC by >2 favors that model, >10 is decisive.

import numpy as np
from scipy.optimize import curve_fit


def michaelis_menten(S, Vmax, Km):
    """Michaelis-Menten kinetics: v = Vmax*S / (Km + S)."""
    return Vmax * S / (Km + S)


def hill_equation(S, Vmax, K_half, n):
    """Hill equation: v = Vmax*S^n / (K_half^n + S^n); n>1 indicates cooperativity."""
    return Vmax * S**n / (K_half**n + S**n)


def aic(n_params, y_obs, y_pred):
    """Akaike Information Criterion from residual sum of squares (Gaussian errors)."""
    n = len(y_obs)
    sse = np.sum((y_obs - y_pred) ** 2)
    return n * np.log(sse / n) + 2 * n_params


S_conc = np.array([0.05, 0.1, 0.2, 0.5, 1.0, 2.0, 5.0, 10.0])
v_obs = np.array([10.6, 18.5, 30.2, 51.8, 66.9, 77.5, 86.4, 90.1])

popt_mm, _ = curve_fit(michaelis_menten, S_conc, v_obs, p0=[np.max(v_obs), np.median(S_conc)])
popt_h, _ = curve_fit(hill_equation, S_conc, v_obs, p0=[np.max(v_obs), np.median(S_conc), 1.0], maxfev=5000)

aic_mm = aic(2, v_obs, michaelis_menten(S_conc, *popt_mm))
aic_h = aic(3, v_obs, hill_equation(S_conc, *popt_h))
print(f"MM:   Vmax={popt_mm[0]:.1f} Km={popt_mm[1]:.3f}  AIC={aic_mm:.2f}")
print(f"Hill: Vmax={popt_h[0]:.1f} K_half={popt_h[1]:.3f} n={popt_h[2]:.2f}  AIC={aic_h:.2f}")

Pitfalls

• Differentiation on noisy data: shrinking h does not improve accuracy — noise dominates below h* ~ sqrt(eps/M2); use central differences with a moderate step, or smooth (spline) first.
• Runge phenomenon: a single high-degree polynomial on a uniform grid oscillates wildly near endpoints — use Chebyshev nodes or piecewise splines instead of raising the polynomial degree.
• Simpson's rule requires an even number of intervals (odd number of points) on a uniform grid; passing an odd interval count raises/produces wrong results.
• curve_fit initial guess (p0): a poor p0 converges to a local minimum, hits maxfev, or silently returns a nonsensical fit — always plot the fit over the data and sanity-check parameter units.
• Coordinate systems: BED is 0-based half-open; VCF/GFF are 1-based inclusive — mixing them causes off-by-one errors when merging interpolated genomic coordinates.
• Multiple testing: when comparing many fitted curves/genes, apply FDR correction (Benjamini-Hochberg) rather than reading raw p-values.

See Also

• bio-applied-enzyme-kinetics — deeper Michaelis-Menten/Hill/allosteric model fitting and validation
• bio-applied-statistics-for-bioinformatics — goodness-of-fit, chi-squared, and model comparison tests
• bio-applied-bayesian-statistics-python — Bayesian alternatives to MLE/curve_fit point estimates
• bio-applied-machine-learning-for-biology — when a parametric curve no longer fits and you need a learned model