Epistasis Field GuideEFG拟南芥第一轮 · 试跑Arabidopsis round 1 · trial run
Step 2 · 互作检测Step 2 · Interaction search

互作检测:流程、候选和误报Interaction search: pipeline, candidates and false positives

真实数据Real data0 对越过门槛pairs over the bar

只指已运行的设置:A1、A2、B r² 0.2(固定方差比工作口径)、旧 C1 的 35 基因小例子。B r² 0.5 和新 C 未在真实数据上运行。Only the settings that have run: A1, A2, B r² 0.2 (fixed-variance-ratio working setting) and the old C1 35-gene example. B r² 0.5 and the new C have not run on real data.

模拟找回 · BSimulated recovery · B0/80 r² 0.2 · 25/80 r² 0.5

同一批 80 个有互作的模拟性状,总体找回率。Overall recovery on the same 80 simulated traits with an interaction.

S0 误报 · BS0 false positives · B8/200 两个面板相同on both panels

95% 区间 1.7%–7.7%。这批 S0 没有检出超出 5% 的证据,不能证明 ≤ 5%。95% interval 1.7%–7.7%. No evidence of exceeding 5% in this S0 batch; not proof of ≤ 5%.

待决定:Open: 选用 r² 0.2 还是 r² 0.5 面板(第 48 周,用户和导师);基因策略 C 的其余规格(与合作者商定)。whether to use the r² 0.2 or r² 0.5 panel (week 48, student and supervisor); the rest of gene strategy C (with a collaborator).

先修Before you start
Step 0、Step 1;Bonferroni 校正(总览第 02 节)Steps 0 and 1; the Bonferroni correction (overview, section 02)
学完能做After this page you can
沿着一条策略的流程说出每一关丢了什么,判断一个位点对为什么没过门槛,区分覆盖、条件检出和总体找回Follow one strategy's pipeline and say what each gate loses, judge why a pair did not cross the bar, and tell coverage, conditional detection and overall recovery apart
预计用时Time
约 30 分钟about 30 minutes
01 · 检验模型01 · The test model

同一面板上的 SNP 对用同一个模型检验SNP pairs on one panel are tested with one model

在同一个 SNP 面板上,A1、A2、B 送来的每一对都用下面的模型检验。C 的情况不同,见后面的表。On a given SNP panel, every pair from A1, A2 or B is tested with the model below. C is different; see the table further down.

现蕾时间 = 截距 + SNP₁ + SNP₂ + SNP₁ × SNP₂ + 亲缘 + 误差
days to buds = intercept + SNP₁ + SNP₂ + SNP₁ × SNP₂ + kinship + error
两个主效应Two main effects

模型先计入每个 SNP 各自的效应。乘积项检验的是两者各自效应之外的部分,所以两个单独就有效应的 SNP 不会被误判为互作。The model first accounts for each SNP's own effect. The product term then tests what is left beyond the two separate effects, so two SNPs with strong separate effects are not mistaken for an interaction.

乘积项Product term

两个位点都为 1 时取 1,否则为 0。它的系数等于总览页的交互对比。检验这个系数是否为 0,得到这一对的 p 值。It is 1 when both sites are 1 and 0 otherwise. Its coefficient equals the interaction contrast on the overview page. Testing whether the coefficient is zero gives the pair's p-value.

亲缘项(K)Kinship (K)

同一地区的品系基因型相近,现蕾时间也相近。不考虑这一点,地理来源会表现为假信号。模型因此加入由全基因组算出的。Lines from the same region have similar genotypes and similar times to buds. Without accounting for this, geographic origin shows up as a false signal, so the model includes a genome-wide .

模型和门槛的适用范围Where the model and bar apply

