文献阅读---玉米干旱响应和耐受性基因表达的调控变异定位

文献阅读—玉米干旱响应和耐受性基因表达的调控变异定位

Mapping regulatory variants controlling gene expression in drought response and tolerance in maize

杂志:
Genome Biology
发表时间:
2020年7月6日

1.摘要背景

(1)基因表达
基因表达是细胞响应的决定因素。
基因表达中的自然变异将遗传变异和表型联系了起来。
(2)干旱胁迫
干旱是全球作物生产的主要环境威胁。因此,干旱胁迫下,控制基因表达的调控变异鉴定尤为重要。

2.研究背景(通常3-4段)

玉米是重要的作物,可以用于粮食、饲料和工业原料。
玉米也是重要的模式作物,具有丰富的遗传多样性和全基因组水平快速的连锁衰减。
干旱是重要的环境胁迫,作为复杂性状,解析困难。全基因组关联分析(GWAS)辅助检测基因和表型关联,但是因果基因和性状响应DNA变异鉴定困难。(GWAS的不足)
前人研究表明( cumulative evidence suggests that),遗传变异通过调控变异调节基因表达和基因变异中的重要功能来影响复杂性状。(后面是前人研究文献举例)

在许多模式生物中研究了基因表达变异的遗传基础,也称为数量性状基因座 (eQTL)。经典的亲本构建群体策略(classical strategy)通常定位区间大,而最近GWAS分析结合玉米材料的特异型(LD衰减快),可以将变异检测定位到单基因水平。另一方面,基因表达是时空依赖性的(spatio-temporal dependent),eQTL结构能够在胁迫和环境刺激时动态变化,基于此在人类组织研究中已经进行了一些探索,在植物中研究较少。

重要的是,基因表达改变是重要的分子表型,将基因型与分子功能和植物适应性联系起来,孟德尔随机化 (MR) 分析和转录组范围关联研究 (TWAS) 等策略来识别表达-性状关联,这些关联优化了基因表达变化与性状相关的候选基因。

最后总结本研究策略。(此部分MR和TWAS文献值得后续阅读)

3.研究材料和方法

3.1 材料种植和取样测序
玉米群体:
224个玉米材料

种植策略:
每个材料种植3个重复,幼苗2周龄开始干旱胁迫处理

干旱胁迫策略:
对两个重复栽培进行干旱处理(WS1、WS2),通过停水进行,同时保持一个重复以正常生长(WW)

取样策略:
WW(正常浇水),WS1(抑制浇水,相对含水量70%),WS2(相对含水量58%)
在干旱处理的第 9 天和第 13 天,分别收获了 WS1 和 WS2 的样品。
WW 样品与WS1在同一天采集 ;WS2 在一天中的相似时间点进行采样,以避免 4 天后基因表达的昼夜节律差异。

RNA测序样品采集:
从每种种质的3株植物中收获第二片叶子中间的 5 厘米叶子切片,汇集并在液氮中冷冻。
在提取 RNA 之前,将样品储存在- 80 °C。使用 TRIzol 试剂(Biotopped)分离总 RNA。

重复设置:
58 个生物重复被随机放置在不同的处理中,以估计实验程序的可重复性。
生物学重复的平均表达相关系数为 0.83,表明广泛采样和表达量化过程的合理重现性。

3.2 reads比对

低质量reads过滤和adapter序列的reads过滤;
过滤测序少于10 million reads的样本(测序未饱和,表达定量不准确)。

测序reads过滤如下:
224个材料的627个测序样本(WW, n = 209; WS1, n = 208; WS2, n = 210),以及58个随机选择的生物学重复样本。

