用 Python 处理生物实验数据的一些经验

生物实验产生的数据量越来越大，用 Excel 手动处理不仅费时还容易出错——尤其在 IVD 试剂研发过程中，分析性能评估（参考 CMDE 2022年第32号《定量检测体外诊断试剂分析性能评估注册审查指导原则》）和临床试验数据处理（参考 2021年第72号《体外诊断试剂临床试验技术指导原则》）往往涉及大量重复性的计算和统计工作。Python 是实验室数据分析的利器，本文分享几个在日常工作中常用的场景和代码。

场景一：qPCR 结果的批量分析

qPCR 仪导出的表格通常包含几十甚至上百个孔的数据，包括 Ct 值、Tm 值等。用 Pandas 可以一键完成 ΔCt 和 ΔΔCt 的计算。这对于评估引物探针的灵敏度和特异性（CMDE 2024年第1号要求对主要原材料进行功能性验证）特别高效：

import pandas as pd
import numpy as np

df = pd.read_excel("qpcr_results.xlsx")
# 计算 ΔCt
df["delta_Ct"] = df["Ct_target"] - df["Ct_reference"]
# 对照组均值作为基线
baseline = df[df["group"]=="control"]["delta_Ct"].mean()
# ΔΔCt 和倍数变化
df["delta_delta_Ct"] = df["delta_Ct"] - baseline
df["fold_change"] = 2 ** (-df["delta_delta_Ct"])
# 按分组汇总
summary = df.groupby("group")[["Ct_target","fold_change"]].agg(["mean","std"])
print(summary)

场景二：标准曲线绘制与定量分析

做 ELISA 或定量 PCR，标准曲线是绕不开的。结合 scipy 做线性回归，计算 R² 和置信区间——这正是定量试剂分析性能评估中"线性范围"和"正确度"评价所需的核心计算：

import matplotlib.pyplot as plt
from scipy import stats
import numpy as np

# 标准品浓度和信号值
conc = np.array([0, 1.56, 3.125, 6.25, 12.5, 25, 50, 100])
signal = np.array([0.02, 0.15, 0.31, 0.58, 1.12, 2.15, 4.28, 8.52])

slope, intercept, r_value, p_value, std_err = stats.linregress(conc, signal)
r_squared = r_value ** 2

x_fit = np.linspace(0, 100, 100)
y_fit = slope * x_fit + intercept
# 95% 预测区间
y_pred = slope * conc + intercept
residuals = signal - y_pred
se = np.sqrt(np.sum(residuals**2) / (len(conc)-2))

fig, ax = plt.subplots(figsize=(6,4))
ax.scatter(conc, signal, c="#2563eb")
ax.plot(x_fit, y_fit, '--', color="#1e3a8a")
ax.fill_between(x_fit, y_fit-1.96*se, y_fit+1.96*se, alpha=0.12, color="#2563eb")
ax.set_xlabel("Concentration"); ax.set_ylabel("OD Value")
ax.text(0.05, 0.92, f"R² = {r_squared:.4f}", transform=ax.transAxes)
plt.tight_layout()
plt.savefig("standard_curve.png", dpi=150)

场景三：临床试验一致性分析

根据 CMDE 2021年第72号，IVD 临床试验需要计算阳性符合率、阴性符合率、总符合率及其 95% 置信区间，并与临床可接受标准（通常≥90%，参考 CMDE 2024年第4号 NGS 基因变异检测审查指导原则中 P₀≥90% 的标准）进行比较。手工做这些统计容易出错，用 Python 写一个标准化的分析函数一劳永逸：

from scipy.stats import norm
import numpy as np

