酵母转录调控网络微阵列数据分析全流程：从原始数据清洗到因果表达网络线性建模（google-research yeast_transcription_network）

酵母转录调控网络微阵列数据分析全流程：从原始数据清洗到因果表达网络线性建模（google-research yeast_transcription_network）

【免费下载链接】google-researchGoogle Research项目地址: https://gitcode.com/gh_mirrors/go/google-research

本篇技术指南以 google-research 仓库 yeast_transcription_network 目录为对象，系统讲解配套论文《Time-resolved genome-scale profiling reveals a causal expression network》的完整数据处理与建模管线：从原始双色微阵列信号出发，经过斑点聚合、异常值修剪、t0 归一化、串扰修复、时间趋势清洗与噪声模型构建，最终得到cleaned/thresholded数据集并用于线性建模与因果网络预测。读完本文，你将掌握如何下载并组织公开数据文件、按顺序运行 YeastNetwork.ipynb 中的各阶段代码、复现论文中的三种数据集（raw / cleaned / thresholded），以及如何用 sklearn 构建设计矩阵、拟合 LARS 模型并把完整预测拆解为每个调控系数（cause）的独立贡献。

一、项目背景与数据文件概览

该目录的核心是一份名为YeastNetwork.ipynb的 Python Notebook（Google Colab 格式，nbformat 4），它是上述论文的配套实现。Notebook 中从原始微阵列读数一直做到因果表达网络的线性模型构建，覆盖面完整，是研究酵母基因调控网络（转录因子 TF 扰动后的时间序列表达）的可复现实战代码。

1.1 配套数据文件

论文发布时配套了四个数据文件（README 中逐一说明），需要从论文作者公开的下载渠道获取后本地使用：

文件名	内容说明
yeast_raw_data_table_20180826.tsv	"原始"微阵列数据（包含红/绿双通道信号值）
yeast_data_table_20180826.tsv	"处理过"的微阵列数据（清洗、阈值化之后的最终表）
insample_coefs_20170601.csv	验证实验之前的预测模型系数
insample_coefs_20180826.csv	加入验证实验之后重训的预测模型系数

从源码可以确认每个文件的实际用法：

• 原始数据表由Process "raw" data阶段的 读取与验证代码 使用，用于从rProcessedSignal/gProcessedSignal两列重新计算r_g_ratio并验证与表内该列一致；
• 处理数据表由Read processed data阶段读取，并用于与本地重建结果做差验证；
• 两个系数文件在Prediction model阶段被pd.read_csv载入，其中insample_coefs_20170601.csv是默认使用的"验证实验前"版本（Notebook 中另一行被注释掉，切换即可使用 20180826 版本）。

1.2 运行依赖

Notebook 的导入代码（首代码单元）说明运行所需的最小依赖集合：

import os import pandas as pd import numpy as np import scipy.stats import sklearn.linear_model

即pandas（表格与 MultiIndex 操作）、numpy（矩阵运算）、scipy.stats（正态分布分位数，用于噪声标定）、scikit-learn（线性建模示例）。

二、论文中定义的三种数据集

README 明确给出了论文使用的三个数据集与数据表中列的对应关系，这是理解整条管线的关键索引：

数据集	对应数据表列
raw	r_g_ratio（位于yeast_raw_data_table_20180826.tsv）
cleaned	log2_cleaned_ratio（位于yeast_data_table_20180826.tsv）
thresholded	log2_cleaned_ratio_zth2dfilt（位于yeast_data_table_20180826.tsv）

从 Notebook 源码可以更精确地还原这些列的来历：

• r_g_ratio是红通道（突变株）与绿通道（对照）处理信号之比，且在2.0处做了下限截断：np.maximum(rProcessedSignal, 2.0) / np.maximum(gProcessedSignal, 2.0)（见 raw 数据验证单元）；
• log2_cleaned_ratio由log_cleaned_ratio块除以log(2.0)转换而来（biologist_table函数），log_cleaned_ratio则是对数比经过 t0 归一化、串扰修复与时间趋势清洗后的结果；
• log2_cleaned_ratio_zth2dfilt是"z 化 + 二维阈值 + 时间序列过滤"后的版本，其中z表示按时间序列零值处理（把不显著的时间过程剔除）、th2d表示使用"基因 × 芯片"二维噪声模型做阈值化、filt表示激进过滤（剔除只有单点检出或存在检测断点的序列）。

