Statistics for Bioinformatics

A two-group test is only meaningful after you decide what counts as an independent sample. This guide uses one numeric feature measured in independent control and treated biological samples. It walks through the input contract, a Welch comparison, its uncertainty, and the separate problem of testing many features.

Start with the experimental unit

Write down the unit that was independently assigned to a condition or independently sampled from the biological population. If treatment was applied to animals, the animal may be the unit; if it was applied to cultures, the independently prepared culture may be. Multiple reads, wells, fields of view, or repeated instrument measurements from one unit do not become independent biological replicates merely because they occupy separate rows. Summarize technical repeats at the biological-unit level or use a model that represents their nesting.

For a single feature, the comparison below targets the difference in group means, defined as treated minus control. It assumes independent observations within and between groups. If the same biological units are measured in both conditions, use a paired analysis with correctly matched units; the independent-samples code here is inappropriate. If samples share a litter, donor, batch, or other dependence, consider a design-aware model rather than treating each row as independent.

Input contract and reproducible setup

Choose either the Python or base R walkthrough below; each runs the same synthetic inputs independently and produces the same core Welch estimate, interval, statistic, degrees of freedom, and p-value. Both expect two CSV files, each with a numeric column named value, one row per independent biological unit, and no duplicated technical measurements. Genuinely blank cells are discarded and counted; nonblank values that fail numeric conversion, infinities, a missing column, fewer than two retained values, or zero variance stop with an error. The output directory is created if needed. The included values are synthetic teaching fixtures, not biological evidence; replace them only with measurements whose experimental-unit structure you have checked.

# Run from a disposable working directory. Do not overwrite real data.
mkdir -p data results
printf "value\n8.9\n9.8\n10.4\n9.5\n" > data/control.csv
printf "value\n11.2\n10.7\n12.1\n11.4\n" > data/treated.csv

Equivalent test in base R

R's built-in t.test uses Welch's two-sample method by default. Save the displayed code as compare_groups.R and run it with the same two CSV files. This branch uses base R only and writes the estimate, interval, statistic, degrees of freedom, and p-value.

#!/usr/bin/env Rscript
args <- commandArgs(trailingOnly = TRUE)
if (length(args) != 3) stop("Usage: Rscript compare_groups.R CONTROL.csv TREATED.csv OUTDIR")
read_values <- function(path) {
  d <- read.csv(path, check.names = FALSE, colClasses = "character",
                na.strings = character(), blank.lines.skip = FALSE)
  if (!("value" %in% names(d))) stop(paste(path, "must contain a value column"))
  raw <- trimws(as.character(d$value))
  blank <- is.na(d$value) | raw == ""
  x <- suppressWarnings(as.numeric(raw[!blank]))
  if (anyNA(x)) stop(paste(path, "contains a nonblank value that is not numeric"))
  cat(path, ": discarded ", sum(blank), " blank value(s)\n", sep = "")
  if (length(x) < 2 || any(!is.finite(x))) stop(paste(path, "needs at least two finite numeric values"))
  if (var(x) == 0) stop(paste(path, "sample variance is zero; inspect these data"))
  x
}
ctrl <- read_values(args[1])
trt <- read_values(args[2])
dir.create(args[3], recursive = TRUE, showWarnings = FALSE)
res <- t.test(trt, ctrl, var.equal = FALSE, conf.level = 0.95)
write.csv(data.frame(
  n_control = length(ctrl), n_treated = length(trt),
  mean_difference_treated_minus_control = mean(trt) - mean(ctrl),
  ci95_low = res$conf.int[1], ci95_high = res$conf.int[2],
  statistic = unname(res$statistic), df = unname(res$parameter), p_value = res$p.value
), file.path(args[3], "comparison_R.csv"), row.names = FALSE)

Rscript --version
Rscript compare_groups.R data/control.csv data/treated.csv results

For Python, use Python 3 and install NumPy, pandas, and SciPy 1.11 or newer in an environment you manage; save the displayed code as compare_groups.py. No R packages are needed for the base R route. The Python script tests one feature. Its Welch confidence interval is for the treated-minus-control mean difference; Cohen's d is a standardized descriptive effect estimate, not a substitute for the difference in the feature's original units.

python -m pip install numpy pandas 'scipy>=1.11'
# After saving the displayed Python code as compare_groups.py:
python compare_groups.py data/control.csv data/treated.csv results

#!/usr/bin/env python3
import sys
from pathlib import Path
import numpy as np
import pandas as pd
from scipy import stats

if len(sys.argv) != 4:
    raise SystemExit("Usage: python compare_groups.py CONTROL.csv TREATED.csv OUTDIR")
