文章：Mapping regulatory variants controlling gene expression in drought response and tolerance

文章整体思路：通过对转录组数据的分析挖掘与抗旱相关的基因。其中的一些数据处理及分析方法值得学习借鉴。

一、数据概览

样本：224 个植株的叶片

条件：3 种浇水处理，正常浇水 (WW) 、、停止浇水9天 (含水量70%，WS1)、停止浇水13天 (含水量58%，WS2)

设计目的：

1. 严格控制基因型变量，即对照组与实验组之间基因型完全相同；
2. 严格控制植株生长阶段。转录组会随着时间变化，处于不同的生长阶段的植株其转录组也存在较大差异。如 本研究中检测到的 eQTL 与授粉后 15 天玉米粒中鉴定的 eQTL 相比，只有约12％的 eQTL 重叠，表明基因表达是组织依赖的和动态的。 植株在干旱等 逆境条件下的生长速率会降低，干旱 13 天和正常生长13天植株所处的生长阶段可能是不同的。所以为了去除生长阶段对转录组的影响，作者使用纵向比对。虽然植株也在生长，但使用同一个植株比较干旱前和干旱后优于两个植株进行相同时间的不同处理。另一方面，三种状态下植株都处于营养生长的 2-3 阶段。从 WW 与 WS1 间隔9天，WS2 与 WS1 间隔 4 天来看，个人认为 WS1 与 WS2 间基因表达的差异更具有说明性。因为 1. 间隔时间短；2. WS1时植物已经处于干旱状态，生长速率明显低于WW到WS1之间的生长速率。综上 WS1 与 WS2 的比较中时间给转录组带来的差异更小。但也许要注意，WW 到 WS1 是出现干旱，而 WS1 到 WS2 是干旱增强，所以 WW 与 WS1 的比较仍然是有重要意义的。

转录组：627，其中 WW (209)、WS1 (208)、WS2 (210)。其中有表达量的 gene 约 30000 个，在不同处理下表达量差异 > 20% 的 gene (DEG) 数量为15369 个。

基因组：129w SNP (筛选条件：缺失率 < 0.6，MAF > 5%)。基因组数据是根据转录组数据 call 出的，所以没有转录本的地方不可能有SNP，SNP总量较少，只有129w。

二、数据处理

1. 不同处理下转录组数据的 PCA 降维结果间下图c ：
   PCA 作图用的是 DEG，增加了处理的组间差异，排除了在组间无差异或差异小的基因，排除了无关变量或混杂有随机变量的基因 (表达量差异较小的基因其表达量的变化难以去除随机因素的干扰)，使PCA结果更具有区分性。

PS：这一点可尝试应用于高光谱数据，将有区分度的波长提取出来进行分析。

2. 对 DEG 使用 kmeans 聚类算法 (k=3) 聚类。聚出了随着实验的进行gene上调 (6546，图d)、下调 (4569，图e)、无规律 (4254，图f) 三类。由于 ABA 对植物响应环境胁迫至关重要，所以作者着重研究了与 ABA 相关的基因，其中 DEG 总计有 56 个，31 种上调，14 种下调、11 种无规律。

4. 通过表达量 GWAS 找寻 QTL。若 QTL 中超过阈值的 SNP<3 则视为噪音，不计入统计。总共挖掘到 73573 个 eQTL (23,771 WW、22,945 WS1、26,857 WS2)，其中约 60% 的 QTL 内只有一个基因。
   作者将 QTL 内的 peaksnp 称为 leadsnp，将 leadsnp 在基因周围 (上下游 20kb 内) 的 eQTL 称为 local eQTL，把非周围的称之为 distant eQTL。下图用基因的起始位点作为纵坐标，eQTL 内 leadsnp 作为横坐标绘制散点图，可以发现对角线上的点十分浓密，说明调控基因的元件主要分布在基因周围，即主要为 local eQTL。
   
5. 作者将在三种处理下均检测到的 eQTL 称为静态 eQTL，在一组或两组中被检测到的 eQTL 称为动态 eQTL。其中 19566 个（26.6％）为静态 eQTL，5400 7个（73.4％）为动态 eQTL。

在 GWAS 中能检测出峰，必须要求表达量在样本间存在差异、基因型在样板间存在差异，基因型和表达量之间存在相关性。由于不同处理下基因型是不变的，变化的只可能是表型，所以对于 动态 eQTL 说明环境改变对基因表达量有影响且影响不是等比例施加在群体中各个样本上的，即不同样本对相同环境的响应不同。而静态 eQTL 则说明环境的变化不同个体的影响是成比例的，个人认为这类 eQTL 调控的基因可能是与基础生理代谢等相关。作者在后续的内容中也进行了佐证：local eQTL 中的 9168 个静态 eQTL 涉及到 3056 个基因，其中最显著的几个 eQTL 关联的基因都是与生长发育相关的。

根据作者先前将表达量划分为上升、下降、不确定三种，动态 eQTL 也可以再细分为三种：信号增强、信号下降、信号不确定。其中信号增强表明随着干旱压力的增加，表达量与 SNP 相关性增加，且表达量也随干旱压力的增加而不断提高，那么该位点大概率与抗旱相关。如果表达量是降低的那就不好说了，虽然可能是抑制因子，随干旱压力的增强而减弱，但由于表达量较低是随机误差占比增加，所以难以说明是否受到随机误差较大的干扰。信号下降也可能受到这一因素的影响。对于不确定的就不好统一分析了。

所以我认为在动态 eQTL 中 leadsnp 的 -log (p-value) 随干旱压力增加而逐渐变强且基因表达量升高的位点尤其值得关注。

5. 孟德尔随机分析
   作者对挖掘到的eQTL进行了生物学功能分析，发现找出来的许多eQTL可能与植物抗旱相关。但不得不说这种方式无法证明作者找到eQTL的有效性，读者更关心实验找到与抗旱相关eQTL的准确性。同时，由于很多基因的功能都是通过同源基因功能推测得到，具体的效应有待实验检验，仅通过这种程度的分析是远远不够的。为了从eQTL中筛选出与可信度更高的基因，作者进行了孟德尔随机化检验（MR test）筛选出了97个候选基因。

6. 实验验证 abh2 基因
   作者通过实验验证了 abh2 基因在植株抗旱过程中的作用。作者首先对 abh2 基因中的 SNP 位点进行了 GWAS 分析，发现其中与表达量显著相关的位点要么是同义突变，要么是非编码区突变，符合预期，因为蛋白结构的改变很难通过表达量 GWAS 挖掘出来。作者以 abh2 基因中的显著 SNP 位点合并绘制单倍型，发现此基因的 1 号单倍型的表达量显著高于其他单倍型，图b。通过基因编辑技术破坏基因制造突变株，发现突变株在干旱条件下的脱落酸含量显著降低 (30%)、叶片失水率显著下降、植株存活率显著提高。验证了 abh2 基因与玉米中的抗旱有关。