Epistasis Field GuideEFG拟南芥第一轮 · 试跑Arabidopsis round 1 · trial run
Step 3 · 互作的贡献Step 3 · What interactions add

互作能解释多少变异,能不能改善预测?How much variation do interactions explain, and do they improve prediction?

3a · 双核模型3a · Two-kernel model0 DTF1 上 K∘K 的方差占比K∘K share on DTF1

估计落在 0 的边界上,与只含 K 的模型相同。只有加性背景的 S0 模拟里,K∘K 也会分到中位 5% 的方差。The estimate sits on the boundary at 0, the same as the K-only model. In S0 simulations with only an additive background, K∘K still takes a median 5% of the variance.

3b · 基准预测 M03b · Baseline M0r 0.822 · RMSE 14.9 天days

GBLUP,5 折平均。之后每个策略的预测都和它比较。GBLUP, mean over 5 folds. Every strategy's prediction is compared with it.

3b · 互作项的增益3b · Gain from interaction termsΔr ≤ 0 五折平均5-fold mean

Δr = M2 − M1 的相关系数之差。A1 −0.014、A2 −0.002、B −0.023、C1(35 基因小例子)−0.003;A2 和 C1 各有一折为正。限于当前候选、固定惩罚和单次 5 折。Δr is the difference in correlation, M2 − M1. A1 −0.014, A2 −0.002, B −0.023, C1 (35-gene example) −0.003; A2 and C1 each have one positive fold. Limited to the current candidates, a fixed penalty and a single 5-fold split.

待决定:Open: 随机对照是否同时匹配端点和基因的复用结构(1b,用户决定);r² 0.5 面板和新规格的 C 还没有在 Step 3 运行。whether random controls should also match how end points and genes are reused (1b, the student decides); the r² 0.5 panel and the new C have not run in Step 3.

先修Before you start
Step 2;知道交叉验证的基本做法Step 2; know the basic idea of cross-validation
学完能做After this page you can
区分全局方差分量(3a)和各策略的预测增益(3b),读出逐折结果和随机对照Tell the global variance component (3a) apart from each strategy's prediction gain (3b), and read the per-fold results against the random controls
预计用时Time
约 20 分钟about 20 minutes
01 · 两个问题01 · Two questions

3a 不分策略,3b 逐个策略比较3a ignores strategies; 3b compares them one by one

3a

全局方差分量Global variance components

在亲缘矩阵 K 之外再加一个核 K∘K(K 与自身逐元素相乘),估计 K∘K 能分走多少方差。这是对全部数据的一次全局估计,不依赖任何策略选出的候选对。Add a second kernel K∘K (K multiplied element by element with itself) next to the kinship matrix K, and estimate how much variance K∘K takes. This is one global estimate on all the data, independent of the pairs any strategy selects.

3b

折外预测增益Out-of-fold prediction gain

在训练折里让每个策略各挑出前 20 对,比较三个模型对测试折的预测:M0 只有 GBLUP;M1 加上这些对的端点主效应;M2 再加上乘积项。M2 − M1 衡量互作项本身带来的增益。In each training fold every strategy picks its top 20 pairs, and three models predict the test fold: M0 is GBLUP alone; M1 adds the pairs' end-point main effects; M2 also adds the product terms. M2 − M1 measures the gain from the interaction terms themselves.

02 · 3a 方差分量02 · 3a variance components

DTF1 上 K∘K 的分量为 0,而只有加性背景的模拟也会分到方差K∘K takes nothing on DTF1, yet takes variance in simulations with only an additive background

y = μ + gA + gAA + e,  gA ~ N(0, σ²A KA),  gAA ~ N(0, σ²AA KAA)

KA 是缩放后的 K,KAA 是缩放后的 K∘K。两个核非对角元的相关是 0.739,所以两个分量很难分开估计。拟合用自己写的带约束 AI-REML,单核特例与 GEMMA 一致。KA is the scaled K and KAA the scaled K∘K. Their off-diagonal elements correlate at 0.739, so the two components are hard to estimate separately. The fit uses a self-written constrained AI-REML, whose single-kernel special case agrees with GEMMA.