策略Strategy检验对象What is tested门槛的分母Bar denominator门槛Bar状态Status
A1、A2、B · r² 0.2SNP 对,乘积项,1 个自由度SNP pair, product term, 1 df面板内四格合格的对:1,742,876,507Eligible pairs on the panel: 1,742,876,5072.87 × 10⁻¹¹1a 已运行Run in 1a
B · r² 0.5同上Same16,449,009,2803.04 × 10⁻¹²只在模拟中运行Simulation only
C1 · 旧试用规格C1 · old trial spec基因对:基因内主成分乘积项的联合检验,自由度可大于 1Gene pair: joint test of products of within-gene components; df can exceed 1≥ 5 个 SNP 的基因两两组合:368,466,231Pairs of genes with ≥ 5 SNPs: 368,466,2311.36 × 10⁻¹⁰oracle(已知答案后直接检验真实的对)与 35 基因小例子Oracle (true pair tested directly once the answer is known) and 35-gene example
C · 新规格C · new spec294 个开花基因之间的基因对,基因 ± 500 bp;检验方法待商定Pairs among 294 flowering genes, gene ± 500 bp; test still to be agreed上界 43,071 对,正式分母待定At most 43,071 pairs; formal denominator open≈ 1.2 × 10⁻⁶(换算) (arithmetic)未运行Not run

只有同一面板上的 A1、A2、B 共用模型和门槛。换面板,分母随之改变;C 的检验对象和检验方法都不同。Only A1, A2 and B on the same panel share the model and bar. A different panel has a different denominator; C differs in both what is tested and how.

四格规则:每格至少 10 个品系Four-cell rule: at least 10 lines per cell

某一格只有两三个品系时,只要它们恰好都现蕾很晚,统计量就会很极端。所以四格都至少有 10 个品系的对才算合格,门槛的分母也只计合格的对(r² 0.2 面板约 17.4 亿对)。If one cell holds only two or three lines that happen to flower late, the statistic becomes extreme. A pair is therefore eligible only if every cell holds at least 10 lines, and the bar's denominator counts only eligible pairs (about 1.74 billion on the r² 0.2 panel).

02 · 策略 B 的流程02 · Strategy B step by step

17.8 亿组合初筛,1,000 对复检,工作口径下 0 对达标1.78 billion pairs screened, 1,000 re-tested, none over the bar in the working setting

① 全部组合① All combinations

r² 0.2 面板有 59,704 个 SNP,两两组合共 1,782,253,956 对。逐对用正式模型检验太慢,所以先做快速初筛。The r² 0.2 panel has 59,704 SNPs, giving 1,782,253,956 pairs. Testing each with the full model is too slow, so a fast screen comes first.

② 快速初筛(PLINK)② Fast screen (PLINK)

PLINK 的 --epistasis 每秒约检验 1,000 万对,但不能加入亲缘项。所以先用零模型从现蕾时间中减去亲缘部分,在上初筛。p < 10⁻⁴ 的 54,158 对被写出。PLINK's --epistasis tests about 10 million pairs a second but cannot include kinship. A null model first removes the kinship part from time to buds, and the screen runs on the . The 54,158 pairs with p < 10⁻⁴ are written out.

③ 四格规则和秩检查③ Four-cell rule and rank check

写出的对中有 26% 不满足每格 ≥ 10 个品系。在残差上,品系很少的格子容易产生极端值。剩下的 40,288 对都通过了,即乘积项确实带来了新的信息。26% of the written pairs fail the ≥ 10 lines per cell rule; on residuals, cells with very few lines readily give extreme values. All 40,288 remaining pairs pass the , meaning the product term adds new information.

④ 预算:前 1,000 对④ Budget: the top 1,000

按初筛 p 值排序,只用正式模型复检前 1,000 对。预算决定计算量和漏检之间的平衡。在一个 500 个 SNP 的测试面板上,初筛与正式模型的排序相关系数为 0.985,正式模型的前 50 名都在初筛的前 1,000 名之内。Pairs are sorted by screening p-value and only the top 1,000 are re-tested with the full model. The budget sets the balance between computing and misses. On a 500-SNP test panel, screen and full-model rankings correlate at 0.985, and the full model's top 50 all fall within the screen's top 1,000.

⑤ 门槛:工作口径下 0 对⑤ The bar: 0 pairs in the working setting

