基于一阶导数同时置信带的基因表达显著变化时间窗参数化识别方法

《Biometrical Journal》：Identification of Changes in Gene Expression

　　

评估基因表达变化是许多研究领域的常见目标，例如在毒理学研究中，此类研究在临床前研究中尤为重要。在实践中，分析通常基于在实验观测时间点上进行的多重t检验，这限制了对基因表达变化发生的精确时间进行判断的准确性。若选择参数化途径，分析往往局限于识别效应的起始点，而非其持续长度。研究人员提出了一种参数化方法，以识别基因表达显著变化的时间窗。该方法通过拟合参数模型并为其一阶导数构建置信带（confidence band）实现。置信带通过两步自助法（two-step bootstrap）获得，并以假设检验的形式概括；拒绝原假设即表示检测到基因表达的显著变化。此外，还开发了一种为感兴趣时间点（例如显著变化的起始点）计算置信区间（confidence interval, CI）的方法。研究人员通过模拟研究验证了该方法的有效性，并将其多种不同应用展示于一项研究西方饮食对非酒精性脂肪性肝病（non-alcoholic fatty liver disease, NAFLD）进展影响的小鼠基因表达数据。

**论文解读**

## 1 研究背景与目的

RNA测序（RNA-seq）能够同时测量大量基因的表达，是毒理学与临床前研究中理解疾病反应机制的重要工具。常见分析目标之一是识别差异表达基因，例如剂量-基因表达数据或时间-基因表达数据。然而，常用的DESeq2等工具将剂量或时间作为离散协变量处理，只能比较已观测剂量水平或时间点，无法确定变化发生的精确时间。类似地，在时间-反应数据分析中，传统方法常通过多个t检验或Dunnett检验识别效应起始点，但无法刻画效应持续的长度。已有研究观察到一类“静止-跳跃基因”（rest-and-jump genes, RJG），即基因在延迟一段时间后才发生失调，因此需要关注整个显著变化时间窗。为此，研究人员提出一种参数化方法，通过拟合参数模型并构造一阶导数的置信带，识别基因表达显著变化的时间区间，从而支持对疾病进展和生物学过程的更完整推断。

## 2 研究方法与关键技术

研究人员首先对时间-基因表达数据拟合参数模型，如4参数逻辑斯蒂模型（4pLL）、beta模型、指数模型和二次模型，或拟合B样条（B-spline）作为更灵活的选择。模型参数通过最大似然法估计，并使用赤池信息准则（AIC）进行模型选择。研究重点不是模型本身，而是其关于时间的第一阶导数，因为一阶导数代表基因表达曲线的斜率，可直接衡量表达变化。进而构造一阶导数的下同时置信带，并建立假设检验：原假设为所有时间点上的一阶导数绝对值不超过预先指定的阈值ρ；若拒绝原假设，则表明存在显著表达变化。置信带通过两步自助法估计临界值。模拟研究的数据生成基于西方饮食小鼠试验中的真实基因参数，包括4pLL模型（基因Cd163）和beta模型（基因Fam83a），并设置六个情景和五个标准差水平。

案例研究数据来自一项西方饮食小鼠研究：共79只雄性小鼠分别喂食西方饮食（Western diet, WD）或标准饮食（standard diet, SD）最多48周，以评估高脂饮食对非酒精性脂肪性肝病（NAFLD）进展的影响。肝脏RNA经RNA-seq测序，使用Salmon、tximeta和DESeq2进行定量与预处理，并使用vst()函数标准化。

## 3 模拟研究结果

### 3.2.1 拒绝H₀

在1000次模拟运行中，情景1（无相关变化）在所有标准差水平下均未出现错误拒绝，I型错误率为0。对于由beta模型生成的情景4–6，方法在中等标准差以下均有很高功效；即使在大标准差下，功效仍分别约为81%（情景4）、95%（情景5）和78%（情景6）。对于4pLL模型生成的情景2和3，检测功效受标准差影响明显：小和中-小变异下较高，而在中-大和大变异下降至约30%–80%；情景3在中等标准差下仅检测到514次拒绝。研究还发现，在4pLL模型拐点附近，估计置信带出现明显凹陷，这在高变异数据中更显著。

### 3.2.2 估计的显著变化时间窗

对于单调的4pLL模型，若H₀被拒绝，通常识别出连续的时间窗；但在少数情况下，置信带凹陷会导致错误地识别出两个不连续子集。beta模型情景中，第二个显著变化时间窗的检测依赖于末期观测；情景4中较浅的下降仅在约120次（小标准差）或约40次（其他标准差）运行中被识别；情景5中更明显的下降在小、中、大标准差下分别被检测到999次、671次和260次。总体来看，显著变化起始点的估计呈正偏，终点估计呈负偏，导致时间窗整体略被低估；偏倚和方差随标准差增大而增大。

## 4 案例应用

### 4.1 应用于单个基因

研究人员将方法应用于WD小鼠中四个基因，每个基因分别用最合适的参数模型拟合：Cd163用4pLL模型，Fam83a用beta模型，Dbp用指数模型，Tm7sf2用二次模型。结果显示，Cd163的显著变化时间窗约从15.1周至20.3周（95%置信区间（CI））；Dbp从35.4周持续至研究结束（45周）；Fam83a从研究开始至24.9周；Tm7sf2从研究开始至20.3周。前两个基因符合典型的RJG形式，后两个则不是。基于B样条拟合的结果与参数模型大体一致，但Dbp被检测出两个显著变化时间窗，而Tm7sf2未检测到显著变化。增加B样条基函数次数和内部节点数会使置信带更波动，但核心活跃时间区域在不同设定下仍能被稳定检测。

### 4.2 GO富集分析

研究人员从Ghallab等研究的9881个差异表达基因中，使用4pLL模型筛选出569个在10至25周内发生显著表达变化的基因，并对这些基因进行基因本体（GO）富集分析。分析采用topGO软件包和Fisher精确检验。在最小校正p值的20个GO组中，包含“炎症反应”组，这与炎症是NAFLD进展中的关键事件一致。该应用说明，所提出方法可与GO富集分析结合，用于推断特定生物学过程发生的时间段。

## 5 讨论与结论

研究讨论指出，该方法存在时间窗估计偏短的问题，尤其在数据变异较大时更明显；B样条基函数的形式会影响置信带形状，过拟合可能使结果解释复杂化。未来可通过P样条惩罚、采用其他同时置信带构造方法或开发多重检验校正策略进行改进。总之，研究人员认为该方法是对现有基因表达数据分析工具的有价值补充，可在整个研究期间进行综合分析，并适用于多种表达曲线。