责编 | 王一
BSA是Bulk Segregation Analysis 的英文首字母缩写,翻译过来就是分离群体分组分析或者是混池分离群体分析。BSA基于遗传分离群体的极端表型构建混池,通过混池中变异的频率差异筛选与表型关联的分子标记,是一种简便高效的性状定位方法。在下一代测序 (NGS,next-generation sequencing) 的背景下,BSA因其快速的定位和超高的性价比逐渐崭露头角并受到遗传和育种科研人员的广泛欢迎。2021年12月21日,The Plant Journal在线刊登了NGS时代的BSA:十年进展综述 (Bulk segregation analysis in NGS era: a review for its teenage years) 。中国农业大学博士许语辉为论文通讯作者,新疆艾德森生物科技有限公司总裁李志强为论文的第一作者。下文为全文翻译,以飨读者。因文章正在校稿过程中,翻译与first online版本会略有差异,更建议读者下载英文版阅读。
1. 前言
Bulk Segregation Analysis (BSA) 最早由加州大学戴维斯分校蔬菜作物系 (Department of Vegetable Crops, University of California) 的Michelmore等人提出 (Michelmore et al., 1991) ,用来鉴定F 2 群体中和Dm5/8 (生菜霜霉病抗性位点) 连锁的分子标记。通过对抗、感病两种极端表型分组,每组内通常14-20个单株的个体混池,再基于分子标记 (AFLP,RAPD,SCAR,SSR等) 观察混池中的条带多态性情况。判别标记和性状连锁与否的基本原理是:如果某个标记在同一个分组内所有混池中的条带有差异,且两极端性状分组混池的条带一致,则说明该标记和目标序列区段不连锁;如果一个混池中条带无差异,且在两极端混池条带呈现多态性,则说明该标记和目标序列区段连锁。
和对单株进行基因分型 (genotyping) 相比,混池的方法能极大地减少工作量和试剂耗材成本,提升连锁标记筛选效率,在基于遗传分离群体的性状基因定位上发挥了重要作用。之后不同研究者又利用DH、RIL等群体进行分子标记和目的基因之间连锁关系的研究,对不同物种的不同性状进行了基因的定位 (Quarrie et al., 1999; LIN et al., 2006; Venuprasad et al., 2009) ,拓展了BSA的应用范围。此外,BSA定位对性状的要求也逐渐从质量性状过渡到有主效基因控制的数量性状 (比如抗病性状) 上来 (Yang et al., 2004; Hyten et al., 2009; Shoba et al., 2012) 。针对复杂的数量性状进行BSA定位,也有研究者做了尝试 (Zhang et al., 2009; Sun et al., 2010; Salunkhe et al., 2011) 。
在NGS (Next Generation Sequencing) 测序流行起来之前,基于遗传分离群体进行QTL定位还是以遗传图谱构建的QTL定位为主。随着NGS测序的成本不断走低,加上三代测序和Hi-C技术的使用 (Belton et al., 2012) ,作物、果树、花卉、林木、水产、昆虫物种的参考基因组 (reference) 序列逐步被释放,这是基于NGS测序的BSA分析的第一个有利因素。此处需要声明的是,虽然无参考基因组物种进行BSA研究的方法有文献报道 (Nordström et al., 2013) ,但是方法复杂,需要研究者有较强生信的基础,很难推广到不同物种中。另外只有小部分的无参考基因组物种BSA研究有文献报道 (Ye et al., 2016) ,但是只能得到大量的跟性状关联的SNP标记,无法知晓染色体上的位置和区间,对后续研究和验证非常不利。第二个有利因素是NGS使高通量、全基因组范围内开发分子标记 (SNP,InDel) 成为可能。BSA和性状关联就是基于全基因组范围的SNP或者InDel标记。目前,微生物酵母 (Ehrenreich et al. 2010; Parts et al. 2011; Liti and Louis, 2012) ,植物 (Ban et al., 2020; Li et al., 2020; Wang et al., 2020) 和动物 (Jiang et al., 2019; Yang et al. 2019) 中均有BSA定位性状的报道。本书后文将以BSA的适用群体类型、算法和如何提升BSA定位效果作为重心,详细描述最近十年来基于NGS的BSA的进展。
图1 BSA的进展史
2. 基于不同测序策略的BSA
基于NGS的BSA需要获得全基因组范围内的SNP标记。当前许多测序技术均可以进行SNP detection。本部分就从测序技术层面介绍BSA。
2.1基于简化基因组测序的BSA
简化基因组技术发展较早,应用也较成熟。简化基因组测序“化石”级的方法有CRoPS (Altshuler et al., 2000) ,RRLs (Van et al., 2008) ,这两种方法建库简单粗暴,现在已经不再使用了。后来陆续出现改良的简化基因组测序方法,比如RAD法 (Davey et al., 2010) ,GBS法 (Elshire et al., 2011) ,2b-RAD法 (Wang et al., 2012) ,ddRAD法 (Peterson et al., 2012) ,SLAF法 (Sun et al., 2013) 等。基于简化基因组开发的SNP能进行BSA研究 (Zhang et al., 2016; Zhao et al., 2016; Watanabe et al.., 2017; Feng et al., 2018; Li et al., 2018) ,关联性状的算法以SNP (或者多态性标签) 的比例,SNP (或者多态性标签) 的密度等较简单的算法为主,下文我们会详细介绍其中的某些算法。基于简化基因组测序的BSA物种往往是参考基因组较大的物种。
2.2基于重测序的BSA
重测序建库与简化基因组测序建库不同之处在于,重测序采用超声波打断基因组DNA (covaris) 的方式建库。基于重测序BSA分析最早是由Schneeberger等于2009 年开发出 SHOREmap 的重测序BSA分析方法 (Scheeberger et al., 2009) ,该方法将 500 个来自于 EMS 突变体和正常植株杂交的 F 2 个体混合成一个样品进行深度重测序,定位到了EMS 诱变的叶色基因突变。此方法又叫Mapping by Sequencing 法 (MBS) 。基于重测序进行BSA的其他方法也陆续发展起来,比如欧氏距离法 (Euclidean Distance,ED算法)(Hill et al., 2013) ,QTL-seq法 (Takagi et al., 2013a) ,InDel-index法 (Singh et al., 2017) ,G' value法 (Magwene et al., 2011) 和GradedPool-seq法 (Wang et al., 2019) 等。参考基因组中等或者较小的物种往往会进行重测序。
2.3基于转录组测序的BSR
BSR是将BSA与RNA-seq结合起来的分析方法,与重测序BSA不同的是,在分离群体中选择极端性状的个体构建两个池,提取两个池的总RNA,进行转录组测序,根据混池个数和物种的基因组大小确定测序的数据量。基于转录组数据同样开发SNP标记,用经典贝叶斯算法找到与突变基因共分离的SNP位点和区间,最后预测候选基因 (Liu et al., 2012) 。由于贝叶斯算法需要花费较长时间和较多的计算资源,仍然可以用ED或者SNP-index的算法进行BSR的性状关联。另外,BSR还能反映亲本 (混池) 的基因表达量信息,可以根据候选区段内相关基因的表达差异筛选候选基因,这个特性使得BSR对抗逆性状 (诱导型表型) 而言更加适用。BSR分析一般需要参考基因组,所以分析结果会受到参考基因组组装质量的影响,这点同重测序BSA一致。由于基因只占基因组的1-3%,BSR也更适用于基因组较大的物种,比如麦类物种。
图2 BSA分析的四个模块
3不同群体类型和BSA
3.1可进行BSA定位的遗传分离群体类型
一般来说,BSA性状定位首先需要构建一个遗传分离群体。遗传分离群体包括F 1 、F 2 或F 2:3 衍生群体、回交BC (back crossing) 群体、重组自交系RIL (recombined inbred lines) 群体、近等基因系NIL (near isogenic lines) 群体、双单倍体DH (doubled haploid) 群体以及多亲本的NAM群体,MAGIC (Multiparent Advanced Generation Intercross) 群体,ROAM (Random-Open-Parent Association Mapping) 群体等。另外,还包括一些次级分离群体,比如染色体片段代换系 (Chromosome Segment Substitution Lines,CSSL) ,剩余杂合系 (Residual Heterozygous Line,RHL) 等。
遗传群体类型丰富多样,目前来看,进行BSA研究的群体为双亲构建的遗传群体。常规能构建遗传图谱的群体,比如F 1 (Dai et al., 2018; Guan et al., 2019) ,F 2 (Klein et al., 2018; Feng et al., 2020) ,BC (Lalli et al., 2008; Khumto et al., 2018; Luo et al., 2018) ,DH系 (Miao et al., 2019) ,RIL (Pandey et al., 2017; Wang et al., 2018) 群体都可以进行BSA分析。那么除了以上典型的遗传群体之外,其他类型群体是否也可以进行BSA分析呢?答案是肯定的,BSA的关联算法并不关注基因/标记频率在整个群体中的分布,只是关注基因/标记在极端混池中的频率差异,只要位于阈值以上的区间都可以认为是潜在的关联区间。所以理论上只要是一个双亲的分离群体,均可以进行BSA分析。BSA在F 3 群体 (或F2:3家系群体)(Park et al., 2019) 和F 4 群体 (Negi et al., 2000) 中均有应用。
需要注意的是,对于某些遗传背景高度一致,目标位点有差异的群体类型,比如NIL群体、RHL群体或者EMS诱变的突变体和野生型亲本构建的分离群体等,不能构建全基因组的遗传图谱,最好就是通过BSA的策略定位,而且这类群体定位的结果往往非常理想。
3.2不同遗传分离群体的BSA定位结果比较
BSA的QTL检测效率和目标性状的遗传力、遗传分离群体数量大小和NGS测序的深度都有关系 (Guo et al., 2016) 。本段讨论群体类型对BSA定位结果的影响。在相同的群体大小和测序指标 (DNA测序,测序数据量) 的前提下,定位的效果决定于极端混池材料携带的重组事件的数量。而重组事件的数量又跟混池的规模和分离群体类型相关 (Zou et al., 2016) 。如果目标性状、群体数量、混池规模、测序方法都一致,那到底哪种群体类型做BSA定位效果最好呢?很早就有文献研究不同群体类型和重组事件的关系,Paran 等 (Paran et al., 1995) 在研究番茄性状遗传定位的时候发现,RIL群体中连锁标记之间的重组是F 2 群体的两倍。这一点也很好理解,因为从F 2 群体到RIL群体,每自交一次,杂合区域重组的概率就高一次,逐渐积累到RIL群体,就会产生比F 2 群体更多的重组交换。F 2 群体和BC群体之间的比较也有研究报道。Mackay等 (Mackay et al., 2000) 比较了相同亲本构建的F 2 群体和BC群体BSA定位的适用性。最终发现,在重组率一定的前提下,BC群体检测QTL能力稍稍优于F 2 群体。另外,BC群体的遗传背景更偏向于轮回亲本,BSA定位时,背景噪音会更小一些。但在实际构建群体时,F 2 只需要套袋自交,而BC则需要人工去雄后辅助杂交,因此,从群体构建工作量的角度,F 2 群体构建更容易,且能短时间产生大量的子代个体。
为了能更直接地展示不同类型群体的QTL检测效力,计算模拟实验结果表明:相比F 2 群体而言,F 7 具有更高的检测能力。F 2 群体中,当p = 0.15 (每个混池选择的个体数占群体个体总数的比例) 和n = 20 (测序深度20X,减少计算资源) 时可以检测贡献约为10% (Qp = 0.1) 的QTL (Takagi, 2014) 。
3.3自然群体/多家系群体BSA
遗传分离群体基于双亲 (bi-parent) 杂交而成,其遗传背景相对简单。而自然群体和多家系群体遗传背景复杂,进行BSA分析会面临一些挑战,其潜在困难和原因包括:
1)自然群体往往存在群体结构,会引入假阳性到分析中。TASSEL (Bradbury et al., 2007) ,EMMAX (Kang et al., 2010) 这些GWAS的分析方法会利用协变量对群体结构进行校正。这样的校正必须基于对个体进行基因分型,而基于混池测序的 BSA无法对群体结构进行校正。
2)自然群体遗传多样性太高,那么在DNA池中必然存在多种基因型,在很小的一段DNA区域存在3种或3种以上的等位基因基因型完全是可能的。DNA混池的高多样性,会导致SNP检测和基因型频率计算可靠性下降。
3)自然群体的QTL效应可能比较弱,影响了QTL被检测到的概率。
虽然自然群体的BSA存在以上的一些困难和风险,但是这并不妨碍一些特定的方法能能推进此类研究。Pool-GWAS (Bastide et al., 2013) 和XP-GWAS (Yang et al., 2015) 就是其中两种代表性的方法。
Pool-GWAS适合于多家系的材料。首先进行种质资源的收集,然后开放性杂交,然后在F 1 代中选择极端表型个体构建混池。在判别SNP与表型关联性上,Pool-GWAS采用分层卡方检验 (Cochran-Mantel-Haenszel test) 法检验 (Kofler et al., 2011) 。首先,依据混池的表型和基因型构建2×2的列联表。每个SNP位点的变异均为独立名义变量,每个等位基因的数量也为独立变量。然后通过分层卡方检验法检测。最后利用FDR校正多重检验值,阈值确定时控制FDR=0.05。该方法需要二等位的SNP,且由于分层卡方检验对测序的深度有要求,因此更适合多家系且亲本较少的群体。从最近5年发表的文献来看,使用该方法的研究并不常见。
XP-GWAS首先需要一个自然群体,然后依据感兴趣的表型进行极端混池的构建,这点类似于遗传分离群体。信息分析上,首先使用R包中的“GLM”函数对每个variant进行广义线性模型分析。对于某个表型的混池而言,与参考基因组一致位点和对应的等位基因位点的reads数均建模为二项式随机变量。将二项式成功概率的对数建模为表型混池数的线性函数。然后计算成功概率和表型混池之间无关联的无效模型似然比检验统计量。引入前人提出的通货膨胀系数λ并用R包计算该值 (https://cran.r-project.org/web/packages/gap/index) 。似然比检验统计量除以λ后得出P值,最后,FDR在5%时的P值作为阈值,阈值线以上的SNP即为和表型显著关联的SNP。
XP-GWAS的检测效力受表型鉴定的准确性、混池大小、选择强度、标记密度、测序深度和连锁不平衡 (linkage disequilibrium,LD) 程度等因素的影响。在选择强度 (5%) 不变的前提下,增加混池规模能提升小效应QTL的检出率。同时,作者建议测序的深度不低于混池的规模,比如100+100的混池规模,测序深度也至少应为100X+100X。如此深度下,标记密度往往不是XP-GWAS的短板,反而是LD对定位的分辨率有较大的影响。XP-GWAS的缺点与BSA类似,都是一个性状需要对应一种混池设计,而且无法估计SNP的遗传效应。
Yang等 (Yang et al., 2015) 利用该方法作者对玉米的每穗籽粒行数 (kernel row number,KRN) 性状进行了定位,并检测到145个性状关联位点 (trait-associated variants,TAVs) 。定位的分辨率与常规GWAS类似,其中17%的位点和已经报道的位点重合,35%的位点和未发表的位点重合。验证了该方法的有效性。目前该方法在水稻 (Xiao et al., 2017) ,印度大麻 (Welling et al., 2020) 等物种上应用。
4. BSA关联分析的算法
BSA一路发展到今天,其关联算法也逐渐演变和分化。根据不同的遗传群体材料和遗传设计,有针对EMS诱变材料MutMap法 (Abe et al., 2012) ,无亲本只有两个极端混池的欧氏距离法 (ED)(Hill et al., 2013) ,基于三个亲本 (一个共有亲本) 构建两个分离群体,每个分离群体都筛选两个极端混池后整合分析的multiple QTL-seq法 (Srivastava et al., 2017) ,基于一个大的F 2 群体并依据表型分组 (三组以上) 的GradedPool-seq法 (Wang et al., 2019) 和较复杂的遗传设计QTG-seq法 (Zhang et al., 2019a) 等;根据数学算法不同,除了以上的SNP-index法、ED法外,还有G' value法 (Magwene et al., 2011) 、贝叶斯 (Liu et al., 2012) 等算法。本部分就从遗传设计和数学算法两个角度切入,尽可能全面的描述BSA中不同的关联算法,最后对不同算法效果作比较。
4.1不同的遗传群体材料和遗传设计
4.1.1 亲本为自交系群体类型 (F2,DH系,RIL等)
亲本自交系类物种典型就是水稻、玉米等等作物。QTL-seq法 (Takagi et al., 2013a) 即是以水稻为研究模型的BSA经典算法。该算法引入了SNP-index的概念,其基本原理如下:SNP-index分别计算两个混池的基因型频率差异,然后用Δ(SNP-index)统计混池之间基因型频率的显著差异。标记SNP与性状关联度越强,Δ(SNP-index)越接近于1。
计算方法如下:
Maa表示aa池来源于母本的深度;
Paa表示aa池来源于父本的深度;
Mab表示ab池来源于母本的深度;
Pab表示ab池来源于父本的深度;
SNP-index(ab)=Mab/(Pab+Mab);
SNP-index(aa)=Maa/(Paa+Maa);
Δ(SNP-index)=SNP-index(aa)-SNP-index(ab)。
由于重测序能获得海量的SNP标记,为了消除假阳性的位点,根据标记在基因组上的位置 (一般需要染色体级别的参考基因组) ,可对同一条染色体上标记的ΔSNP-index值进行拟合,并根据关联阈值,选择阈值以上的区域作为与性状相关的区域。
QTL-seq基于SNP计算频率,而NGS测序除了检测SNP,还能检测InDel (1-5 bp) ,理论上也能基于InDel进行InDel-index计算。2017年Singh等 (Singh et al., 2107) 提出了InDel-seq的方法并成功利用InDel-index法进行了木豆的枯萎病抗性和不育花叶病抗性的QTL定位。该方法在番茄抗叶霉病的QTL定位中得到了应用 (Liu et al., 2019) 。
除了以上经典的基于SNP-index的QTLseq法,其他算法我们将在下文的数学算法内容中专门讨论。
4.1.2 EMS诱变类材料构建的群体
在Mutmap (Abe et al., 2012) 策略出来之前,Austin等报道了一种名为Next Generation Mapping (NGM) 的定位方法 (Austin et al., 2011) ,该方法利用EMS诱变的隐性突变体 (Columbia型拟南芥) 和比对系 (Landsberg型拟南芥,参考基因组已经公布的品种) 构建F 2 分离群体,F 2 群体中隐性性状个体构建混池后测序。测序后的数据与比对系对应的参考基因组比对和变异检测 (SNP calling) 。通过全基因组范围内SNP的频率判别非重组区域和目标位点区域。又基于修改过的chastity statistic统计量Termed discordant chastity (ChD) 判别真正的目标位点区域达到定位的目的。作者用该方法成功定位到拟南芥中3个和细胞壁合成相关的基因,显示出该方法的有效性。但是也有几点不足之处:1)该方法为了省去F 2 群体中非突变体亲本的测序成本,在设计之初就要求突变体和对应该物种已经发表了的参考基因组的品种进行杂交,这就增加了实验的沟通成本和时间成本;2)该方法形成的分析流程 (http://bar.utoronto.ca/NGM/) 位于网页上,对基因组较小的物种比如拟南芥而言实现起来较为容易,但是对于某些基因组较大的物种 (比如玉米) 而言,网页版分析流程不具备较大数据量计算的能力。3)该文发表时主要面向拟南芥 (只支持小基因组物种拟南芥,果蝇,秀丽隐杆线虫和酵母) ,对其他物种的支持度不够。
之后有更具普适性的方法被开发出来,即Mutmap法 (Abe et al., 2012; Takagi et al., 2015) ,Mutmap+法 (Fekih et al., 2013) ,Mutmap-Gap法 (Takagi et al., 2013b) 。下面我们就重点介绍这三种算法。
Mutmap:
Mutmap的材料是由EMS诱变的可育突变体和其与对应的野生型杂交得到F 1 ,再自交得到的F 2 群体。F 2 分离群体会产生突变型表型和野生型表型。突变表型为隐性性状的话,只需要取野生型亲本样品和F 2 代中突变体表型构建混池。同样计算混池中的SNP-index值,
SNP-index(Mut)=ρx/(ρX+ρx)
其中,Mut为子代的突变池,ρX和ρx分别为野生型亲本以及突变体的等位基因在突变池中出现的SNP的测序深度。只有SNP-index=1的位点才可能和表型关联。因此,Mutmap法实际上是QTLseq的一个极端的特例。
Mutmap+:
EMS诱变也会产生大量不育或者早期发育致死类型的突变,无法进行杂交而限制了MutMap的应用,而通过MutMap+策略则可以解决这个问题 (Fekih et al., 2013) 。其基本原理如下:用EMS处理刚受精后的种子胚胎,发育成熟后种植形成M1代植株。大多数M1代植株突变位点为杂合位点。M1自交获得M2种子。M2代会产生分离,选择M2中野生型:突变型=3:1的群体作为研究的群体,在符合上述条件的M2中随机选择多个M2的后代,并测定表型。对于发生分离的感兴趣的表型,2/3的野生型表型姊妹系目标区域为杂合。这些野生型表现的M2代自交后结的种子单株收获后获得M3。超过80个M3代种子种下用于表型观察,M3代突变型表型和野生型表型两个混池。计算两个混池的SNP-index,SNP-index=1时有两个可能,一个可能就是突变性状的位点,另一个可能就是和性状不相关,但是却固定在M2的纯合区域最后又遗传到所有的M3植株中。通过计算ΔSNP-index即可过滤掉假阳性区域,最终只留下SNP-index=1且与目标性状关联的SNP。
Mutmap-Gap:
当我们所用的野生型亲本材料是该物种参考基因组组装的品种时,通过Mutmap策略能定位到我们的目标基因。可是如果我们用的野生型材料和参考基因组材料差异较大时,就可能出现这种情况:SNP calling时需要比对参考基因组,如果ΔSNP-index定位区域内不存在跟突变性状相关的基因,那么该区域可能存在结构变异 (SV) 、大的InDel或者存在缺失变异 (PAV) ,那么用上述的MutMap就无法找出该基因及其序列。MutMap-Gap法便在MutMap的基础上多了一个de novo组装的过程 (Takagi et al., 2013b) 。De novo组装的具体过程是,首先收集野生型亲本un-mapped的reads和位于这个ΔSNP-index peak区域的野生型亲本mapped reads进行de novo组装,然后将组装的序列与之前的参考基因组重新mapping,最后再次计算SNP-index值并筛选SNP-index=1的位点和对应的序列。利用MutMap-Gap方法,该课题组成功鉴定了水稻Hitomebore品系的耐盐基因Pii上的一个SNP,而该基因Pii在日本晴参考基因组上是不存在的 (Takagi et al., 2013b) 。
由于本文研究的是水稻这个模式物种,推广到其他物种时候往往会出现的问题。笔者总结如下:重测序采用的是Illumina短片段测序,文库大小一般是是350 bp,而且不管是亲本还是混池,测序的深度都不高,一般在几十X (depth) ,会导致de novo组装完成的contig片段非常碎;对于参考基因组组装质量不高的物种 (尤其是杂合度较高的物种,比如林木花卉水产) 而言,Un-mapped的reads来源不明,组装后会引入较多的错误。针对这种情况,笔者提出除了构建BAC文库之外的三种解决方案:
1)目前很多物种组装了不止一套参考基因组,比如水稻就有日本晴 (Kawahara et al., 2013) 、93-11 (Lang et al., 2020) 、明恢63 (Zhang et al., 2016) Zhenshan97 (Zhang et al., 2016) 、R498 (Du et al., 2017) 、Basmati (Choi et al., 2020) 、33种不同品种 (Qin et al., 2021) 等等数十种不同籼粳品种的参考基因组。玉米也释放了B73 (Schnable et al., 2009) 、Mo17 (Sun et al., 2018) 、W22 (Springer et al., 2018) 、黄早四 (Li et al., 2019) 和K0326Y (Li et al., 2020) 等自交系的参考基因组。针对这些物种而言,我们可以先利用其中一套组装较完整的参考基因组分析,如果达不到目的就换一个参考基因组重新分析,有可能能达到意想不到的效果。
2)基于定位区间的序列信息设计探针,特异性地捕获 (capture) 该区域的序列,并结合三代 (PacBio,Nanopore) 测序,之后进行de novo组装是解决这类问题的潜在途径。关于捕获测序后de novo组装也有相关文献报道 (Steuernagel et al., 2016; Arora et al., 2019) 。
3)QTL定位完成后,对亲本进行Nanopore或者PacBio重测序,通过长读长覆盖BSA定位的区域,再进行变异检测或者进行de novo组装,识别该区间的大片段插入缺失和结构变异。比如Soyk等 (Soyk et al., 2019) 利用QTL-seq的方法发现在1号与3号染色体上存在2个控制番茄分枝的主效QTLs。接着基因组测序覆盖率的变化显示亲本fla.8924在EJ2位点周围有一个83 kb的缺口。之后使用ONT长读长测序在fla.8924中检测到相同的结构变异,证实了在这个品系中相同位置处有一个长约 83 Kb 的重复片段。
4.1.3 亲本高杂合类群体 (F1)
上文的遗传群体类型中我们介绍过F 1 类型的群体,这类群体主要在林木、果树、花卉和水产物种中。该类群体的特点是双亲均为高杂合个体,F 1 代即出现性状分离。最早开发出的 MMAPPR 法,利用斑马鱼构建的F 1 代分离群体,构建表型的极端混池,基于转录组数据检测SNP,利用欧氏距离法 (ED) 计算混池之间SNP的距离差异,成功定位到了斑马鱼两个心血管突变体基因 (Hill et al., 2013) 。
ED算法,是利用测序数据寻找混池间存在显著差异标记,并以此评估与性状关联区域的方法。理论上,BSA项目构建的两个混池间除了目标性状相关位点存在差异,其他位点均趋向于一致,因此非目标位点的ED值应趋向于0。ED方法的计算公式如下所示,ED值越大表明该标记在两混池间的差异越大。
其中:
A mut 为A碱基在突变混池中的频率,A wt 为A碱基在野生型混池中的频率;
C mut 为C碱基在突变混池中的频率,C wt 为C碱基在野生型混池中的频率;
G mut 为G碱基在突变混池中的频率,G wt 为G碱基在野生型混池中的频率;
T mut 为T碱基在突变混池中的频率,T wt 为T碱基在野生型混池中的频率。
在分析时,检测两混池基因型之间的SNP位点,统计各个位点在不同混池中的深度,并计算每个位点ED值。如果某些亲本遗传差异过大,对原始ED值进行乘方处理,然后采用局部线性回归的LOESS方法对ED值进行拟合。
最近,欧氏距离法最近又有新的发展,比如Shen 等 (Shen et al., 2019) 利用Nadaraya-Watson核回归作为欧氏距离smoothing的函数。其计算公式为:
这个smoothing处理源于Magwene等 (Magwene et al., 2011) 介绍的方法。
欧氏距离法解决了分离群体的亲本丢失或者死亡的问题,只需要两个极端混池即可开展定位的工作。这里也有例外,无亲本进行SNP-index计算的处理方法也有报道 (Clevenger et al., 2018) ,但是太小众,笔者在此就不展开了。
F 1 群体BSA分析除了以上的ED法,SNP-index的算法同样适用 (需要双亲测序数据) 。Xue等 (Xue et al., 2017) 在利用QTL-seq方法定位梨皮颜色性状QTL时,便用到了一种修改版的SNP-index法。对于杂合型亲本而言,相邻SNP位点的Δ(SNP-index)值可能发生正值与负值交替出现的情况,当对一个滑窗内的Δ(SNP-index)取平均值时,候选区域的值趋向于0,这种情况会干扰候选区域鉴定。所以该研究以Δ(SNP-index)的绝对值|Δ(SNP-index)|作为指标,发现其在杂合亲本构建的F 1 群体中具有较好的基因定位效果。
除了|Δ(SNP-index)|外,还有其他F 1 分离群体进行BSA研究的方法。笔者将在下文的数学算法部分继续揭晓。
4.1.4 特殊遗传设计
以上的是常规的遗传分离群体,但也有更复杂的遗传群体设计,下面我们就看看以下三种特殊遗传设计的方法。
multiple QTL-seq法
multiple QTL-seq (Srivastava et al., 2017) 中的multiple指的是有两个或者两个以上的分离群体,每个分离群体都单独进行QTLseq分析,然后整合结果。分离群体构建一般会共用一个亲本,其他两个亲本在该性状上和共有亲本有显著差异。同样基于分离群体的表型差异分别构建两个群体的极端混池。每个混池单独进行QTL-seq分析,然后根据定位结果取交集,以此来缩小定位区间。由于用了多个亲本,该方法还能找到新的QTL位点。
GradedPool-seq法
GradedPool-seq法由韩斌和黄学辉课题组提出 (Wang et al., 2019) 。该方法的基本原理如下:性状有较大差异的双亲组配,形成一个较大规模的F 2 群体 (最好数千株) ,依据F 2 分离群体的表型分成三组 (也可以三组以上) ,分别为“高值组”,“中值组”,“低值组”,然后进行个体的DNA提取和等量混合,形成三个混池 (或者三个以上的混池) 。这里面的“高值组”和“低值组”和普通的BSA一样,不同的是增加了一个“中值组”。与性状连锁的SNP在三个混池中会呈现出不同的频率分布,基于Ridit算法的三个混池之间SNP差异显著性比较,获得每个SNP的P值,为了消除背景噪音的影响,通过滑窗法 (slide window) 处理背景噪音。GradedPool-seq最终实现了水稻400Kb左右的定位分辨率。Liu等 (Liu et al., 2020) Gradedpool-seq法进行了水稻抽穗期表型的定位。共定位到2个主效QTL位点,分别位于第6号和第8号染色体上,这两个位点分别包含了已经克隆的Hd1基因 (Yano et al., 2000) 和Ghd8基因 (Yan et al., 2011) 。但是,GradedPool-seq并未声明对除了F 2 群体之外,其他群体的支持。
除了GradedPool-seq法,Yang等 (Yang et al., 2019) 提出了多组混池测序成对比较分析 (Pair-wise Comparison Analysis for Multiple Pool-seq,PCAMP) 。该方法需要将具有目标性状差异明显的亲本杂交产生较大规模的F 2 代分离群体,依据表型分级分成n个group (n≥3) ,再从每个group中选择表型相同的30个个体构建DNA混池,两两之间检测混池的SNP-index和△SNP-index值确定候选区间。最后对候选区间取交集,即为最终的定位区间。该方法也可以拓展到不同的分离群体,比如DH系群体,RIL群体等。
QTG-seq法
QTG-Seq (Zhang et al., 2019a) 基本原理较复杂,简述如下:
1)群体构建:F 1 , F 2 ,F 2:3 家系,BC 1 F 1 ,BC 1 F 1:2 群体;
2)QTL定位:F 2 群体的遗传图谱QTL定位,F 2:3 家系验证QTL定位结果;
3)QTL分离:F 1 和轮回亲本回交后得到BC 1 F 1 ,利用以上QTL内分子标记筛选目标QTL杂合,其它QTL纯合的BC 1 F 1 单株 (BC1F1数量要大,保证足够的重组交换事件) ,然后自交得到BC 1 F 1:2 家系;
4)极端表型选择:从BC 1 F 1:2 中选2个极端表型组 (两极端混池各占群体总量比例20%且总共至少1000个单株) ;
5)极端表型混池:提取DNA并等量混池,形成2个极端池;
6)测序并变异检测:对两个混池测序,与参考基因组比对检测变异;
7)关联分析:结合smoothing LOD (准确关联靶标位点) ,ED (Hill et al., 2013) 与G' value (Magwene et al., 2011) 定位目标基因/QTL。
作者以玉米株高为模式性状展开QTG-seq研究:利用4代玉米群体材料,通过QTL-seq精细定位株高主效qPH7,并确定一个编码NF-YC结构域蛋白的基因Zm00001d020874为候选基因,该基因的拟南芥同源蛋白与开花时间有关。
QTG-Seq特点总结如下:
1)该方法操作较复杂,仅构建符合要求的分离群体就需要较多的时间和精力。但是相对于构建RIL群体,NIL群体又比较省时,因此该方法具备一定的优势。因为交换单株数量与定位区间大小密切相关,所以确保BC 1 F 1:2 群体有较大的混池规模至关重要;
2)将多个QTL分解为单个QTL分析 (数量性状质量化) ,在每个QTL区域选择3个 (或多个) 多态性标记筛选目标QTL杂合,其它QTL纯合的交换单株,确保交换单株后续自交发生表型分离以对目标QTL进行解析;
3)测序策略:在高的混池个体数目和测序深度之间做个权衡的话,优先保证混池的个体数量;
4)关联分析采用新方法smoothing LOD,smoothing LOD可准确关联靶标位点,并结合ED与G' value,定位精度和准确度较高,可识别候选基因。该方法综合了几种应用广泛的算法,确保了最终定位结果的可靠性。
4.2 从数学算法的角度
BSA关联方法从最开始的通过分子标记看混池条带多态性情况逐渐发展和二代测序相结合阶段,和二代测序结合也从简单的看测序后多态性标记深度的比例 (Xia et al., 2015; Zhang et al., 2015; Yu et al., 2017) 到SNP频率 (Takagi, et al., 2013a) ,G' value统计量 (Magwene et al., 2011) 以及多种方法综合 (Zhang et al., 2019a) ,再到基于多个混池的 GradedPool-Seq 中Ridit算法和滑窗降噪法 (Wang et al., 2019) 阶段。算法逐渐从简单到复杂。本部分就从算法角度探索BSA进化之旅 (以上欧氏距离法,SNP-index法本部分不再赘述) 。
Ratio_ab
Ratio_ab的公式如下:Ratio_aa = Paa/Maa; Ratio_ab = Mab/Pab。假设父本P,母本M;构建的分离群体有两个极端混池,一个是aa,一个ab。Maa指的混池aa中是来自于母本的SNP测序深度,Paa指的混池aa中是来自于父本的SNP测序深度;Mab指的混池ab中是来自于母本的SNP测序深度,Pab指的混池ab中是来自于父本的SNP测序深度。当Ratio_ab≥2,Ratio_aa≥1时 (阈值可调整) ,则认为这个SNP属于差异性的位点,这些差异位点就是潜在的跟性状关联的位点。早期的BSA就有采用这种算法 (Xia et al., 2015; Zhang et al., 2015; Yu et al., 2017) 。
Allele frequency difference
该方法是在以上的Ratio_ab的方法上改进而来,通过计算每个等位基因上SNP的reads深度,并在混池之间进行Fisher exact检验,以此判定等位基因频率差异,一般认为P-value <1e -10 ,allele frequency difference (AFD) > 0.6 的SNP位点和性状关联。由于该方法引入了Fisher exact检验,因此得到的关联SNP比Ratio_ab法更加可靠。该方法在BSR中应用较普遍 (Wang et al., 2017; Cui et al., 2019; Hua et al., 2020; Wang et al., 2020) 。
贝叶斯法
以上的Allele frequency difference在BSR中应用较为普遍,BSR中另外一种算法是贝叶斯法 (Liu et al., 2012) ,算法过程过于复杂,笔者就不详述。由于贝叶斯法需要消耗较多的计算资源,目前基于贝叶斯算法的BSA相关文献并不多。
Fst法
针对亲本高杂合构建的F 1 群体,除了欧氏距离法 (Zhang et al., 2019c) 和SNP-index法 (Huang et al., 2020a) ,Fst法也是一种备选方法 (Gu et al., 2018; Jiang et al., 2019) 。Fst又称固定系数,反映群体等位基因杂合度。其计算公式是:Fst=(H T -H S )/H T ,H T 是种群 (所有亚群) 的期望杂合度,H S 是亚群 (每个亚群) 预计杂合度。Jiang等 (Jiang et al., 2019) 在罗非鱼耐盐性QTL定位时便用到了Fst法,作者构建了一个全同胞家系,依据耐盐性差异构建了两个极端混池。基于两混池的Fst值和a=0.001阈值 (Fst=0.297) 筛选,将该性状的QTL定位在罗非鱼第18号染色体的3.3 Mb的区间,该区间也能通过全基因组关联分析 (GWAS) 重复检测到,证明了该方法的有效性。Knorst 等 (Knorst et al., 2019) 在定位黑麦草青枯病抗性QTL时也用到了该算法。目前Fst法只在家系群体和F 1 群体中有文献报道,在其他群体中的应用效果还有待观察。
Allele frequency directional difference and density法
2018年,Dougherty等 (Dougherty et al., 2018) 依据F 1 群体的遗传特点,提出了allele frequency directional difference and density (AFDDD) 法。该方法在混池、测序和SNP检测上与其他方法并无太大区别。区别就在于作者将SNP进行二等位遗传编码,编码成以下6种形式:ab×cd,ef×eg,hk×hk,lm×ll,nn×np,qq×qq。首先删除qq×qq这种类型的变异,因为该种类型在混池之间无多态性。其次,双亲杂交群体中某个性状的等位基因一般只含有2个,所以含有三个或者四个等位基因的类型ef×eg和ab×cd需要被删除。最终剩下hk×hk,lm×ll,nn×np三种类型的变异。另外,如果显性且显性杂合混池排序在前的情况下,nn×np类型的变异也需要被舍弃。这样就只剩下hk×hk和lm×ll两种编码类型的SNP被留下来用于后续的分析。然后依据混池中不同的编码类型分别计算MAFD值 (Mutant allele frequency and density) 和AFDDD (Allele frequency directional difference and density) 值。如果某个基因组区域的变异密度和平均值相比有一个显著的增长,则该区间被认为与目标性状连锁。为了确定与目标性状显著连锁的位点/区间,显著性采用标准化后的Z检验 (Standard score (z) test) 。通过Z检验即可最终确定该混池的显著连锁区间。
该方法首次将F 1 群体的变异依据遗传分离情况进行分类,选择合适的变异用于表型和性状的连锁分析,这样能大大降低高杂合型物种背景的复杂度,进而减少无关变异对结果的干扰。其次,AFDDD法与传统的△SNP-index法相比,考虑了变异的密度,似乎更有利于性状的定位。最后,该方法不需要对亲本进行测序。需要强调的是,该方法虽然不需要亲本的数据,但是由于分析时去除了一部分变异类型,因此笔者建议两个混池的测序深度更高一些,以便能有足够多的分类后的变异可以利用。本方法的遗憾之处在于,分析步骤略显复杂,作者并未公布分析流程的生信代码,对于无生信背景的科研人员而言,使用起来有较大难度,因此可能会限制该方法的使用。目前该方法在苹果的柱状生长表型 (Dougherty et al., 2020) 和苹果的高果酸表型 (Ban et al., 2020) 性状定位上取得了进展。
G' value
G' value是对G统计量进行smoothing平滑处理后的一个新的统计量,G' value的BSA法 (Magwene et al., 2011) 计算过程较为复杂,笔者在此也不详述。应用方面,自该方法提出以后,有较多的研究使用了G' value进行BSA性状定位 (Win et al., 2017; Selby et al., 2018; Dong et al., 2019; Liao et al., 2019; Shen et al., 2019) ,群体类型有F 1 ,F 2 等。有研究表明杂合型的QTL计算G' value时具有明显的无偏性 (Unbiasedness) 和低假阳性率 (low false positives)(Shen et al., 2019) ,因此该方法在F 1 类型群体里有更广阔的应用空间。
BRM
We noticed that “区块回归定位” (Block Regression Mapping,BRM) 的方法 (Huang et al., 2020b) 。Briefly,该方法首先计算每个SNP位点的频率以及数学期望值,并进行fisher exact检验。然后通过不同block区块内SNP的算术平均值或者是深度加权的平均值计算最优的block大小。最优区块大小则使用在给定最少标记的前提下,有效区块数最多时的大小。之后使用显著率 (significant rate,SR) 作为度量的统计量。通过计算区块 (block) 内显著标记占总标记数的比例来获得SR值 (significant rate) ,ANLP值 (average negative logarithm of p value) ,AAFD值 (average allele frequency difference) 和AAFU值 (average allele frequency in unselected pool) 。下一步采用LOESS拟合区块的散点,并用最优参数由赤池信息准则 (Akaike information criterion,AIC) 来决定确定最优的拟合参数。基于LOESS拟合后的峰值,使用QTL支撑置信区间 (supported confidence interval,SCI) 估计最终的QTL置信区间,即为最终的性状关联区间。
其他方法
以上几种算法是是进行BSA分析用到较普遍的算法,另外还有一些比较小众的算法,比如基于动态贝叶斯网络的MULTIPOOL法 (Edward et al., 2012) ,基于隐马尔科夫链的EXPLoRA法 (Pulido-Tamayo et al., 2016) ,index scores法 (Watanabe et al., 2017) 和PyBSASeq (zhang et al., 2019b) 等,笔者就不一一详解了。
图3. 不同群体不同算法和不同分析软件包的对应关系
4.3 不同算法结果比较
BSA算法非常多,有研究对欧氏距离法 (ED) ,|△SNP-index|,G' value三种应用广泛的方法的性能进行了评测 (Shen et al., 2019) 。其主要结果如下:ED,|△SNP-index|和G' value法均显示出了对主效QTL位点检测的一致性,说明这三种方法均具备主效QTL稳健的检测能力。|△SNP-index|法得到的区间最大,显著性低于欧氏距离法和G' value法,说明|△SNP-index|法对QTL的敏感性相对较低。欧氏距离法检测出了一些假阳性的QTL位点。G' value法检测的QTL具有明显的无偏性和较低的假阳性率,尤其是对于亲本均杂合的QTL,同时G' value的定位区间也要明显小于欧氏距离法。结论显示G' value法综合性能最优。
5. 关联定位结果阈值决定
BSA的阈值标准众多。首先,由于BSA的算法众多,每种算法都会对应不同的阈值;其次,即便是同一种算法,阈值有时候也是标准不一。下文笔者就来介绍下BSA的阈值情况。
SNP-index算法的阈值
SNP-index算法阈值由Abe等 (Abe et al., 2012) 描述SNP-index算法时一并提出。在无QTL的零假设下生成SNP-index的置信区间。具体操作如下:基于计算模拟,首先通过随机抽样建立两个具有给定数量个体的后代池。从每个池中采取给定数量的等位基因和对应的测序reads深度信息,然后计算SNP-index和△SNP-index值,如此重复10000次后可生成90%,95%和99%的置信区间。此方法每个SNP位点都会对应3个阈值点,因此该阈值线是动态变化的。应用的案例也很多 (Kumar et al., 2018; Cao et al., 2019; Hao et al., 2019; Zhang et al., 2019d; Tan et al., 2020) 。
为了计算起来更简便快捷,SNP-index法阈值确定又做了一些细微的调整。在以上10000次模拟计算后,把每次的计算结果做一个排序,分别取位于95%,99%处的值作为阈值。这个时候的阈值就是一个恒定值。该阈值确定方法也有较多的文献采用 (Sun et al., 2018; Wen et al., 2019) 。
除了以上利用计算机模拟10000次得到的阈值之外,还有更加简单的阈值确定方法。比如Feng等 (Feng et al., 2018) 定位酵母抗冻性QTL时,取△SNP-index值的99%分位数作为阈值。也有其他研究采用了99%的分位数作为阈值 (Zhou et al., 2017; Zheng et al., 2018; Xu et al., 2018) 。当然有了99%分位数,也有95%的分位数作为阈值的案例 (Zhao et al., 2016) 。
其他阈值
1)利用遗传分离比确定阈值,比如在一个F 2 群体中,△SNP-index阈值取0.67 (Li et al., 2018; Liu et al., 2019; Xu et al., 2019) 。
2)利用滑窗法 (slide window) 拟合SNP-index时,每个window内通过卡方检测确定标记和性状的连锁性 (Shu et al., 2018) 。
3)所有SNP-index的平均值+标准差作为阈值 (Xu et al., 2015) 。
ED算法的阈值
ED算法阈值也众多:有平均值+3 X标准差 (Wang et al., 2018; Yin et al., 2018; Zhao et al., 2018; Du et al., 2019; Shen et al., 2019; Li et al., 2020) 。有ED 4 >0.1 (非染色体版本参考基因组)(Zhang et al., 2019c) 。也有ED 5 的99%分位数 (Su et al., 2016) 。
Allele frequency difference (AFD) 算法的阈值
AFD算法目前出现较多的有三种不同的阈值:allele frequency difference (AFD) > 0.8 and Fisher's exact test P-value <1e-10 (Li et al., 2018) 。allele frequency difference (AFD) >0.5 and Fisher’s exact test P-value < 1e-5 (Zhao et al., 2020) 。allele frequency difference (AFD) >0.8 and Fisher’s exact test P-value <1e-10 (Wang et al., 2020) 。
BRM法的阈值
BRM法通过Fisher确切概率测验计算每个SNP标记的P值,然后设定总体显著性水平,比如0.05。即可在全基因组上获得显著的关联标记。然后以一定大小区间作为间隔 (比如0.5kb) ,选择不同block大小 (比如1kb,1.5kb,2kb一直到10kb) ,当有效block达到最大的窗口即为最佳窗口。之后计算每个block内的SR值 (significant rate) ,ANLP值,AAFD值和AAFU值,通过LOESS拟合形成平滑曲线,最后使用QTL支撑置信区间 (supported confidence interval,SCI) 估计最终的QTL置信区间目前已有文献报道水稻抗稻瘟病的性状定位上用到该方法 (Liang et al., 2020) ,证明了该阈值方法的有效性。
BSA的阈值标准不一,笔者认为阈值选择更偏实用主义pragmatism和empirical。不管何种标准,最终能定位到目标区间即可,毕竟多数情况下BSA定位就是初定位,后续对初定位区间的验证以及精细定位还需要进行其他实验设计。
6. 影响BSA定位结果的干试验因素
影响BSA定位结果的因素很多 (Guo et al., 2016) 。上文我们也讨论了不同算法对BSA定位结果的影响。但是本部分聚焦在当前已有的数据和确定的分析方法下,影响关联分析结果的两个重要因素。
6.1 参考基因组质量
我们有提到过,进行BSA分析虽然有一些针对无参考基因组物种的解决办法 (Nordström et al., 2013) ,但是一般还是需要有参考基因组。事实上,随着测序技术的不断进步,测序的质量上升而成本极大地降低,不同物种的参考基因组逐渐被释放,某些常见物种还释放了多个参考基因组。Shen等 (Shen et al., 2019) 在定位苹果果实轮纹病抗性QTL时,便使用了两个版本的基因组比较BSA的定位结果。作者分别利用2010年 (Velasco et al., 2010) 和2017年 (Daccord et al., 2017) 释放的两版苹果参考基因组,基于同一套混池数据 (抗感轮纹病混池) 和同一套分析流程,分别对抗轮纹病性状的QTL进行了定位。2010版参考基因组和2017版参考基因组分别鉴定出14个和12个显著的QTL位点。2017版参考基因组中7个可靠的QTL和2010版参考基因组恰好重合或者有重叠区间,5个可靠的QTL与2010版参考基因组定位结果完全不匹配。相反,2010版参考基因组中有2个可靠的QTL (G09.1和G09.2) 在2017版参考基因组中并未检测到。序列比对后发现,两版参考基因组在这两个QTL区域存在明显的gaps,该区域2010版参考基因组有组装好的序列,而2017版没有,这可能就是2017版参考基因组检测不到这两个QTL的原因。另外,2010版参考基因组检测到的一个QTL (J05.1) 能比对到2017版参考基因组的第10号染色体,并和QTL J.10.2重合。这些结果表明了QTL检测的可靠性,同时也表明组装完整度更高的2017版参考基因组能获得更好的结果。但是,由于2017版参考基因组使用的是DH系材料,会因为重组而丢失一些遗传片段,增加了可能丢失一些QTL位点的风险。
在另一篇GradedPool-seq报道中,Wang等 (2019) 也测试了不同参考基因组版本对性状定位的影响。首先基于同样的混池数据,作者利用粳稻日本晴2013年释放的IRGSP 1.0 (Kawahara et al., 2013) 和2009年释放的IRGSP 4.0 (http://rice.plantbiology.msu.edu/pub/data/Eukaryotic_Projects/o_sativa/annotation_dbs/pseudomolecules/version_4.0/) 两个参考基因组版本定位籽粒宽度 (GW) 性状,发现定位的结果几乎一致。作者又对千粒重性状进行了定位,同样两个版本参考基因组定位的结果几乎一致。考虑到水稻基因组组装材料是自交系,两个不同的版本参考基因组的组装水平差异并没有上文中两版苹果参考基因组大。本文的材料均为籼稻,于是作者又用了基于三代测序平台 (PacBio) 组装的籼稻参考基因组Minghui63重新定位籽粒宽度性状位点,仍然发现定位的结果别无二致。说明参考基因组对GradedPool-seq定位结果影响不大。
同样的研究为啥会得到似乎不同的结论呢?笔者认为有以下原因:首先,苹果属于高杂合物种,基于二代测序组装的参考基因组会比较碎 (contig N50=13.4kb) ,三代测序使组装的contig长度提升了接近50倍,这种参考基因组质的飞跃极大地提升了SNP calling的效率。其次,水稻参考基因组材料为自交系品种,杂合度极低。三代参考基因组和二代参考基因组相比,虽然质量上也有较大提升,但是更多的是提升着丝粒等重复序列较高区域的组装完整度。由于作者测试的是有主效基因控制的数量性状 (GW) ,潜在的QTL位点不多,而且主效位点位于3号染色体的末端区域,所以不同染色体版本对水稻的GW性状定位差异并不大。最后,被选择的水稻的几个参考基因组中均存在GW的等位基因,而如果某些稀有变异存在于不同的参考基因组上,则针对该稀有变异的表型的BSA定位结果大概率会不相同。
因此综上,更完整的参考基因组更有利于获得更好的BSA定位结果。
6.2 滑窗法(slide window)拟合窗口大小
在进行BSA关联分析时有一个细节往往会被忽略,即滑动窗口的大小。进行欧氏距离或者SNP-index等算法计算每个SNP位点的值之后,需要用滑窗法对SNP值进行拟合以降低噪音。那么滑窗法中的滑动窗口大小对BSA定位影响如何呢?Shen等 (Shen et al., 2019) 对滑窗大小对BSA定位结果的影响也进行了评估。结果发现:选择不同滑动窗口后,部分窗口下四种病原体抗性的QTL定位区间被显著缩小。其中若干主效QTL平均从1.69 Mb减少到744 Kb,G' value值也从5.80上升到6.25。而几个区间较大的QTL被分成两个或者更多较窄的子峰区间。区间大小从4.8 Mb缩小到1 Mb和500 Kb,平均G' value值也从8.93 增加到最高9.71。以上结果表明通过优化滑动窗口的大小可以缩小定位区间、提升定位的可靠性,但是遗憾的是,以上的研究并没有继续深入研究最优滑动窗口如何设置。Huang等 (Huang et al., 2020) 对最优窗口大小判定进行了探索。他们以固定长度0.5kb为间隔,并设定block区间1kb-10kb,调试后发现block在2.5kb时,有效区块数量最大,以此判定2.5kb是最合适的block大小。笔者认为由于不同物种,不同群体类型,不同混池规模和不同测序深度都会影响定位区间大小和可靠性,在滑动窗口大小选择上也不是一成不变的,后续研究的一个方向是如何让最优窗口大小的判定更加智能化。
7. 如何快速实现精细定位
精细定位最经典的策略是先进行初定位,然后在定位区间内开发分子标记,通过回交 (自交) 扩大遗传分离群体后寻找重组交换单株实现精细定位,或者构建NIL系后寻找重组交换单株实现精细定位。该途径耗时漫长,甚至有些基因克隆耗时长达10年之久 (Fan et al., 2016) 。NGS技术的大规模应用极大地提升了基因精细定位的进程。mutmap (Takagi et al., 2015) , Gradedpool-seq (Liu et al., 2020) 可以实现精细定位和候选基因获得,但是目前绝大多数BSA相关报道只是对性状的初定位,通常定位的分辨率在若干Mb。所以如何进行遗传设计或者在BSA定位的基础上继续进行精细定位是我们将要讨论的内容。笔者提出以下几个策略:
1)加大F 2 群体的构建和极端混池的规模。F 2 群体构建时间短,花费也较少,适合达到数千甚至上万个单株的规模。更大规模的子代能携带更多的重组交换,从而能提升重组活跃区的定位分辨率。Xu等 (Xu et al., 2018) 构建了一个包含949个子代的F 2 大群体定位黄瓜耐水淹性状QTL,表性鉴定后选择50+50的极端池,综合欧氏距离法 (ED) 和SNP-index法将目标位点定位在301 Kb的区间。
这类案例的特点是物种参考基因组较小 (通常1 Gb以内) ,构建的分离群体数量多 (数千株) ,挑选的极端个体也多 (50+50以上) ,此时往往能达到意想不到的效果。
2)第二种思路是进行BSA定位获得一个峰值较高 (可信度高) 的定位区间,然后在这个区间设计均匀分布的SNP-based KASP标记 (或者CAPs标记等) ,再在未入极端混池组的分离群体中随机选择几百个个体进行基因分型,然后构建这个定位区域的遗传图谱和QTL定位 (Huang et al., 2017; Liu et al., 2019) 。如果BSA定位的区间已经比较小 (一般1 Mb以内) ,则可以在群体中选择未BSA混池的极端表型个体,并基因分型后寻找重组交换单株来缩小定位区间。
该思路同样需要构建一个较大规模的分离群体或者重新构建新的大规模的遗传分离群体,利用KASP快速基因分型,理想状态下,能在较短的时间内 (1-2年) 实现一年多熟 (multiple cropping in one year) 类作物的精细定位。
3)联合其他组学。比如,BSA和BSR联合进行共定位和定位区间内基因差异表达分析 (Borovsky et al., 2020; Zhao et al., 2020) ,BSA和转录组联合筛选定位区间内差异表达基因 (Kodama et al., 2018; Hao et al., 2019; Wen et al., 2019; Ou et al., 2020) 。BSA和GWAS一起联合分析,通过自然群体的加持提升定位的分辨率 (Du et al., 2018; Jiang et al., 2019; Su et al., 2019) 。BSA和多组学联合分析也是实现精细定位和候选基因筛选的有力途径。
参考文献请查看论文原文。
论文链接:
https://onlinelibrary.wiley.com/doi/abs/10.1111/tpj.15646
特别声明:以上内容(如有图片或视频亦包括在内)为自媒体平台“网易号”用户上传并发布,本平台仅提供信息存储服务。
Notice: The content above (including the pictures and videos if any) is uploaded and posted by a user of NetEase Hao, which is a social media platform and only provides information storage services.