门槛为 2.9 × 10⁻¹¹。最接近的一对在工作口径(固定方差比)下 p = 5.2 × 10⁻¹¹,没有越过。逐对重估方差比后 p = 9.0 × 10⁻¹⁴,越过了门槛。但它的信号来自一个只有 13 个品系的格子,其中 6 个来自意大利;后 p = 1.3 × 10⁻⁸,品系后最差 p = 3.8 × 10⁻⁹,都明显变弱。它属于在敏感性分析下过阈值、稳健性存疑的探索性候选。是否把稳健性作为发现的附加条件,待导师讨论。下一节可以查看这一对。The bar is 2.9 × 10⁻¹¹. Under the working setting (fixed variance ratio) the closest pair has p = 5.2 × 10⁻¹¹ and does not clear it. With a per-pair variance ratio, p = 9.0 × 10⁻¹⁴, which does. But the signal comes from a cell of only 13 lines, 6 of them from Italy, and it weakens clearly after (p = 1.3 × 10⁻⁸) and (worst p = 3.8 × 10⁻⁹). It is an exploratory candidate that clears the bar in a sensitivity analysis but has doubtful robustness. Whether robustness should be a condition for a finding is open for discussion with the supervisor. The next section shows this pair.

03 · 位点对浏览器03 · Pair browser

查看每一对在 936 个品系中的四格分布See each pair's four cells across the 936 lines

每条弧线是真实分析中检验过的一对 SNP。点击一条弧线,右侧用 936 个品系的真实现蕾时间画出它的,并列出 p 值、稳健性检查和附近的基因。下方四个按钮是有代表性的例子。Each arc is a SNP pair tested in the real analysis. Click an arc and the panel on the right plots its from the real times to buds of 936 lines, with p-values, robustness checks and nearby genes. The four buttons below are representative examples.

A1、A2 按单个 SNP 的信号挑选,弧线集中在 5 号染色体已知的主效区域(FLC、DOG1 附近)。B 不参考单 SNP 信号,弧线分布在全基因组。这里列出的是真实数据上检验过的对:A1 为严格版,A2 为筛选版。Step 5 比较的是 A1、A2 筛选版在模拟上的结果。A1 and A2 select by single-SNP signal, so their arcs concentrate in the known main-effect regions on chromosome 5 (near FLC and DOG1). B does not use single-SNP signals and its arcs spread across the genome. The pairs shown are those tested on real data: the strict version of A1 and the screened version of A2. Step 5 compares the screened versions of A1 and A2 on simulations.

04 · 面板与模拟04 · Panels and simulation

面板越小,计算越快,覆盖到的因果位点对也越少A smaller panel is faster to scan and covers fewer causal pairs

相邻的 SNP 倾向于一起遗传,称为(LD),用 r² 衡量,取值 0 到 1。r² 高的两个 SNP 几乎总是同时出现,保留其中一个就能代替另一个,这个 SNP 叫。面板是从 181 万个常见 SNP 中精简出来的:彼此 r² 超过阈值的相邻 SNP 只保留一个。Neighbouring SNPs tend to be inherited together. This is (LD), measured by r² from 0 to 1. Two SNPs with high r² almost always occur together, so one can stand in for the other; it is then called a . A panel is thinned from the 1.81 million common SNPs: of neighbours whose r² exceeds the cut-off, only one is kept.

阈值越低,面板越小,计算越快。但因果 SNP 可能在面板里没有标签,B 就无法检验到它。The lower the cut-off, the smaller and faster the panel. But a causal SNP may then have no tag in the panel, and B cannot test it.

第 48 周由用户和导师决定Decided in week 48 by student and supervisor

覆盖率用 2,000 个随机抽取的因果位点对估计;初筛时间为每个性状在 32 线程下的实测墙钟时间。r² 0.8 只测了覆盖率和时间,没有运行完整流程。Coverage was estimated from 2,000 randomly drawn causal pairs; screening time is the measured wall time per trait on 32 threads. For r² 0.8 only coverage and time were measured, not the full pipeline.

同一批 80 个模拟性状在两个面板上的结果The same 80 simulated traits on both panels

图中每一行是一关,数字是埋入的互作对通过这一关的个数:覆盖、初筛写出、资格过滤与预算、越过门槛。Each row is one gate, and the number is how many planted pairs pass it: covered, written by the screen, eligibility and budget, over the bar.