每个点是一个性状的 K∘K 方差占比 σ²AA ÷ (σ²A + σ²AA + σ²e),竖线是中位数。S0、S2、S3 由 Step 0 的生成方式得到:加性背景来自随机抽取的 1,000 个 SNP,不是从核模型生成的;S2、S3 另含一对互作 SNP(解释 10% 和 5% 的方差)。LRT 的 p 值按边界混合分布(χ²₀ 与 χ²₁ 各半)近似。Each dot is one trait's K∘K variance share, σ²AA ÷ (σ²A + σ²AA + σ²e); the vertical tick is the median. S0, S2 and S3 come from Step 0's generator: the additive background uses 1,000 randomly drawn SNPs and is not generated from the kernel model; S2 and S3 also contain one interacting SNP pair (10% and 5% of the variance). LRT p-values use the boundary mixture approximation (an equal mix of χ²₀ and χ²₁).

  • S0 的 K∘K 占比中位 5%,97.5% 分位约 31%,最大约 55%;40 次中 3 次 LRT p < 0.05。所以一组数据上 K∘K 不为 0,本身不能说明存在上位性。In S0 the K∘K share has a median of 5%, a 97.5% quantile of about 31% and a maximum of about 55%; 3 of 40 have LRT p < 0.05. So a non-zero K∘K on one data set does not by itself show epistasis.
  • S2、S3 的分布和 S0 差别不大,各 5/40 次 LRT p < 0.05。在当前的核、生成方式和检验设置下,双核模型对单个互作对的功效很低;这不能推广为"936 个品系的互作检验都没有功效",也不能证明不存在互作。S2 and S3 look much like S0, each with 5/40 LRT p < 0.05. Under the current kernels, generating model and test, the two-kernel model has low power for a single interacting pair; this does not mean every interaction test on 936 lines lacks power, nor does it show that interactions are absent.
  • 另一组检查直接从双核模型本身生成数据(与上图的 S0 不同):真值 σ²AA = 0 时,LRT 拒绝 2/50(4%,95% 区间 0.5%–13.7%);真值 0.2 时平均估计 0.175,偏差 −0.025,标准差 0.116。两个分量估计值在重复之间的经验相关是 −0.60 和 −0.78。A separate check generates data directly from the two-kernel model (unlike the S0 above): with true σ²AA = 0 the LRT rejects 2/50 (4%, 95% interval 0.5%–13.7%); with a true value of 0.2 the mean estimate is 0.175, a bias of −0.025 with SD 0.116. Across replicates the two components' estimates correlate at −0.60 and −0.78.
03 · 3b 折外预测03 · 3b out-of-fold prediction

本轮单次五折中,四组设置的平均 Δr 都没有增加In this round's single 5-fold split, the mean Δr did not increase for any of the four settings

5 折划分(种子 20261102,每折 187–188 个品系)对所有模型共用。表型相关的步骤都在训练折内完成:边际筛选、候选提名、四格和秩检查。每折取前 20 对,岭惩罚固定为该块的列数。对照是每折 10 组匹配的随机对,匹配 MAF、端点显著程度和最小格。One 5-fold split (seed 20261102, 187–188 lines per fold) is shared by all models. Every phenotype-dependent step happens inside the training fold: marginal screening, nominating candidates, and the four-cell and rank checks. Each fold takes the top 20 pairs, with the ridge penalty fixed at the block's number of columns. The control is 10 matched random groups per fold, matched on MAF, end-point significance and smallest cell.

真实候选(每折一行)Real candidates (one row per fold) 同折的 10 个随机组The same fold's 10 random groups 零对随机组(没有可匹配的对,M2 = M1 = M0)Empty random group (no matchable pair; M2 = M1 = M0)