def concordance_analysis(tp, tn, fp, fn):
    """计算诊断试剂与对比方法的一致性统计"""
    total = tp + tn + fp + fn
    ppa = tp / (tp + fn) if (tp + fn) > 0 else 0  # 阳性符合率
    npa = tn / (tn + fp) if (tn + fp) > 0 else 0  # 阴性符合率
    opa = (tp + tn) / total                     # 总符合率

    # Wilson Score 法计算 95% 置信区间
    def wilson_ci(p, n):
        z = norm.ppf(0.975)
        denom = 1 + z**2 / n
        center = (p + z**2/(2*n)) / denom
        margin = z * np.sqrt(p*(1-p)/n + z**2/(4*n**2)) / denom
        return max(0, center - margin), min(1, center + margin)

    ppa_ci = wilson_ci(ppa, tp + fn)
    npa_ci = wilson_ci(npa, tn + fp)

    print(f"阳性符合率(PPA): {ppa:.2%}  (95%CI: {ppa_ci[0]:.2%}-{ppa_ci[1]:.2%})")
    print(f"阴性符合率(NPA): {npa:.2%}  (95%CI: {npa_ci[0]:.2%}-{npa_ci[1]:.2%})")
    print(f"总符合率(OPA): {opa:.2%}")

# 示例: 申报试剂 vs Sanger测序金标准
concordance_analysis(tp=147, tn=93, fp=3, fn=5)

场景四：序列数据处理

引物探针设计是核酸检测试剂研发的核心环节。CMDE 2019年第83号（CYP2C19 检测试剂指导原则）和 2017年第52号（NIPT 指导原则）均对引物探针的靶向性和特异性验证提出了明确要求。Biopython 可以高效处理这些任务：

from Bio import SeqIO
from Bio.Seq import Seq

# 批量读取 FASTA 文件
for record in SeqIO.parse("sequences.fasta", "fasta"):
    seq = record.seq
    print(f">{record.id}")
    print(f"  长度: {len(seq)} bp")
    print(f"  GC含量: {(seq.count('G')+seq.count('C'))/len(seq)*100:.1f}%")
    print(f"  Tm值 (2+4法): {2*(seq.count('A')+seq.count('T')) + 4*(seq.count('G')+seq.count('C'))}°C")

# 翻译和酶切位点查找
coding_seq = Seq("ATGCGTAAGCTGTCGTCG")
protein = coding_seq.translate()
print(f"蛋白序列: {protein}")

场景五：数据可视化与报告生成

在注册申报中，分析性能数据需要以清晰、规范的图表呈现。Matplotlib + Seaborn 的组合可以满足大部分需求——从精密度评价的箱线图到方法学比对的 Bland-Altman 图（用于评估两种方法检测结果的一致性，是正确度评价的经典方法）：

import matplotlib.pyplot as plt
import numpy as np

def bland_altman(method_a, method_b):
    """Bland-Altman 一致性分析"""
    means = (np.array(method_a) + np.array(method_b)) / 2
    diffs = np.array(method_a) - np.array(method_b)
    bias = np.mean(diffs)
    sd = np.std(diffs)

    fig, ax = plt.subplots(figsize=(6,4))
    ax.scatter(means, diffs, alpha=0.6)
    ax.axhline(bias, color='#2563eb', linestyle='-')
    ax.axhline(bias + 1.96*sd, color='#94a3b8', linestyle='--')
    ax.axhline(bias - 1.96*sd, color='#94a3b8', linestyle='--')
    ax.set_xlabel("Mean of two methods")
    ax.set_ylabel("Difference (A - B)")
    ax.legend(["Bias", "±1.96 SD"])
    return fig

一个实用的工作习惯

所有数据分析脚本统一放在一个目录下，文件名用日期开头（如 2025-03-08_qpcr_analysis.py），方便以后查找。每个脚本开头写注释说明输入文件格式、依赖包版本和输出内容。关键计算步骤加上断言（assert）验证中间结果。这样半年后——或者注册发补时——回看自己写的代码能迅速理解上下文。

Python 在生物数据分析领域的生态已经非常成熟，入门门槛不高。无论你是要做试剂的性能评估统计分析，还是处理常规实验的批量数据，花一两个周末把基本功学会，回报是长期的效率提升。

注：本文中的统计分析方法和公式参考了 CMDE 发布的定量试剂分析性能评估指导原则（2022年第32号）和临床试验技术指导原则（2021年第72号），以及 CLSI-EP 系列标准。涉及注册申报的数据处理请以最新版法规和指导原则为准。