三、运行 Notebook 的完整步骤

3.1 准备数据目录

将四个数据文件下载到本地后，编辑 Notebook 中如下代码单元，把占位路径替换为实际位置：

# Replace this with the location of the downloaded data and models datadir = '/path/to/datafiles'

所有后续的读取代码（yeast_raw_data_table_20180826.tsv、yeast_data_table_20180826.tsv、insample_coefs_*.csv）都是通过os.path.join(datadir, ...)定位文件的，因此该变量是唯一需要修改的配置点。

3.2 按顺序执行各阶段

README 给出了官方推荐的执行顺序：

1. 运行Common functions中的代码单元：加载全部辅助函数与常量定义，这是后续所有阶段的前置条件；
2. （可选）运行Process the "raw" data：从原始读数重建处理数据，用于验证处理流程可被完整复现；
3. 运行Read processed data：读取官方发布的数据表并进入建模阶段；只有在前一步Process the "raw" data确实运行过时，才需要（也才有意义）执行其中的Verify the processing procedure验证单元。

需要特别强调的是，Process the "raw" data与Read processed data是二选一的两条入口：前者从零重建（耗时、占内存），后者直接读官方结果（快速）。两者的验证逻辑是：把本地重建的df_bio与读入的df_bio_read逐列做差，检查min(v), max(v)是否为 0（见 验证单元）。

3.3 验证实验过滤（复现第一轮结果的关键开关）

在build_cleaned_table(df_dm)调用之前，Notebook 提供了一个"时间过滤"开关，用于复现验证实验之前的原始训练集：

validation_f = pd.to_datetime(df_dm.index.get_level_values('date')) > '2017-07-01' # 复现验证实验之前的训练集（清洗与噪声模型输出会略有不同）： # df_dm = build_cleaned_table(df_dm[~validation_f]) # 复现完整数据集： df_dm = build_cleaned_table(df_dm)

由于数据清洗与噪声模型依赖全量数据，而 2017-07-01 之后新增的验证实验轮次会影响统计量，README 与源码注释均提醒：若要严格复现论文第一轮结果，应启用被注释掉的df_dm[~validation_f]过滤分支。

四、数据清洗管线源码深度解析

Common functions定义了整条管线的核心算子。以下按数据流顺序逐个拆解。

4.1 差分算子：construct_difference_operators

该函数根据时间序列（若干条从 0 开始的严格递增序列拼接而成）构造三个len(times) × len(times)矩阵算子（源码）：

• x0：作用在表达向量上，返回"每个时间过程 t=0 点取值复制到整条序列"的向量；
• diff：有限差分算子，对时间过程内部使用对称差分（相邻时间点ddt = 0.5 / (shifted_times - times)），对实验边界处置 0；
• diff_inv：差分算子的"伪逆"——由于差分算子不可逆，实现上删去每条序列的第一列与最后一行后求逆，满足diff_inv . diff . X == X - X0。

这三个算子是后续所有归一化、求导与"积分还原"操作的数学基础：x0用于 t0 归一化与截距校正，diff用于把表达水平转化为瞬时变化率（即d/dt ln(gene)），diff_inv则用于从变化率还原表达曲线（预测模型阶段）。

函数同时包含严格的输入校验：首时间点必须为 0，否则抛出ValueError('First time point is non-zero, exiting.')；时间序列必须严格递增排序，否则抛出ValueError('Timeseries not sorted properly.')。

4.2 设计表索引体系

Notebook 定义了四级 MultiIndex 常量（源码），理解它们等于理解了数据的组织方式：

TIMECOURSE_IDX = ['TF', 'strain', 'date', 'restriction', 'mechanism'] CHIP_IDX = TIMECOURSE_IDX + ['time'] TIMESERIES_IDX = ['GeneName', 'SystematicName'] + CHIP_IDX DESIGN_IDX = CHIP_IDX + ['seq', 'rseq']

• TIMECOURSE_IDX：标识一条独立的时间过程（转录因子、菌株、实验日期、限制性内切酶、机制 GEV/ZEV）；
• CHIP_IDX：在时间过程之上叠加time，标识单张微阵列芯片；
• TIMESERIES_IDX：进一步叠加基因名，标识每个"基因 × 芯片"数据点；
• DESIGN_IDX：在芯片基础上追加seq（时间点在该时间过程内的序号）与rseq（倒序序号），供差分/过滤逻辑使用。