横轴是 M2 − M1 的相关系数差,竖线是 0。每一行是一折,真实结果只和同一行的 10 个随机组比较。右侧是真实组高于随机组的次数(共 50 次比较,描述性计数,不是 50 次独立的交叉验证)。A2 的 50 个随机组里有 10 个是零对(第 4 折全部为空),它们不是有效的对照。C1 是 35 个基因的小例子,真实组和随机组的主效应列数相差约一倍;它的随机池是 3.6 版,PCA 回退修正后随机池可能改变,18/50 尚未在新池上重验。The x-axis is the difference in correlation, M2 − M1; the vertical line is 0. Each row is one fold, and the real result is compared only with the 10 random groups in the same row. The number on the right counts how often the real group beats a random group (50 comparisons; a descriptive count, not 50 independent cross-validations). 10 of A2's 50 random groups are empty (all of fold 4), so they are not valid controls. C1 is the 35-gene example, and its real and random groups differ about twofold in main-effect columns; its random pool is the version from run 3.6, the pool may change after the PCA fallback fix, and 18/50 has not been re-checked on the new pool.

逐折读数Fold-by-fold readings

折Fold真实组 ΔrReal Δr10 个随机组的 Δr:最小 · 中位 · 最大Δr of the 10 random groups: min · median · max真实组高于几个随机组Random groups below the real one零对随机组Empty random groupsΔRMSE

与上图同一组数据,按折列出精确值。用上面的按钮选策略,按钮可以用键盘的 Tab 和回车操作。Δr 显示到 4 位小数;“高于几个随机组”按完整精度计算,所以个别看起来相等的值也可能分出高低。零对随机组的 M2 与 M1 相同,Δr 为 0,计数时照常参与比较,所以 A2 第 4 折的计数要谨慎解读。各折计数加起来等于上图右侧的总数。The same data as the figure above, with exact values per fold. Choose a strategy with the buttons; they work with Tab and Enter. Δr is shown to 4 decimals; the count of random groups below the real one uses full precision, so values that look equal can still be ordered. An empty random group has M2 equal to M1 and Δr of 0; it still enters the count, so read A2's fold-4 count with care. The per-fold counts add up to the total on the right of the figure.

策略Strategy端点列 / 乘积列End-point / product columns平均 Δr(M2 − M1)Mean Δr (M2 − M1)Δr 为正的折Folds with Δr > 0ΔRMSE(天)ΔRMSE (days)高于随机组Above random

ΔRMSE = M2 − M1 的均方根误差差,正值表示 M2 更差。A2 的随机组有 280/950 对匹配不上、71 对放宽了最小格分箱,部分随机组为空,所以 A2 的 26/50 不是在严格匹配下得到的。ΔRMSE is the RMSE of M2 minus that of M1; a positive value means M2 is worse. For A2, 280 of 950 random pairs could not be matched and 71 had the smallest-cell bins relaxed, leaving some random groups empty, so A2's 26/50 was not obtained under strict matching.

04 · 结论的范围04 · Scope of the results

"没有增益"只适用于这一组设置"No gain" holds only for this set of settings

  • 候选空间:r² 0.2 面板上的 A1、A2、B,以及旧 C1 的 35 个基因。r² 0.5 面板和新规格的 C 没有运行。Candidate space: A1, A2 and B on the r² 0.2 panel, and the old C1's 35 genes. The r² 0.5 panel and the new C have not run.
  • 惩罚固定为列数,没有调过;要调参,必须在训练集内部再划验证集。The penalty is fixed at the number of columns and was not tuned; tuning would need an inner validation split within the training set.
  • 只有一次 5 折划分,各折共享样本,没有做显著性检验。There is a single 5-fold split; folds share samples, and no significance test was run.
  • 随机对照没有匹配端点和基因的复用结构。真实的前 20 对会反复用到同一些 SNP 或基因,随机组几乎不会,这对 C1 的 M1 比较影响最大。The random controls do not match how end points and genes are reused. The real top 20 pairs reuse the same SNPs or genes; random groups hardly do, which affects the C1 M1 comparison most.
  • "没有增益"不能推广为"互作对预测没有价值"。"No gain" does not generalise to "interactions have no value for prediction".
05 · 常见错误05 · Common mistakes

复核中发现、已更正的三个问题Three problems found in review and corrected

1

双核模型的对数似然比单核还低The two-kernel model had a lower log-likelihood than the single-kernel one

第一版在分量变负时把它截到很小的值,其余分量却仍按完整的牛顿步更新,结果收敛到错误的点。DTF1 上双核的 logL 比嵌套的单核低 30.3,这在最优点上不可能出现。The first version clipped a component that turned negative to a tiny value but still updated the others with the full Newton step, so it converged to the wrong point. On DTF1 the two-kernel logL was 30.3 below the nested single-kernel model, which is impossible at an optimum.

