酵母转录调控网络微阵列数据分析全流程从原始数据清洗到因果表达网络线性建模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 NotebookGoogle 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 明确给出了论文使用的三个数据集与数据表中列的对应关系这是理解整条管线的关键索引数据集对应数据表列rawr_g_ratio位于yeast_raw_data_table_20180826.tsvcleanedlog2_cleaned_ratio位于yeast_data_table_20180826.tsvthresholdedlog2_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 给出了官方推荐的执行顺序运行Common functions中的代码单元加载全部辅助函数与常量定义这是后续所有阶段的前置条件可选运行Process the raw data从原始读数重建处理数据用于验证处理流程可被完整复现运行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作用在表达向量上返回每个时间过程 t0 点取值复制到整条序列的向量diff有限差分算子对时间过程内部使用对称差分相邻时间点ddt 0.5 / (shifted_times - times)对实验边界处置 0diff_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/ZEVCHIP_IDX在时间过程之上叠加time标识单张微阵列芯片TIMESERIES_IDX进一步叠加基因名标识每个基因 × 芯片数据点DESIGN_IDX在芯片基础上追加seq时间点在该时间过程内的序号与rseq倒序序号供差分/过滤逻辑使用。4.3 从原始信号到设计表Process the raw data阶段的数据整形流程为读取并验证读取 TSV验证r_g_ratio列可由双通道信号在 2.0 截断后相除重建跨斑点聚合用pivot_table以TIMESERIES_IDX为索引对r_g_ratio、rProcessedSignal、gProcessedSignal分别计算len、median、min、max、std微阵列通常每个基因有 2 个重复斑点取中位数聚合见 聚合单元异常值修剪trim_timeseries_outliers用两类规则修正离群点源码斑点不一致MAX_SPOT_RATIO 4.0同一基因在芯片上多个斑点amax/amin 4时用前后时间点的几何平均端点用邻居值替换中位数时间趋势尖峰JUMP_THRESHOLD 4.0非端点位置出现先上后下或先下后上且每步超过 4 倍的跳变时同样以相邻点几何平均替换构建设计表design_table以CHIP_IDX为行索引、GeneName为列抽取r_g_median、red_median、green_median三个块随后用timeseries_sequence计算seq/rseq列并按DESIGN_IDX排序最后dropna(axis1)剔除 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_dmnormalize_to_t0以基因表达中位数r_g_median为基准构造log_ratio log(r_g_median / (x0_op · r_g_median))即每个时间点的对数比相对于该时间过程 t0 的比值repair_crosstalk_outliers针对绿通道信号暴涨8 倍且未校正对数比同步上涨log 2的基因把log(green/green0)加回log_ratio从而把串扰造成的假信号纠正为以 t0 绿通道为基准的表达式修复时打印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_THRESHOLDING0.1)的置 0标记显著时间过程并剥离对剩余非显著数据逐芯片计算中央五分位误差t0 点芯片误差置 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软阈值thmedian(xnoise, x−noise, 0)即经典的软收缩一维阈值使用 t0 基因噪声广播threshold_error二维阈值2d使用完整的基因 × 芯片噪声矩阵threshold_error_2d序列过滤filt激进模式额外剔除只有一次检出且不在末点或检出存在断点的劣质时间过程apply_hard_and_soft_thresholding(..., aggressiveTrue)。阈值化过程中会打印检出统计量例如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 thresholdedu为目标基因列二次项uu u * u_m设计矩阵d_m hstack([u_m, uu])即一阶 二次双线性的调控模型稳态参考d_m - 1.0使截距表示t0 偏离稳态delta_modelTrue时假设稳态截距为 0对数模型log_weightsTrue时设计矩阵除以目标基因u因变量切换为log_cleaned_ratio对应拟合d/dt ln(gene)差分/积分模式inverse_designFalse时用diff_op对因变量求差分并过滤每条序列的末点因末点导数线性相关权重weights 1/(diff_op**2).sum(axis1)inverse_designTrue时则用inv_diff_op对设计矩阵积分、过滤 t0 行返回(设计矩阵, 因变量, 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, dtypebool) out_of_sample_f np.full(len(df_dm.index), False, dtypebool) 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_modelTrue, # 假设稳态截距为 0 thresholded_zth2dfilt, # 使用过滤后的二维阈值数据 model_quadratic_termsTrue, log_weightsTrue, # 拟合 d/dt ln(gene) inverse_designFalse, ...)拟合使用 sklearn 的 LARS最小角回归larsmodel sklearn.linear_model.Lars(fit_interceptFalse, normalizeFalse) 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 正则化路径选取第lambdacol40步的非零系数作为调控关系集合。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选定模型强制截距为 0delta_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(TFAFT1) hits np.sum(np.abs(np.log(aft1)) 0.2) 0 aft1[aft1.columns[hits]]即仅展示在 AFT1 时间过程中对数表达变化超过0.2约 22% 变化的调控因子列用于快速定位该时间过程中的主导调控因子。七、复现注意事项与适用前提结合 README 与 Notebook 源码复现或改造时需注意以下几点数据与模型文件必须齐备datadir下需同时具备 2 个 TSV 与 2 个 CSV若只跑Read processed data之后的阶段TSV 是必需的CSV 在Prediction model阶段必需两条入口不可混淆直接读官方处理表时不要运行Verify the processing procedure否则拿本地重建与官方表比对属于无效验证而重建路径下该验证是检验管线正确性的唯一手段差值min/max应为 0 量级GEV/ZEV 芯片基因集差异dropna(axis1)会把两套芯片共有的转录本之外的基因剔除因此最终公共基因集小于两套芯片的并集时间序列格式约束所有时间过程必须以 0 起始且严格递增construct_difference_operators会在违反时抛异常模型数值与论文的差异论文使用 Fortranglmnet封装Notebook 的 sklearnLars仅用于流程演示权重方案无法完全等价重现论文精确数值需使用论文原始工具链验证实验开关若目标是复现论文第一轮验证实验之前的结果需启用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创作声明：本文部分内容由AI辅助生成（AIGC），仅供参考 酵母转录调控网络微阵列数据分析全流程：从原始数据清洗到因果表达网络线性建模（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 NotebookGoogle 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 明确给出了论文使用的三个数据集与数据表中列的对应关系这是理解整条管线的关键索引数据集对应数据表列rawr_g_ratio位于yeast_raw_data_table_20180826.tsvcleanedlog2_cleaned_ratio位于yeast_data_table_20180826.tsvthresholdedlog2_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 给出了官方推荐的执行顺序运行Common functions中的代码单元加载全部辅助函数与常量定义这是后续所有阶段的前置条件可选运行Process the raw data从原始读数重建处理数据用于验证处理流程可被完整复现运行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作用在表达向量上返回每个时间过程 t0 点取值复制到整条序列的向量diff有限差分算子对时间过程内部使用对称差分相邻时间点ddt 0.5 / (shifted_times - times)对实验边界处置 0diff_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/ZEVCHIP_IDX在时间过程之上叠加time标识单张微阵列芯片TIMESERIES_IDX进一步叠加基因名标识每个基因 × 芯片数据点DESIGN_IDX在芯片基础上追加seq时间点在该时间过程内的序号与rseq倒序序号供差分/过滤逻辑使用。4.3 从原始信号到设计表Process the raw data阶段的数据整形流程为读取并验证读取 TSV验证r_g_ratio列可由双通道信号在 2.0 截断后相除重建跨斑点聚合用pivot_table以TIMESERIES_IDX为索引对r_g_ratio、rProcessedSignal、gProcessedSignal分别计算len、median、min、max、std微阵列通常每个基因有 2 个重复斑点取中位数聚合见 聚合单元异常值修剪trim_timeseries_outliers用两类规则修正离群点源码斑点不一致MAX_SPOT_RATIO 4.0同一基因在芯片上多个斑点amax/amin 4时用前后时间点的几何平均端点用邻居值替换中位数时间趋势尖峰JUMP_THRESHOLD 4.0非端点位置出现先上后下或先下后上且每步超过 4 倍的跳变时同样以相邻点几何平均替换构建设计表design_table以CHIP_IDX为行索引、GeneName为列抽取r_g_median、red_median、green_median三个块随后用timeseries_sequence计算seq/rseq列并按DESIGN_IDX排序最后dropna(axis1)剔除 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_dmnormalize_to_t0以基因表达中位数r_g_median为基准构造log_ratio log(r_g_median / (x0_op · r_g_median))即每个时间点的对数比相对于该时间过程 t0 的比值repair_crosstalk_outliers针对绿通道信号暴涨8 倍且未校正对数比同步上涨log 2的基因把log(green/green0)加回log_ratio从而把串扰造成的假信号纠正为以 t0 绿通道为基准的表达式修复时打印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_THRESHOLDING0.1)的置 0标记显著时间过程并剥离对剩余非显著数据逐芯片计算中央五分位误差t0 点芯片误差置 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软阈值thmedian(xnoise, x−noise, 0)即经典的软收缩一维阈值使用 t0 基因噪声广播threshold_error二维阈值2d使用完整的基因 × 芯片噪声矩阵threshold_error_2d序列过滤filt激进模式额外剔除只有一次检出且不在末点或检出存在断点的劣质时间过程apply_hard_and_soft_thresholding(..., aggressiveTrue)。阈值化过程中会打印检出统计量例如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 thresholdedu为目标基因列二次项uu u * u_m设计矩阵d_m hstack([u_m, uu])即一阶 二次双线性的调控模型稳态参考d_m - 1.0使截距表示t0 偏离稳态delta_modelTrue时假设稳态截距为 0对数模型log_weightsTrue时设计矩阵除以目标基因u因变量切换为log_cleaned_ratio对应拟合d/dt ln(gene)差分/积分模式inverse_designFalse时用diff_op对因变量求差分并过滤每条序列的末点因末点导数线性相关权重weights 1/(diff_op**2).sum(axis1)inverse_designTrue时则用inv_diff_op对设计矩阵积分、过滤 t0 行返回(设计矩阵, 因变量, 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, dtypebool) out_of_sample_f np.full(len(df_dm.index), False, dtypebool) 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_modelTrue, # 假设稳态截距为 0 thresholded_zth2dfilt, # 使用过滤后的二维阈值数据 model_quadratic_termsTrue, log_weightsTrue, # 拟合 d/dt ln(gene) inverse_designFalse, ...)拟合使用 sklearn 的 LARS最小角回归larsmodel sklearn.linear_model.Lars(fit_interceptFalse, normalizeFalse) 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 正则化路径选取第lambdacol40步的非零系数作为调控关系集合。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选定模型强制截距为 0delta_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(TFAFT1) hits np.sum(np.abs(np.log(aft1)) 0.2) 0 aft1[aft1.columns[hits]]即仅展示在 AFT1 时间过程中对数表达变化超过0.2约 22% 变化的调控因子列用于快速定位该时间过程中的主导调控因子。七、复现注意事项与适用前提结合 README 与 Notebook 源码复现或改造时需注意以下几点数据与模型文件必须齐备datadir下需同时具备 2 个 TSV 与 2 个 CSV若只跑Read processed data之后的阶段TSV 是必需的CSV 在Prediction model阶段必需两条入口不可混淆直接读官方处理表时不要运行Verify the processing procedure否则拿本地重建与官方表比对属于无效验证而重建路径下该验证是检验管线正确性的唯一手段差值min/max应为 0 量级GEV/ZEV 芯片基因集差异dropna(axis1)会把两套芯片共有的转录本之外的基因剔除因此最终公共基因集小于两套芯片的并集时间序列格式约束所有时间过程必须以 0 起始且严格递增construct_difference_operators会在违反时抛异常模型数值与论文的差异论文使用 Fortranglmnet封装Notebook 的 sklearnLars仅用于流程演示权重方案无法完全等价重现论文精确数值需使用论文原始工具链验证实验开关若目标是复现论文第一轮验证实验之前的结果需启用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创作声明：本文部分内容由AI辅助生成（AIGC），仅供参考