4.3 从原始信号到设计表

Process the "raw" data阶段的数据整形流程为：

1. 读取并验证：读取 TSV，验证r_g_ratio列可由双通道信号在 2.0 截断后相除重建；
2. 跨斑点聚合：用pivot_table以TIMESERIES_IDX为索引，对r_g_ratio、rProcessedSignal、gProcessedSignal分别计算len、median、min、max、std（微阵列通常每个基因有 2 个重复斑点，取中位数聚合，见 聚合单元）；
3. 异常值修剪：trim_timeseries_outliers用两类规则修正离群点（源码）：
   • 斑点不一致（MAX_SPOT_RATIO = 4.0）：同一基因在芯片上多个斑点amax/amin > 4时，用前后时间点的几何平均（端点用邻居值）替换中位数；
   • 时间趋势尖峰（JUMP_THRESHOLD = 4.0）：非端点位置出现"先上后下或先下后上且每步超过 4 倍"的跳变时，同样以相邻点几何平均替换；
4. 构建设计表：design_table以CHIP_IDX为行索引、GeneName为列，抽取r_g_median、red_median、green_median三个块；随后用timeseries_sequence计算seq/rseq列，并按DESIGN_IDX排序；最后dropna(axis=1)剔除 GEV 与 ZEV 芯片转录本集合的差异基因，收敛到公共基因集（见 注释单元）。

4.4 清洗三步曲：归一化、串扰修复、时间趋势

build_cleaned_table把清洗流程串成固定顺序（源码）：

def build_cleaned_table(df_dm): x0_op, _, _ = construct_difference_operators(...) # 仅需要 x0 df_dm = normalize_to_t0(df_dm, x0_op) # 1. 相对 t0 归一化 df_dm = repair_crosstalk_outliers(df_dm, x0_op) # 2. 红绿串扰修复 df_dm = clean_time_trends(df_dm) # 3. 中位时间趋势清洗 df_dm = construct_noise_model(df_dm) # 4. 构建噪声模型 return df_dm

• normalize_to_t0：以基因表达中位数r_g_median为基准，构造log_ratio = log(r_g_median / (x0_op · r_g_median))，即每个时间点的对数比相对于该时间过程 t=0 的比值；
• repair_crosstalk_outliers：针对"绿通道信号暴涨（>8 倍）且未校正对数比同步上涨（>log 2）"的基因，把log(green/green0)加回log_ratio，从而把串扰造成的假信号纠正为以 t=0 绿通道为基准的表达式（修复时打印Repairing %s, %s timecourse for gene %s.）；
• clean_time_trends：按TIME_BINS时间窗（(0,6), (6,12), (12,17.5), (17.5,25), (25,35), (35,85), (85,1000)，对应典型时间过程 0/6/12/17.5/25/35/85 分钟取样）分别对 GEV / 非 GEV 样本减去窗内中位数，消除系统性的时间批次效应，输出log_cleaned_ratio块。

4.5 噪声模型与阈值化

噪声模型（construct_noise_model）是一个"基因 × 芯片"的二维误差矩阵，核心思想是把误差分解为基因层面与芯片层面（源码）：

• 基因级误差：log_cleaned_gene_error用中央五分位（60% 分位 − 40% 分位）乘以CENTRAL_QUINTILE_SCALE（正态分布下约等于 1/0.506，即1/(norm.ppf(0.6) − norm.ppf(0.4))）换算成标准差估计；
• 阈值因子：threshold_factor = sqrt(2 * log(n_eff))，n_eff为time > 0的样本数（来自极大值统计：N 次正态抽样最大值的期望水平）；
• 芯片级误差：先对log_cleaned_ratio做硬阈值（小于max(gene_error × factor, MINIMUM_THRESHOLDING=0.1)的置 0），标记显著时间过程并剥离，对剩余"非显著"数据逐芯片计算中央五分位误差，t=0 点芯片误差置 0；
• 合并：noise = sqrt(gene_error² + CHIP_WEIGHT × chip_error²)，其中CHIP_WEIGHT = 0.2；非零时间点再除以sqrt(1 + CHIP_WEIGHT)做归一化，最终输出log_noise_model块。