第三关同时包含四格规则、秩检查和只复检前 1,000 对,现有记录无法把三者拆开。第四关的"旧门槛"数字是固定候选的诊断对照;按本轮预设门槛,强互作检出 20/40,弱互作检出 5/40。r² 0.5 的初筛时间约为 r² 0.2 的 7 倍,门槛约严 10 倍。The third gate combines the four-cell rule, the rank check and the top-1,000 re-test budget, and the records cannot separate the three. The "old bar" figures at the fourth gate are a fixed-candidate diagnostic; with this round's pre-set bar, 20/40 strong and 5/40 weak interactions are detected. Screening on r² 0.5 takes about 7 times as long as on r² 0.2, and the bar is about 10 times stricter.

逐个查看 80 个性状The 80 traits one by one

每一列是一个,左边 40 个强互作,右边 40 个弱互作;每一行是一条策略。颜色表示埋入的那一对走到了哪一关。点击一列或用下拉框选择,可以看到这个性状的答案和各策略的结果。Each column is a , with 40 strong interactions on the left and 40 weak ones on the right; each row is a strategy. Colour shows how far the planted pair got. Click a column or use the drop-down to see that trait's answer and each strategy's outcome.

05 · 误报05 · False positives

200 个无互作性状上的假发现False findings on 200 traits without interaction

每个点代表一个不含互作的模拟性状(S0)。有颜色的点表示这个性状在整套搜索中至少有一对越过门槛,算一次假发现。判定规则在运行前写定,看到结果后没有修改。点击一个点或用下拉框,可以看到它的最小 p 值。Each dot is a simulated trait without interaction (S0). A coloured dot means at least one pair in the whole search cleared the bar on that trait, which counts as one false finding. The rules were fixed before the run and not changed afterwards. Click a dot or use the drop-down to see its smallest p-value.

r² 0.2 面板的假发现False finding, r² 0.2 panel r² 0.5 面板的假发现False finding, r² 0.5 panel
8 / 200r² 0.2:4.0%(95% 区间 1.7%–7.7%)r² 0.2: 4.0% (95% interval 1.7%–7.7%)
8 / 200r² 0.5:4.0%(1.7%–7.7%)r² 0.5: 4.0% (1.7%–7.7%)
0两个面板上都有假发现的性状Traits with a false finding on both panels

目标是不超过 5%。95% 包含 5%,说明这批 S0 没有检出超出 5% 的证据;区间上界为 7.7%,所以不能证明一定 ≤ 5%。两个面板配对来看,都有假发现的 0 个,只有 r² 0.2 有的 8 个,只有 r² 0.5 有的 8 个,都没有的 184 个;这张表只作描述,不能说明假发现是怎样产生的。S0 性状仍含加性背景效应,只是没有互作。A1、A2 在另一批 S0(0.7 批次 10 个、周末批次 40 个)上是 0/50(区间 0–7.1%),它们检验的对少得多,不能据此认为误报更低。The target is a of at most 5%. The 95% contains 5%, so this S0 batch shows no evidence of exceeding 5%; its upper bound is 7.7%, so it does not prove the rate is at most 5%. Paired across panels: 0 traits have a false finding on both, 8 only on r² 0.2, 8 only on r² 0.5 and 184 on neither. This table is descriptive and does not show how false findings arise. S0 traits still carry additive background effects; they lack only the interaction. A1 and A2 scored 0/50 on a different S0 batch (10 from the 0.7 batch and 40 from the weekend batch; interval 0–7.1%); they test far fewer pairs, so this does not show a lower error rate.

06 · 复核中发现的问题06 · Problems found in review

四个在复核中发现并已更正的问题Four problems found in review and corrected

分析由一个 AI(Claude)执行,另一个 AI(Codex)独立复核。以下问题都已在运行记录中更正。One AI (Claude) ran the analyses and another (Codex) reviewed them independently. All the problems below have been corrected in the run records.

1

按规则算作假发现的 13 对,一端几乎就是因果 SNPThirteen pairs counted as false by the rule had one end almost at the causal SNP

一个模拟性状上有 13 对越过门槛。按事先的规则(标签须在因果 SNP ± 20 kb 内),它们全部算作假发现。逐个检查发现,每一对的一端几乎就是一个因果 SNP,另一端离另一个因果 SNP 78–920 kb,但仍与它连锁;最强的一对相距 634 kb,r² 0.84。这与 LD 代理的解释相容,即远处的 SNP 通过连锁反映真实信号,但机制没有验证。On one simulated trait, 13 pairs cleared the bar. Under the pre-set rule (a tag must lie within ± 20 kb of the causal SNP) all counted as false. Checked one by one, each pair had one end almost at a causal SNP and the other 78–920 kb from the second causal SNP yet still linked to it; the strongest was 634 kb away with r² 0.84. This is compatible with an LD-proxy explanation, in which a distant SNP reflects the true signal through linkage, but the mechanism was not verified.

