Epistasis Field GuideEFG拟南芥第一轮 · 试跑Arabidopsis round 1 · trial run
Step 1 · 单 SNP 扫描Step 1 · Single-SNP scan

哪些 SNP 单独就和现蕾时间相关?Which SNPs are associated with time to buds on their own?

混合模型 MLMMixed model (MLM)0 个 SNP 过门槛SNPs over the bar

λ = 0.995。两个探索性峰:DOG1 内,以及 FLC 下游约 12 kb,都没有过门槛。λ = 0.995. Two exploratory peaks, inside DOG1 and about 12 kb downstream of FLC; neither clears the bar.

GLM · 10 个主成分GLM · 10 PCs3,642 个 SNP 过门槛SNPs over the bar

λ = 2.54,QQ 图明显偏离。存在结构或亲缘校正不足的疑虑,膨胀的来源目前无法分解。λ = 2.54 with a clearly deviating QQ plot. Structure or relatedness may be under-corrected; the source of the inflation cannot yet be separated.

SNP 遗传力SNP heritability0.910 (0.868–0.952)

当前单核加性模型下的估计,95% 区间按 Wald 近似。不代表已找到解释 91% 方差的 SNP。An estimate under the current single-kernel additive model; 95% interval by Wald approximation. It does not mean SNPs explaining 91% of the variance have been found.

待决定:Open: GLM 用几个主成分(与导师讨论,不以"加到 λ ≈ 1 为止"作为规则)。how many PCs the GLM should use (to discuss with the supervisor; "add PCs until λ ≈ 1" is not the rule).

先修Before you start
Step 0;知道 p 值是什么Step 0; know what a p-value is
学完能做After this page you can
解释 GLM 和 MLM 的差别,读懂 Manhattan 图和 QQ 图,说明 λ 能说明什么、不能说明什么Explain how the GLM and MLM differ, read a Manhattan and a QQ plot, and say what λ can and cannot tell you
预计用时Time
约 20 分钟about 20 minutes
01 · 两个模型01 · Two models

GLM 用主成分校正群体结构,MLM 用亲缘矩阵The GLM corrects for structure with principal components, the MLM with the kinship matrix

两个模型都对 P1 的 1,812,578 个 SNP 逐个检验,门槛同为 0.05 ÷ 1,812,578 = 2.76 × 10⁻⁸(−log₁₀ p = 7.56)。区别在于怎样扣除"亲缘相近的品系现蕾时间也相近"这部分。它们的结果分别为 A1 和 A2 提供候选 SNP。Both models test each of the 1,812,578 P1 SNPs, with the same bar of 0.05 ÷ 1,812,578 = 2.76 × 10⁻⁸ (−log₁₀ p = 7.56). They differ in how they remove the part where related lines have similar times to buds. Their results supply candidate SNPs to A1 and A2 respectively.

GLM · A1

y = μ + PC₁ … PC₁₀ + SNP + e

PLINK 2,10 个主成分(工作值),另跑了 0、3、5 个作对照。每次约 1.5 秒。PLINK 2 with 10 PCs (working value); 0, 3 and 5 PCs were also run for comparison. About 1.5 s per run.

MLM · A2

y = μ + SNP + g + e,  g ~ N(0, σ²ᵍ K)

GEMMA 0.98.5 -lmm 4,只含亲缘矩阵 K,不另加主成分,以 Wald 检验 p 为主。约 33 分钟。GEMMA 0.98.5 -lmm 4 with the kinship matrix K only, no extra PCs; the Wald p-value is the main one. About 33 minutes.

02 · 全基因组扫描02 · Genome scan

MLM 没有 SNP 过门槛,GLM 有 3,642 个No SNP clears the bar in the MLM; 3,642 do in the GLM

打开局部视图:Open a close-up:
也可以点击图上方框标出的区域。You can also click a boxed region on the plot.

QQ 图QQ plot

斜线是全部 p 值都服从均匀分布时的期望。MLM 的主体贴近斜线,只能说明中位统计量和主体没有明显的整体膨胀,不能说明尾部已经校准。The diagonal is the expectation if all p-values were uniform. The MLM body follows it, which shows only that the median and the bulk are not inflated overall; it does not show that the tail is calibrated.