阈值化（apply_all_thresholding）在噪声模型之上生成 9 个后缀版本（源码），对应LOG_SUFFIXES中cleaned_ratio_zth/hth/th/zth2d/hth2d/th2d/zth2dfilt/hth2dfilt/th2dfilt的排列：

• 硬阈值（hth）：|x| < noise的项置 0；
• 软阈值（th）：median(x+noise, x−noise, 0)，即经典的软收缩；
• 一维阈值：使用 t=0 基因噪声广播（threshold_error）；二维阈值（2d）：使用完整的"基因 × 芯片"噪声矩阵（threshold_error_2d）；
• 序列过滤（filt，激进模式）：额外剔除"只有一次检出且不在末点"或"检出存在断点"的劣质时间过程（apply_hard_and_soft_thresholding(..., aggressive=True)）。

阈值化过程中会打印检出统计量，例如timecourses with signal、timecourses with one detection、timecourses with only final detection、timecourses with no detection gaps，便于核对数据质量。

4.6 输出 log2 生物学家格式表

biologist_table把设计表"堆叠"（stack）为长表，并把所有log_前缀列除以log(2.0)转为 log2 单位，列名改为log2_noise_model / log2_ratio / log2_cleaned_ratio / log2_cleaned_ratio_*等，同时保留green_median / red_median / r_g_median原始通道列（源码）。这正是发布文件yeast_data_table_20180826.tsv的列结构来源。

五、线性建模：从设计矩阵到 LARS 拟合

5.1 设计矩阵构建

construct_full_design是建模核心（源码），其要点：

• 筛选目标基因自身的过表达实验：对基因gene建模时，剔除TF == gene的实验（"自我促进实验"，避免自调控混淆）；
• 列块设计：u_m为全部基因的表达矩阵（列块val = 'cleaned_ratio' + thresholded），u为目标基因列，二次项uu = u * u_m，设计矩阵d_m = hstack([u_m, uu])，即一阶 + 二次（双线性）的调控模型；
• 稳态参考：d_m -= 1.0使截距表示"t=0 偏离稳态"；delta_model=True时假设稳态（截距为 0）；
• 对数模型：log_weights=True时设计矩阵除以目标基因u，因变量切换为log_cleaned_ratio，对应拟合d/dt ln(gene)；
• 差分/积分模式：inverse_design=False时用diff_op对因变量求差分（并过滤每条序列的末点，因末点导数线性相关），权重weights = 1/(diff_op**2).sum(axis=1)；inverse_design=True时则用inv_diff_op对设计矩阵积分、过滤 t=0 行；
• 返回(设计矩阵, 因变量, t0 基值, 权重, 行索引, 系数标签, 样本内过滤, 样本外过滤)。

5.2 以 FKH1 为例的 sklearn 拟合

Notebook 以转录因子FKH1作为建模示例（示例代码）：