control_path, treated_path, outdir = map(Path, sys.argv[1:])

def read_values(path):
    frame = pd.read_csv(path, dtype=str, keep_default_na=False, skip_blank_lines=False)
    if "value" not in frame.columns:
        raise ValueError(f"{path}: required column 'value' is missing")
    raw = frame["value"]
    blank = raw.isna() | raw.astype("string").str.strip().eq("")
    try:
        converted = pd.to_numeric(raw[~blank], errors="raise")
    except (ValueError, TypeError) as exc:
        raise ValueError(f"{path}: nonblank value is not numeric") from exc
    values = converted.to_numpy(dtype=float)
    print(f"{path}: discarded {int(blank.sum())} blank value(s)")
    if not np.isfinite(values).all():
        raise ValueError(f"{path}: values must be finite")
    if len(values) < 2:
        raise ValueError(f"{path}: need at least two non-missing units")
    if np.var(values, ddof=1) == 0:
        raise ValueError(f"{path}: sample variance is zero; inspect these data")
    return values

ctrl = read_values(control_path)
trt = read_values(treated_path)
outdir.mkdir(parents=True, exist_ok=True)

# Welch's test allows the two population variances to differ.
res = stats.ttest_ind(trt, ctrl, equal_var=False)
ci = res.confidence_interval(confidence_level=0.95)
mean_diff = float(np.mean(trt) - np.mean(ctrl))
pooled_sd = np.sqrt(((len(ctrl)-1)*np.var(ctrl, ddof=1) +
                     (len(trt)-1)*np.var(trt, ddof=1)) /
                    (len(ctrl)+len(trt)-2))
cohens_d = mean_diff / pooled_sd

pd.DataFrame([{
    "n_control": len(ctrl), "n_treated": len(trt),
    "mean_control": np.mean(ctrl), "mean_treated": np.mean(trt),
    "sd_control": np.std(ctrl, ddof=1), "sd_treated": np.std(trt, ddof=1),
    "mean_difference_treated_minus_control": mean_diff,
    "welch_t": res.statistic, "welch_df": res.df, "p_value": res.pvalue,
    "ci95_low": ci.low, "ci95_high": ci.high, "cohens_d_pooled_sd": cohens_d
}]).to_csv(outdir / "comparison.csv", index=False)
print(f"Wrote {outdir / 'comparison.csv'}")

Inspect the data and interpret the comparison

• Plot the raw observations by biological unit, alongside group summaries. With few samples, histograms and formal normality tests have little power to establish that a distribution is normal. Look for data-entry problems, severe skew, and influential observations, then check them against records instead of deleting them because they change a p-value.
• Welch's two-sample t-test targets equality of population means for independent groups and does not assume equal population variances. It still relies on a suitable sampling/design structure and a mean-based model; very small samples, strong skew, or influential outliers can make inference fragile. No single diagnostic certifies assumptions.
• The p-value describes how extreme a statistic at least this large would be under the null model and test assumptions. It is not the probability that the null is true, the probability the result is “due to chance,” or a measure of biological importance.
• Read the mean difference and interval before the p-value. The interval shows values compatible with the data and model at its stated confidence level; it is not a probability distribution for the fixed true effect. Decide in advance what difference would matter biologically, using domain knowledge rather than generic small/medium/large labels.
• A non-significant result does not establish no effect. It may reflect imprecision, limited independent replication, a small effect, or model mismatch. Increasing technical measurements alone does not replace additional independent biological units.

When thousands of features are tested

The example above produces one raw p-value. If a separate, appropriately specified analysis generates one valid p-value per feature, adjust the prespecified family together rather than interpreting thousands of unadjusted values as independent discoveries. The Benjamini-Hochberg procedure controls the expected false-discovery proportion under its conditions; dependence structure and analysis design matter, and the adjustment does not fix invalid p-values, confounding, or inappropriate feature filtering.

Choose either language below. Each adjusts the same four synthetic p-values in their original order and writes an adjusted-value file plus a 0/1 rejection file. For real analyses, obtain p-values from a design-appropriate model for every feature, document the tested family and filtering rules, then adjust that complete vector. These four numbers are only a software demonstration.

First, create the shared input below in a Bash terminal, not an R or Python console. This step applies to both routes. Then save and run only the script for your chosen language; you do not need the other language's installation or run commands.

printf "0.01\n0.04\n0.20\n0.80\n" > pvalues_all_features.txt

Base R with p.adjust

This equivalent route uses the base R p.adjust function with method BH; it needs no add-on packages. Save as adjust_features.R and run it in the same disposable working directory as the shared input above; skip the Python branch entirely.