处理:改为带约束的 AI-REML,每次拟合都检查 logL双核 ≥ logL单核。Correction: a constrained AI-REML, with every fit checked for logLtwo ≥ logLone.

2

失败的拟合碰巧给出了正确的答案A failed fit happened to give the right answer

加入明确的失败处理后,DTF1 第一次被判为失败。旧版在线搜索耗尽时照样更新参数,碰巧回到了最优点,失败被掩盖了。Once explicit failure handling was added, DTF1 was flagged as failed at first. The old version updated the parameters even when the line search ran out, happened to land on the optimum, and hid the failure.

处理:步长只在评估后才接受,改为投影线搜索;失败时不报告 LRT。修正后 221 个拟合全部成功,数值与之前相同。Correction: a step is accepted only after evaluation, with a projected line search; failed fits report no LRT. After the fix all 221 fits succeed with unchanged values.

3

把 C1 的 M1 增益归因于基因清单C1's M1 gain was attributed to the gene list

最初写"C1 的 M1 增益主要反映基因清单本身"。核对列数后发现,真实组每折 87–92 个主效应列,随机组 170–180 个,差别可能来自列数而不是基因。It was first written that "C1's M1 gain mainly reflects the gene list". Checking the columns showed 87–92 main-effect columns per fold for the real group and 170–180 for the random groups, so the difference may come from the column count rather than the genes.

处理:这句收回;1b 建议让随机对照保持真实对的复用结构。Correction: the sentence was withdrawn; for 1b, random controls that keep the real pairs' reuse structure are proposed.

06 · 练习06 · Exercises

先算两核的相关First, compute the correlation between the two kernels

下面的代码用随机基因型算一个亲缘矩阵 K,再算 K∘K,比较两者非对角元的相关。The code builds a kinship matrix K from random genotypes, computes K∘K, and correlates their off-diagonal elements.

python3 -c " import numpy as np rng = np.random.default_rng(1) n, m = 300, 2000 Z = rng.binomial(1, rng.uniform(0.1, 0.5, m), (n, m)).astype(float) Z -= Z.mean(0) K = Z @ Z.T / m KK = K * K iu = np.triu_indices(n, 1) print('corr of off-diagonals, K vs K*K:', round(np.corrcoef(K[iu], KK[iu])[0, 1], 3)) "

随机基因型没有群体结构,相关很弱(这段代码约为 −0.20);真实数据有群体结构,936 个品系上是 0.739。再结合 02 节的检查:两核相关为 0.739 时,同一组数据的两个分量估计在重复之间的相关是 −0.60 到 −0.78,其中一个偏高时另一个往往偏低。这就是两者难以分开的原因。Random genotypes have no population structure and give a weak correlation (about −0.20 with this code); the real data have structure, and on the 936 lines it is 0.739. Combine this with the checks in section 02: with the kernels correlated at 0.739, the two component estimates correlate at −0.60 to −0.78 across replicates, so when one comes out high the other tends to come out low. That is why the two are hard to separate.

DTF1 上 K∘K 的方差估计为 0。能不能说 DTF1 没有上位性?The K∘K variance on DTF1 is estimated at 0. Does that show DTF1 has no epistasis?

本页产物Your output for this page

写下你自己的解释Write your own explanation

用两三句话分别说明 3a 的 K∘K 方差占比和 3b 的 M2 − M1 回答什么问题,各引用本页的一个数字,并说明为什么 DTF1 上 K∘K 为 0 不能说明没有上位性。In two or three sentences each, say which question the 3a K∘K variance share and the 3b M2 − M1 answer, quote one number from this page for each, and explain why a K∘K share of 0 on DTF1 does not show that there is no epistasis.

状态Status

只保存在这个浏览器里,换设备或清除浏览器数据后会丢失。所有页面的记录可以在总览页的“我的学习记录”里一起导出。Saved in this browser only; it is lost on another device or after clearing browser data. The records of all pages can be exported together under “My learning record” on the overview page. 去导出Go to export