gene = 'FKH1' promoted = df_dm.index.get_level_values('TF') == gene in_sample_f = np.full(len(df_dm.index), True, dtype=bool) out_of_sample_f = np.full(len(df_dm.index), False, dtype=bool) in_sample_f = in_sample_f[~promoted] df = df_dm[~promoted] times = np.array(df.index.get_level_values('time')) x0_op, diff_op, inv_diff_op = construct_difference_operators(times) (design_m, y_dep, y_base, weights, row_index, col_index, in_sample_f, out_of_sample_f) = construct_full_design( df, in_sample_f, out_of_sample_f, gene, delta_model=True, # 假设稳态，截距为 0 thresholded='_zth2dfilt', # 使用过滤后的二维阈值数据 model_quadratic_terms=True, log_weights=True, # 拟合 d/dt ln(gene) inverse_design=False, ...)

拟合使用 sklearn 的 LARS（最小角回归）：

larsmodel = sklearn.linear_model.Lars(fit_intercept=False, normalize=False) larsmodel.fit(design_m, y_dep) f = larsmodel.coef_path_[:, lambdacol] != 0 # lambdacol = 40，取正则化路径第 40 步 coefs = (larsmodel.coef_path_)[f, lambdacol] for gene, coef in zip(col_index[f], coefs): print('%15s % .4e' % (gene, coef))

即沿 LARS 正则化路径选取第lambdacol=40步的非零系数作为调控关系集合。README 特别说明：论文原文使用的是一众与 sklearn 不完全兼容的 Fortranglmnet封装之一，Notebook 中的 sklearn 示例用于展示建模流程，而非逐位复现论文数值。

六、预测模型：把完整预测拆解为单系数贡献

6.1 读取系数并构造预测

Prediction model阶段从系数文件恢复调控网络，并以YHB1基因为例演示预测（示例代码）：

coeffs = pd.read_csv(os.path.join(datadir, 'insample_coefs_20170601.csv')) design = df_dm['cleaned_ratio_zth2dfilt'] x0op, diff, diff_inv = construct_difference_operators(np.array(design.index.get_level_values('time'))) gene = 'YHB1' tf_filt = design.index.get_level_values('TF') != gene # 不预测自身过表达实验中的该基因 actualg = design[gene] actualg_diff = np.matmul(diff, np.log(actualg)) # 实际 d/dt ln(gene) f = coeffs['effect'] == gene # 影响该基因的系数行

对每条影响gene的系数（cause列指定调控因子，coef列指定系数值）：

if cause.startswith('quad'): marginal = (design[cause[5:]] - 1.0 / design[gene]) * coef else: marginal = (design[cause] - 1.0) / design[gene] * coef

• 非二次项：(design[cause] − 1) / design[gene] × coef，对应d/dt ln(gene)中"因子与目标之比"贡献；
• 二次项：(design[cause] − 1/gene) × coef；
• 选定模型强制截距为 0（delta_model稳态假设）。

所有边际贡献求和得到prediction，再用diff_inv算子"积分"并取指数还原为表达水平：

data_pred = np.exp(np.matmul(diff_inv, prediction))

Notebook 同时构造marginal_diff_df（差分域边际）与marginal_df（积分还原后的表达域边际），把完整预测逐列拆分为['data', 'prediction'] + causes，实现"完整预测 = 各调控系数贡献之和"的可视化分解。

6.2 单个时间过程调查

最后，示例以AFT1转录因子的实验为例筛选边际贡献表，并保留贡献超过阈值的列：

aft1 = marginal_df.query('TF=="AFT1"') hits = np.sum(np.abs(np.log(aft1)) > 0.2) > 0 aft1[aft1.columns[hits]]

即仅展示在 AFT1 时间过程中对数表达变化超过0.2（约 22% 变化）的调控因子列，用于快速定位该时间过程中的主导调控因子。

七、复现注意事项与适用前提

结合 README 与 Notebook 源码，复现或改造时需注意以下几点：

1. 数据与模型文件必须齐备：datadir下需同时具备 2 个 TSV 与 2 个 CSV；若只跑Read processed data之后的阶段，TSV 是必需的，CSV 在Prediction model阶段必需；
2. 两条入口不可混淆：直接读官方处理表时不要运行Verify the processing procedure，否则拿本地重建与官方表比对属于无效验证；而重建路径下该验证是检验管线正确性的唯一手段（差值min/max应为 0 量级）；
3. GEV/ZEV 芯片基因集差异：dropna(axis=1)会把两套芯片共有的转录本之外的基因剔除，因此最终公共基因集小于两套芯片的并集；
4. 时间序列格式约束：所有时间过程必须以 0 起始且严格递增，construct_difference_operators会在违反时抛异常；
5. 模型数值与论文的差异：论文使用 Fortranglmnet封装，Notebook 的 sklearnLars仅用于流程演示，权重方案无法完全等价，重现论文精确数值需使用论文原始工具链；
6. 验证实验开关：若目标是复现论文第一轮（验证实验之前）的结果，需启用df_dm[~validation_f]过滤（以2017-07-01为界），此时清洗与噪声模型的统计量与全量版本略有差异。

八、小结

yeast_transcription_network提供了一条从原始双色微阵列信号到因果调控网络模型的端到端可复现管线：差分算子支撑的数学框架（construct_difference_operators）贯穿归一化、求导与积分还原；基因 × 芯片二维噪声模型配合软/硬阈值与时间序列过滤，产出论文定义的cleaned与thresholded数据集；construct_full_design与 LARS 拟合完成因果系数学习；Prediction model把网络预测还原为每个调控因子的边际贡献。对于希望复现论文结果、学习微阵列时间序列清洗方法或构建基因调控网络线性模型的读者，这份 Notebook 是一份结构清晰、可直接运行的完整参照实现。

【免费下载链接】google-researchGoogle Research项目地址: https://gitcode.com/gh_mirrors/go/google-research