处理:主规则不变,看到结果后不改规则;另在事先写定一条信号归并的敏感性规则,两种口径都报告。Correction: the main rule stays, since rules are not changed after seeing results; a signal-merging sensitivity rule was written in advance, and both versions are reported.

2

最接近门槛的候选依赖一个 13 个品系的格子The candidate closest to the bar depends on a 13-line cell

真实数据上最接近门槛的一对,信号来自一个只有 13 个品系的格子,其中 6 个来自意大利,而意大利品系只占全部的 7.8%。信号可能来自少数特殊品系,而不是两个位点之间普遍存在的互作。The pair closest to the bar on real data draws its signal from a cell of only 13 lines, 6 of them from Italy, although Italian lines are only 7.8% of the total. The signal may come from a few unusual lines rather than a general interaction between the two sites.

处理:每个候选都报告稳健性检查(秩变换、逐个剔除品系)。在浏览器中打开这一对Correction: every candidate report now includes robustness checks (rank transformation, leave-one-out). Open this pair in the browser

3

弱互作的漏检原因最初归错了The losses of weak interactions were first misattributed

最初的记录写 r² 0.5 的漏检主要来自预算和更严的门槛。复核指出,"写出到进入复检"这一关同时包含四格过滤和预算。固定候选、只换门槛后发现,弱互作的 12 次漏检中有 11 次在旧门槛下也通不过,门槛变严只多漏了 1 次。The first record said the losses on r² 0.5 came mainly from the budget and the stricter bar. Review pointed out that the step from written to re-tested combines the four-cell filter with the budget. Holding the candidates fixed and changing only the bar showed that 11 of the 12 weak-interaction losses fail even at the old bar; the stricter bar adds just 1.

做法:逐关归因,并用每次只改一个因素的诊断来验证。Practice: attribute losses gate by gate, and check with a diagnostic that changes one factor at a time.

4

验收脚本会漏掉放错目录的文件The acceptance script missed files in the wrong folder

200 个结果文件分布在三台服务器的目录里。验收脚本先把所有文件名合并成一个名单,所以放错目录的文件即使含有坏数据,也会被跳过,脚本仍报告通过。复核者用一份小数据复现了这个问题。The 200 result files were spread over three servers' folders. The acceptance script first merged all file names into one list, so a file in the wrong folder was skipped even if it held bad data, and the script still reported a pass. The reviewer reproduced this with a small data set.

处理:改为逐个目录核对,缺少文件即判失败;用反例测试后重新验收,结果仍为 8/200。检查脚本本身也要用故意构造的坏数据测试。Correction: the script now checks folder by folder and fails on any missing file; after testing on the counterexamples it was re-run, still giving 8/200. A checking script should itself be tested on deliberately broken data.

07 · 学习清单07 · Study list

理解这一步需要的七个概念Seven concepts behind this step

勾选只保存在你自己的浏览器里。Ticks are saved only in your own browser.

练习 1:模拟一对互作,再检验它Exercise 1: simulate an interaction, then test it

在 Terminal 里运行,需要 numpy 和 scipy(Anaconda 自带)。代码模拟 936 个品系和 20 天的互作效应,打印四格的品系数和平均值,再拟合带乘积项的模型。这个简化版本没有亲缘项。Run in Terminal; it needs numpy and scipy, both included with Anaconda. The code simulates 936 lines with a 20-day interaction, prints the line count and mean of each cell, and fits the model with the product term. This simplified version has no kinship term.