λ 是怎么算的How λ is computed

把每个 p 值换成 1 个自由度的 χ² 统计量,取中位数,再除以 0.4549(理论中位数)。没有膨胀时 λ ≈ 1。Convert each p-value to a 1-df χ² statistic, take the median and divide by 0.4549, the theoretical median. Without inflation λ ≈ 1.

0.995MLM
2.54GLM,10 个主成分GLM, 10 PCs

多基因的遗传结构、LD 和样本量也会让 λ 升高(Yang et al. 2011),所以 λ 只能提示校准有问题,不能单凭它决定主成分个数。Polygenic architecture, LD and sample size also raise λ (Yang et al. 2011), so λ can only flag a calibration problem; it cannot decide the number of PCs on its own.

03 · 局部区域03 · Regions

三个已知开花相关基因附近的 ±100 kb± 100 kb around three known flowering-related genes

04 · GLM 的膨胀04 · GLM inflation

加到 10 个主成分,GLM 的 λ 仍是 2.54With 10 PCs the GLM's λ is still 2.54

主成分个数PCsλ过门槛的 SNPSNPs over the barp < 10⁻⁴ 的 SNPSNPs with p < 10⁻⁴峰(其中过门槛)Peaks (over the bar)
010.13191,279433,015368 (367)
34.1717,974110,981330 (286)
54.0819,211111,113328 (297)
102.543,64236,526216 (109)

峰指 p < 10⁻⁶ 的 SNP 按 ±250 kb 贪心分组,只作描述,不是 LD clumping。A peak is a greedy ±250 kb grouping of SNPs with p < 10⁻⁶, for description only, not LD clumping.

作对照,用随机 SNP 作加性背景的模拟性状,在同样 10 个主成分的 GLM 下 λ 为 1.2–1.5。真实 DTF1 是 2.54。这些模拟的随机背景没有在遗传力、LD 结构等条件上与 DTF1 匹配,所以两个 λ 的差别不能用来判断原因。膨胀的来源目前无法分解:结构校正不足、多基因效应和 LD 都可能有贡献,也不能逐个判断 GLM 多出的峰是真是假。A1 仍按事先规定用这个 GLM 做边际筛选,GLM 显著的 SNP 不称为可靠发现;所有互作都在 Step 2 用共同模型复检。For comparison, simulated traits with a background of random SNPs give λ of 1.2–1.5 under the same 10-PC GLM. The real DTF1 gives 2.54. The simulated random backgrounds were not matched to DTF1 in heritability, LD structure or similar conditions, so the difference between the two λ values cannot identify a cause. The source of the inflation cannot yet be separated: under-corrected structure, polygenic effects and LD may all contribute, and the GLM's extra peaks cannot be judged true or false one by one. A1 still uses this GLM for marginal screening as pre-specified; SNPs significant in the GLM are not called reliable findings, and every interaction is re-tested with the common model in Step 2.

05 · 与 AraGWAS 对照05 · Check against AraGWAS

两个探索性区域与 AraGWAS 的同一表型大体对应The two exploratory regions broadly match AraGWAS for the same trait

AraGWAS Catalog 的研究 703(doi:10.21958/gwas:703)对同一个 DTF1 做过 GWAS。它用的是 2,029 个品系的填补基因型和 AMM 模型,位点的 MAF 低到约 0.001。按 MAF 是否 ≥ 0.05 分开看:Study 703 in the AraGWAS Catalog (doi:10.21958/gwas:703) ran a GWAS on the same DTF1. It used imputed genotypes for 2,029 lines and the AMM model, with SNPs down to MAF ≈ 0.001. Split by MAF ≥ 0.05:

AraGWAS 位点AraGWAS SNPsSNP 数SNPs门槛 −log₁₀ pBar −log₁₀ p峰Peaks领头 SNP 的 MAF < 0.05Lead SNP MAF < 0.05
全部All10,374,3598.32114114
MAF ≥ 0.051,819,5087.561(5:18592588,DOG1)1 (5:18592588, DOG1)0
本项目的峰Our peak本项目 pOur pAraGWAS 同一 SNP(−log₁₀ p)AraGWAS, same SNP (−log₁₀ p)AraGWAS ±250 kb 内最好AraGWAS best within ±250 kb
MLM 5:18590501 · DOG13.38 × 10⁻⁸6.188.18,5:18592588,相距 2,087 bp8.18, 5:18592588, 2,087 bp away
MLM 5:3161477 · FLC 下游MLM 5:3161477 · downstream of FLC7.31 × 10⁻⁸7.31同一个 SNPsame SNP
GLM 5:3188327 · FLC 上游GLM 5:3188327 · upstream of FLC2.57 × 10⁻²⁵6.707.31,5:3161477,相距 26,850 bp7.31, 5:3161477, 26,850 bp away

MAF ≥ 0.05 的位点是全部位点的子集,但门槛随位点数重算后,得到的峰集合不是原来 114 个峰的子集:全部位点下的 114 个峰,领头 SNP 都是稀有变异;只留 MAF ≥ 0.05 时门槛放宽,才出现 1 个常见变异的峰。两份分析用的是同一个表型和相近的基因型来源,所以这只是大体对应,不是独立样本的重复验证。The MAF ≥ 0.05 SNPs are a subset of all SNPs, but once the bar is recomputed for the smaller set, the resulting peaks are not a subset of the original 114: all 114 peaks on the full set have rare-variant lead SNPs, and the single common-variant peak appears only because the bar loosens with MAF ≥ 0.05. Both analyses use the same trait and related genotype sources, so this is broad agreement, not replication in an independent sample.

06 · A1、A2 的候选06 · Candidates for A1 and A2

A1 的候选对是 A2 的约 360 倍A1 has about 360 times as many candidate pairs as A2

候选 SNP 只从 P2 面板(r² 0.2)中挑。严格版取过门槛的 SNP;筛选版取 p < 10⁻⁴ 的 SNP,按 r² 0.2、250 kb 做 clumping 后最多保留 200 个,不足不补。挑出的 SNP 两两配对后交给 Step 2,在那里先查四格和秩,再用共同模型检验。Candidate SNPs are taken only from the P2 panel (r² 0.2). The strict version keeps SNPs over the bar; the screened version keeps SNPs with p < 10⁻⁴, clumped at r² 0.2 and 250 kb, up to 200 and not topped up. The selected SNPs are paired and passed to Step 2, where the four-cell and rank checks come first and then the common-model test.

策略Strategy严格版:SNP → 对(四格合格)Strict: SNPs → pairs (four cells OK)筛选版:p < 10⁻⁴ → clumping → 保留Screened: p < 10⁻⁴ → clumped → kept筛选版对数(四格合格)Screened pairs (four cells OK)
A1(GLM,10 个主成分)A1 (GLM, 10 PCs)52 → 1,326 (1,290)879 → 841 → 20019,900 (19,443)
A2(MLM)0 → 011 → 11 → 1155 (54)

A2 严格版没有候选,工作流事先预见了这种情况,主分析用筛选版。A1 的候选来自膨胀的 GLM,可能含有较多与群体结构相关的信号;这是比较 A1 和 A2 时要解释的现象,没有事先删除。表中的对数还不是正式检验数,秩检查在 Step 2 进行。A2's strict version has no candidates; the workflow anticipated this, and the main analysis uses the screened version. A1's candidates come from the inflated GLM and may be rich in structure-related signal; this is something to explain when comparing A1 and A2, and it was not removed in advance. The pair counts are not yet the formal test counts; the rank check happens in Step 2.

SNP 遗传力 0.910 指什么What the SNP heritability of 0.910 means

这是 GEMMA 零模型的 REML 估计(pve = 0.910,标准误 0.021,Wald 近似 95% 区间 0.868–0.952),条件是当前 936 个品系、原尺度的 DTF1、中心化的 K 和只含截距的单核加性模型。它说明全基因组 SNP 合起来能解释大部分方差,但效应分散在很多位点上,所以没有单个 SNP 过门槛并不矛盾。它不是折外预测的 R²,也不是互作的贡献。This is GEMMA's REML estimate from the null model (pve = 0.910, standard error 0.021, Wald 95% interval 0.868–0.952), given the current 936 lines, DTF1 on its original scale, a centred K and a single-kernel additive model with an intercept only. It says the genome-wide SNPs together explain most of the variance, but the effects are spread over many loci, so having no single SNP over the bar is no contradiction. It is not an out-of-fold prediction R², nor the contribution of interactions.