#!/usr/bin/env Rscript
lines <- trimws(readLines("pvalues_all_features.txt", warn = FALSE))
if (length(lines) == 0 || any(lines == "")) {
  stop("Expected one nonblank p-value per line")
}
pvals <- suppressWarnings(as.numeric(lines))
if (anyNA(pvals)) stop("Every input line must be numeric")
if (any(!is.finite(pvals)) || any(pvals < 0 | pvals > 1)) {
  stop("P-values must be finite and between 0 and 1")
}
alpha <- 0.05
p_adjusted <- p.adjust(pvals, method = "BH")
write.table(p_adjusted, "pvalues_adjusted.txt", row.names = FALSE,
            col.names = FALSE, quote = FALSE)
write.table(as.integer(p_adjusted <= alpha), "discoveries_fdr05.txt",
            row.names = FALSE, col.names = FALSE, quote = FALSE)
cat("adjusted:", paste(signif(p_adjusted, 8), collapse = " "), "\n")
cat("reject:", as.integer(p_adjusted <= alpha), "\n")

Rscript --version
Rscript adjust_features.R

Python with statsmodels

This optional Python branch requires Python 3, NumPy, and statsmodels. Save the code as adjust_features.py. After saving it, run these commands in the same disposable working directory as the shared input.

python -m pip install numpy statsmodels
python adjust_features.py

#!/usr/bin/env python3
from pathlib import Path
import numpy as np
from statsmodels.stats.multitest import multipletests

lines = Path("pvalues_all_features.txt").read_text(encoding="utf-8").splitlines()
if not lines or any(not line.strip() for line in lines):
    raise ValueError("Expected one nonblank p-value per line")
try:
    pvals = np.asarray([float(line.strip()) for line in lines], dtype=float)
except ValueError as exc:
    raise ValueError("Every input line must be numeric") from exc
if not np.isfinite(pvals).all() or ((pvals < 0) | (pvals > 1)).any():
    raise ValueError("P-values must be finite and between 0 and 1")
reject, p_adjusted, _, _ = multipletests(pvals, alpha=0.05, method="fdr_bh")
np.savetxt("pvalues_adjusted.txt", p_adjusted)
np.savetxt("discoveries_fdr05.txt", reject.astype(int), fmt="%d")
print("adjusted:", p_adjusted)
print("reject:", reject.astype(int))

For this fixture, BH-adjusted values are [0.04, 0.08, 0.2666667, 0.8] (approximately), and rejection flags at alpha 0.05 are [1, 0, 0, 0]. Adjusted values quantify multiplicity-adjusted evidence; a rejection is not proof of a biological effect. Both scripts preserve input order, so keep feature IDs in a parallel file and verify row alignment before joining results. Report the correction procedure, alpha, number and definition of tests, raw and adjusted p-values, effect estimates, and uncertainty where available. Non-numeric, non-finite, or out-of-range p-values are errors, not values to silently discard.

Choose a method that matches the design

• Paired measurements: analyze within-unit differences with a paired method, or fit a model that represents the repeated structure. Do not pair groups by row order.
• Count-based sequencing measurements: use methods designed for count distributions and their mean–variance structure, with the experimental design represented. A t-test on raw counts is not the default differential-expression workflow.
• Strong skew or influential values: inspect the measurement process and estimand. A transformation changes the scale and often the question; rank-based tests do not automatically test equality of means. Choose an alternative based on the scientific target and design.
• Batch, donor, litter, or covariates: include them in a model only when the design supports estimating their effects. If condition is completely confounded with batch, statistical adjustment cannot separate the two.
• Power and planning: define the smallest meaningful effect, outcome, design, and expected variability before choosing sample size. A post-hoc “observed power” calculation does not resolve an imprecise result; use a prospective design-specific calculation and document its assumptions.

Record an analysis plan

Before looking at group p-values, write a short plan naming the outcome/feature, comparison and direction, experimental unit, inclusion rules, independent or paired structure, covariates or blocks, primary model, diagnostic checks, multiplicity family and correction, and the effect size that would matter biologically. In the results, report the exact independent sample counts, missing-data handling, estimate with interval, test and settings, raw/adjusted p-values as appropriate, software versions, and limitations. Preserve reference-genome and annotation release details when feature mapping depends on them.

Use the checklist to decide whether this simple comparison fits. If you cannot identify independent units, if samples are paired or clustered, or if the feature is a sequencing count, stop before running it and choose a design-appropriate method.

References

• SciPy documentation: independent-samples t-test
• statsmodels documentation: multiple-testing corrections
• R documentation: t.test
• Benjamini Y, Hochberg Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. Journal of the Royal Statistical Society Series B. 1995.