python3 - <<'EOF' import numpy as np from scipy import stats rng = np.random.default_rng(1) n = 936 x1 = rng.binomial(1, 0.3, n) x2 = rng.binomial(1, 0.4, n) y = 60 + 10*x1 + 15*x2 + 20*x1*x2 + rng.normal(0, 25, n) for a in (0, 1): for b in (0, 1): m = (x1 == a) & (x2 == b) print(f"cell {a}{b}: {m.sum():4d} lines, mean {y[m].mean():5.1f} days") X = np.column_stack([np.ones(n), x1, x2, x1*x2]) beta, rss, *_ = np.linalg.lstsq(X, y, rcond=None) se = np.sqrt(rss[0] / (n - 4) * np.linalg.inv(X.T @ X)[3, 3]) t = beta[3] / se print("interaction:", round(beta[3], 1), "days, p =", f"{2 * stats.t.sf(abs(t), n - 4):.2g}") EOF

然后把 20*x1*x2 改成 5*x1*x2,再把 0.3 改成 0.05,使其中一格变得很小。观察 p 值的变化,以及它离 10⁻¹¹ 的门槛有多远。Then change 20*x1*x2 to 5*x1*x2, and 0.3 to 0.05 so that one cell becomes very small. Watch how the p-value changes and how far it is from a 10⁻¹¹ bar.

练习 2:门槛和 8/200 的置信区间Exercise 2: the bar and the 8/200 interval

python3 -c "from scipy.stats import binomtest; print(0.05/1742876507); print(binomtest(8, 200).proportion_ci(method='exact'))"

结果应约为 2.87e-11 和 (0.0174, 0.0773)。再运行 binomtest(0, 50),想一想为什么 0 次假发现的区间上界仍有 7.1%。You should get about 2.87e-11 and (0.0174, 0.0773). Then run binomtest(0, 50) and consider why zero false findings still leaves an upper bound of 7.1%.

08 · 待决定的事项08 · Open questions
第 48 周 · 用户和导师Week 48 · student and supervisor

选用 r² 0.2 还是 r² 0.5 面板The r² 0.2 or the r² 0.5 panel

现有证据:当前模拟中 r² 0.5 检出更多;这批 S0 没有检出超出 5% 的证据,但不能证明 ≤ 5%。代价是计算量更大、门槛更严。r² 0.5 还没有在真实数据和预测中运行。Current evidence: r² 0.5 detects more in the current simulations; this S0 batch shows no evidence of exceeding 5%, but does not prove ≤ 5%. The cost is more computing and a stricter bar. r² 0.5 has not yet run on real data or in prediction.

与合作者商定With a collaborator

基因策略 C 的其余规格The rest of gene strategy C

暂定:只用开花基因,基因 ± 500 bp。待定:开花基因名单、基因内频率下限、每个基因的主成分数、检验方法、相邻基因的处理。Provisional: flowering genes only, gene ± 500 bp. Open: the gene list, the within-gene frequency cut-off, components per gene, the test, and how to treat neighbouring genes.

第 48 周 · 导师Week 48 · supervisor

是否增加因果位点对取自开花基因的模拟Whether to add simulations with causal pairs in flowering genes

新 C 只检验 294 个开花基因,而现有模拟随机抽取因果 SNP,80 对中没有一对的两端都在开花基因里。这是使用先验的代价,照常报告。如果新增这类模拟,所有策略都要用同一批性状。The new C tests only 294 flowering genes, while the current simulations draw causal SNPs at random; none of the 80 pairs has both ends in a flowering gene. This is the cost of using a prior and is reported as such. Any new scenario of this kind must be shared by all strategies.

原因未明Unexplained

两个原因未明的现象Two unexplained observations

真实数据上随机对的检验统计量偏低(λ ≈ 0.89)。在周末批次 r² 0.2 的有信号模拟中,11/80 = 13.75% 的性状出现了按事后规则归类为弱 LD、远距离的达阈值对。是否在冻结前排查,待定。On real data, test statistics for random pairs are slightly deflated (λ ≈ 0.89). In the weekend-batch r² 0.2 signal simulations, 11/80 = 13.75% of traits had over-the-bar pairs classed after the fact as weak-LD and distant. Whether to investigate before the freeze is open.

本页产物Your output for this page

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

在位点对浏览器里选一对,写出:它来自哪个策略,p 值离门槛多远,最小的格子有多少品系,稳健性检查的结果如何,以及要把它当作发现还缺哪些证据。Pick a pair in the pair browser and write down: which strategy it comes from, how far its p-value is from the bar, how many lines the smallest cell has, what the robustness checks show, and what evidence is still missing before it could count as a finding.

状态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