07 · 常见错误07 · Common mistakes

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

1

按坐标大小判断上下游,方向写反了Upstream and downstream judged by coordinates, and got backwards

FLC 在负链上,5′ 端在高坐标一侧。MLM 峰 5:3161477 在低坐标一侧,所以是 FLC 的 3′ 下游约 11.9 kb,最初写成了"上游约 12 kb"。GLM 最强的 5:3188327 在高坐标一侧,才是上游约 8.9 kb。FLC is on the minus strand, so its 5′ end is at the high-coordinate side. The MLM peak 5:3161477 lies on the low side, about 11.9 kb downstream (3′) of FLC; it was first written as "about 12 kb upstream". The strongest GLM SNP, 5:3188327, is on the high side and is upstream by about 8.9 kb.

做法:上下游一律按基因所在的链判断。靠近 FLC 也不能说明这个位点通过 FLC 起作用。Practice: decide upstream or downstream from the gene's strand. Being near FLC does not show that a site acts through FLC.

2

把两组峰当成嵌套关系Two sets of peaks treated as nested

最初写成"去掉稀有变异后,AraGWAS 保留了 1 个峰、删掉了 113 个"。实际上全部位点下的 114 个峰都由稀有变异领头;MAF ≥ 0.05 时门槛放宽,才出现那 1 个常见变异的峰。It was first written that "after removing rare variants, AraGWAS keeps 1 peak and loses 113". In fact all 114 peaks on the full set are led by rare variants; the single common-variant peak appears only because the bar loosens with MAF ≥ 0.05.

做法:位点集合是子集,不等于峰集合也是子集。门槛随检验数变化时,两种口径的峰要分开报告,不能相减。Practice: a subset of SNPs does not give a subset of peaks. When the bar changes with the number of tests, report the peaks for each view separately; do not subtract them.

3

作图时的抽稀被读成数据特征Thinning for plotting read as a data feature

Manhattan 图为了加快绘图,把较大的 p 值抽稀,图上会在抽稀边界处出现一条密度分界线。这是作图方式造成的,不是数据本身的特征。本页的图同样做了抽稀,图注写明了边界。To draw quickly, the Manhattan plot thins the larger p-values, which leaves a visible density edge at the thinning cut-off. It comes from the plotting, not the data. The plot on this page is also thinned, and its caption states where.

08 · 练习08 · Exercises

自己算一次 λCompute λ yourself

下面的代码生成 100 万个均匀分布的 p 值,算出 λ;再把它们整体"放大",看 λ 怎么变。The code generates one million uniform p-values and computes λ, then inflates them and shows how λ changes.

python3 -c " import numpy as np from scipy import stats rng = np.random.default_rng(1) p = rng.uniform(size=1_000_000) chi2 = stats.chi2.isf(p, 1) print('lambda, null ', round(np.median(chi2) / 0.4549, 3)) print('lambda, chi2 x 2.5', round(np.median(chi2 * 2.5) / 0.4549, 3)) "

第一行应接近 1,第二行接近 2.5。真实数据里 λ 升高可能来自结构、多基因效应或 LD,单看 λ 分不出来。The first line should be close to 1 and the second close to 2.5. In real data a raised λ can come from structure, polygenic effects or LD, and λ alone cannot tell them apart.

FLC 在负链上,坐标 3,173,382–3,179,448。SNP 5:3161477 在它的哪一侧?FLC is on the minus strand at 3,173,382–3,179,448. Which side of it is SNP 5:3161477?

本页产物Your output for this page

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

用三到五句话说明:在 DTF1 上,为什么 GLM 有 3,642 个 SNP 过门槛而 MLM 一个也没有;这对 A1、A2 的候选意味着什么;为什么不能只看 λ 决定主成分个数。In three to five sentences: why does the GLM put 3,642 SNPs over the bar on DTF1 while the MLM puts none; what does that mean for the A1 and A2 candidates; and why can λ alone not decide the number of principal components?

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