比对软件:
HISAT2 ver.2.0.5,使用参数 (with parameters: --dta; --score-min L,-0.6,-0.6),
比对玉米基因组为:B73_RefGen_v4
(ftp://ftp.ensemblgenomes.org/pub/plants/release-38/fasta/zea_mays/dna/)

最后,仅保留单端读取的比对质量 (MQ) > 25 的比对和双端reads的唯一比对数据进行后续分析(这部分数据的比例是多少??)。

3.3 表达量计算和聚类分析
所有具有 MQ > 25 的比对reads都根据基因组位置排序。通过featureCounts计算每个基因,通过edgeR包识别不同种质间的差异表达基因。(根据基因组位置和基因名字排序的区别)
差异基因的鉴定
任何两种处理之间至少有20%的种质存在显着差异的基因(错误发现率(FDR)< 0.05)被认为是差异表达基因(DEGs)。

3.4 FPKM计算和批次效应剔除
StringTie ver.1.3.4d(参数:-G -e -B -o -A)计算FPKM,为消除实验技术引起的批次效应和系统变异,通过 R 中 LIMMA 包的 normalizeQuantiles 函数对每个样本的所有表达基因的 FPKM 进行归一化。

3.5 聚类
主成分分析是使用 R 中的 prcomp 函数进行的,设置 scale = TRUE。根据肘法算法,通过pheatmap包绘制表达式矩阵,kmeans_k = 3

3.6 eQTL计算转化FPKM
我们使用正态分位数变换(R中的qqnorm函数)对每个基因的FPKM值进行归一化。在至少 20% 的种质中,标准化 FPKM 高于 0.05 的那些基因被计算用于通过混合线性模型检测 eQTL。
MR因果分析为什么要使用qqnorm进行转换?因果分析不能使用正态分布的数据?

3.7 SNP检测和过滤
为了检测SNP,将3种处理的比对文件(MQ>25)合并到一个BAM文件进行SNP鉴定(GATK 4.0)。
第一,使用Picard 包 v1.115 MarkDuplicates (http://broadinstitute.github.io/picard/) 标记那些重复的读取,并使用 GATK 的 SplitNCigarReads 包过滤掉拼接映射读取。通过 GATK 的 HaplotypeCaller 包和 SAMtools 的 mpileup 包进行 SNP变异检测。

第二,将上面两个软件检测的SNP提取交集,并进行如下过滤获得可信SNP。
质量深度(QD)≥ 2.0 和 ReadPosRankSum ≥ − 8.0 和 FS ≤ 60.0 和 Qual ≥ meanQual。

第三,使用 GATK 的 BaseRecalibrator 包根据定义的可信SNP重新校准基础质量得分。

第四,通过 GATK 的 HaplotypeCaller 和 GenotypeGVCFs 包执行了最终的 SNP 变体调用。所有检测到的双等位基因 SNP 均用于进一步分析。

对于检测到的SNP,认为MAF小于0.05为错误检测。将minor allele替换为major allelel,将杂合子位点掩盖为缺失数据,获得了 1,899,223 个 SNP。

与368 个近交系的高质量1.25M SNP 数据集结合,共获得2.9M SNP,过滤MAF >= 0.05和max missing rate < 0.6,最后得到1,288,889个SNP用于后续关联分析。

3.8 eQTL的GWAS鉴定
使用混合线性模型计算每个表达基因的 eQTL 分析,其表达水平 (FPKM) 在整体样本的20%材料中表达高于0.05(WW 中的 29,614 个基因,WS1 中的 30,019 个基因,WS2 中的 30,272 个基因)。通过结合亲属系数、种群结构和隐藏的混杂因素,由 TASSEL(5.0 版)和 -mlm -mlmVarCompEst P3D -mlmCompressionLevel None 执行。

亲缘关系通过PHYLIP均匀选择81,612个SNP,群体结构信息使用 fastSTRUCTURE,混杂因子使用PEER推断,保留14个混杂因子放入混合模型中

显著阈值为:
Benjamini-Hochberg (BH) rejection thresholds were set for WW, P < 7.64E−7; WS1, P < 7.20E−7; and WS2, P < 8.36E−7 (FDR < 0.05)。

eQTL检测

三个连续步骤确定特定性状的 eQTL,
(1) 如果它们的物理距离 <5 kb,我们将所有重要的 SNP 分组为一个 eQTL; 
(2) 如果来自不同 eQTLs 的 SNP 在一个 LD 块内 (r2 >0.1),由 Haploview 计算,
参数为:-n -maxdistance 250 -minMAF 0.05 -hwcutoff 0 –dprime,将 eQTLs 合并,最重要的 SNP 被保留;
 (3) 如果在一个候选调控基因中检测到多个eQTL,则将它们合并,保留最显着的SNP。
 最后,最显着的 SNP(lead SNP)代表检测到的 eQTL(`该lead snp可以用来代替显著性,效应值等需要一个数值表示snp的结果`)。

3.9 层级网络推断
整套可用的表征玉米转录因子(www.grassius.org/grasstfdb.php 和 http://plntfdb.bio.uni-potsdam.de/v3.0/)用于鉴定某些 TF 家族。如前所述计算 TF 的层次结构高度:h = (O − I)/(O + I),其中 O 和 I 分别是通过调节关系检查的 TF 的出度和入度,如所述 (前人一篇science文章)。

3.10 孟德尔随机化分析
进行 MR 分析以揭示基因表达水平与表型之间的关系,如前所述(MR分析文章值得阅读)。

4.研究结果

Figure1 种植示意图和表达数据聚类
a 表型示意图
b 相对含水量取样点(每次测量三个生物学重复?这个生物学重复是如何获得的?推测是每个盆里面种植了很多水稻苗子,作为生物学重复,测量时来源不同的生物学个体即为生物学重复。)
c 使用转录组数据进行群体聚类
d 56个ABA相关基因聚类研究其表达模式
在这里插入图片描述

Figure2 eQTL整体展示
a 三个不同干旱条件下eQTL检测汇总,对角线为local eQTL,其他为distant eQTL;
b 不同处理中local 和distant eQTL的 marker R方概率密度,local eQTL面积更大,表明,可能作用更强;;
c 不同处理中划分动态和静态eQTL,将三个处理中均能检测到的eQTL定义为静态eQTL,否则为动态eQTL
在这里插入图片描述

Figure3 动态和静态eQTL中检测到的基因
图中a,b为生长发育相关的eQTL,在三个处理中均能检测到;
图中c,d为仅在胁迫下检测到的,表明这些eQTL是盐胁迫响应的;
图e为TF基因的调控关系(文献佐证),证明这类基因的重要性
在这里插入图片描述

Figure4 构建层级网络,寻找调控关系
盐胁迫基因表达调控,作者关注Dynamic eQTL基因,所以关注TF基因和target gene(etrait)。
根据文献,作者选择44个已知基因的etrait(分散不同的生物学过程),共检测45个TF基因。
其中一些基因的互作在拟南芥的ChIP研究中有研究记录。
在这里插入图片描述
Figure5 通过孟德尔随机化分析,缩小候选基因集
图a 基于测量的群体表型和基因表达量,进行MR分析,97个基因显著。
51个基因正调控,46个基因负调控。其中该可以组之前研究的ZmVPP1基因也被检测到了。(Nature genetics)
图b 9个候选基因单倍型表达分析和表型展示,并用于后续验证工作
在这里插入图片描述

Fgirue6
基于MR分析,挖掘到的候选基因中,共得到9个候选基因,Figure5所示。
将上面9个候选基因中Zm00001d051554 (abh2) 基因敲除验证。该基因缺失增加植物耐旱性。


问题如下:

  1. survival rate 是如何调查的,实际上Figure6中植株并未死亡,我推测是根据植物的枯萎百分比人工调查。
  2. 玉米在V2-V3期间开始处理,这个生育阶段是怎么记录的,去查一下。
  3. RNAseq的unique mapped reads 被提取用于后续表达分析。
  4. 单倍型分析中的数字是什么含义,这个不能解决。
  5. 转基因和基因敲掉的操作步骤了解一下
  6. 孟德尔随机化分析,TWAS
  7. qqnorm的操作问题
  8. eQTL的检测问题,三个显著SNP之间的距离怎么划分,是根据群体的LD,0.2为1.6kb,还是0.1的~20kb(方法中使用5kb合并形成eQTL
  9. eQTL中lead snp 上下20kb的etrait为local ,否则为distant

参考:
https://genomebiology.biomedcentral.com/articles/10.1186/s13059-020-02069-1#Sec2

  • 2
    点赞
  • 8
    收藏
    觉得还不错? 一键收藏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值