遗传统计|BSA-seq混合分组鉴定功能基因
分离体分组混合分析法（Bulked Sergeant Analysis, BSA）也称为集群分离分析法或混合分组分析法，是从近等基因系分析法演变而来的。常规的育种通常从近等基因系（NIL）的构建开始。近等基因系指一组遗传背景相同或相近，只在个别染色体区段上存在显著差异的株系。但是并非所有植株都能够简单构建近等基因系。因此衍生出了集群分离分析法或混合分组分析法，即BSA。Michelmore等人在1991年首先提出分离体分组混合分析法[1]。
Michelmore老爷子1991年画的，简朴中透露着一丝牛逼。91年测序技术尚未普及。那个年代书信很远，车马很慢，漫漫育种路只够爱RAPD技术。
随后，随着群体分组技术的改进和二代测序(NGS)成本的显著降低，全基因组重测序已经与BSA相结合。BSA与NGS的结合(BSA-seq)加快了重要性状紧密连锁标记的鉴定，提高了基因鉴定和QTL定位图谱的分辨率。更重要的是BSA-seq在定位数量性状位点（Quantitative trait locus, QTL）时能够节省大量构建群体的时间。因此成为了快速高效的功能基因或者数量性状位点鉴定手段。总而言之，BSA-seq基于群体分组混合以及二代测序技术，能够鉴定出跟表型关联的单核苷酸多态性位点（SNP）。目前已经被广泛用于定位QTL和功能靶基因[2]。
BSA-seq技术的原理
BSA法克服了许多作物没有或难以创建NIL的限制。其原理是从具有相对性状的一对亲本杂交，在其任一分离后代群体中（一般选取F2代交配群体），选择具有目标性状极端表型的20～50个单株，将其DNA提取并且等量混合，构建两个基因池（pool）。这两个基因池应在感兴趣的性状方面存在差异，此外的所有的位点均随机化。对两个池筛选出的多态性标记可能与感兴趣的某个功能基因或QTL连锁。有关基因组重测序和BSA结合的方法最权威的是来自于Ryohei Terauchi大佬在2011年发表的NBT，介绍了mutmap软件来鉴定功能性状位点。Mutmap亲测简单易用，需要配置好python3环境。这篇NBT涉及的具体BSA方法和实验流程图也放在了下方以供参考。
Nature Biotechnology, 2011., Akira Abe et al.
BSA-seq数据分析
为了鉴定出功能基因位点，需要对混池并测序的两个群体进行SNP层面的比较。最常用的方法是SNP-index法。其原理是利用测序reads对每个碱基位点的碱基进行统计，选择某一亲本或已有的参考基因组为参考，统计子代池中和亲本或者已有参考基因组在某一个碱基位点相同或者不相同的reads条数，计算不相同reads条数占总条数的比例，即SNP-index=不同的reads条数/该位点reads总数。在计算时，通常用1Mb作为窗口滑动计数SNP-index，每次滑动增量为10kb。在此基础上，Δ(SNP-index)用于衡量两个基因池SNP指数的差值，即两池各个位点的差值。最终绘制成曼哈顿图以显示与目标性状显著关联的基因区域。如下图，横轴代表染色体位置，SNP-index图的纵轴代表每个SNP位点计算出的reads比例。Δ(SNP-index)图的纵轴代表两个混池的SNP-index的差值。Δ(SNP-index)的绝对值越大代表与极端性状关联的位点越显著。
SNP-index和Δ(SNP-index)计算结果示例
BSA-seq技术的应用实例
（1）基于简化基因组的BSA-seq技术用于鉴定辣椒果色基因
辣椒中c1、c2和y位点是控制果实颜色的关键调控基因。JeMin Lee等人将具有橙色果实的辣椒栽培种Habanero Orange与具有桃红色果实的Habanero Peach杂交得到381株F2群体[3]。并在F2群体中选取20-25个单株进行DNA混池及简化基因组测序（GBS）。下机得到的数据以C.chinense基因组作为参考，构建了亲本之一Habanero Orange的SNP图谱。再以后者为参考将F2子代中的橙色果实和桃色果实的数据进行比对，从而分别得到橙色果实和桃色果实的SNP图谱。此时，将得到的两个F2子代的SNP图谱与亲本Habanero Orange的SNP图谱进行比较计算。得到F2两子代各自的SNP-index，那么此时与桃色相关的基因所关联的SNP的SNP-index值应该在子代桃色果实和亲代桃色果实中接近1，在亲代橘色果实中小于0.5。该研究以此为标准计算Δ(SNP-index)，定位出11个可能的候选SNP位点，并将SNP位点扩大到CAPS位点并进行跑胶验证。辅以后继生化实验验证。
Figure1. BSA-seq用于鉴定辣椒果色基因过程中，子代极端性状混池测序群体的构建。
（2）BSA-seq与RNA-seq联合用于鉴定水稻耐冷基因
水稻耐冷性状由一系列数量性状基因位点调控。Guo Z等人利用耐冷水稻Longjing25 和冷敏感水稻Longjing11作为亲本杂交得到F2子代[4]。在F2子代群中各选择50株耐冷和冷敏感水稻构建了两个基因混池进行测序。通过计算混池的SNP-index和混池间的Δ(SNP-index)定位到了6个候选基因，分布在6号和9号染色体，跨度3.67Mb。Δ(indel-index)也作为另一个检测指标，与的Δ(SNP-index)共同定位候选基因。将候选基因最终定位在0.82 Mb区域包含98个注释基因。进一步作者对水稻进行了冷胁迫处理并且开展RNA-seq以精细定位耐冷基因。RNA-seq显示在这98个注释基因中，87个基因在冷胁迫下表达，而这其中50个基因差异表达。GO和KEGG分析定位到9个通路包含10个差异表达基因，QRT-PCR进一步验证其功能。
同样使用RNA-seq辅助BSA-seq鉴定功能基因的案例还有鉴定花菜早熟基因。Honghui Gu等人利用早熟和非早熟花菜杂交后的F2代进行混池测序，计算两池间的Δ(SNP-index)，将候选基因定位到了染色体4Mb的区域。RNA-seq鉴定出二者的差异表达基因，在Δ(SNP-index)显著差异区搜索差异表达基因。最终定位到了MADS-box蛋白SOC1[5]。
（3）BSA-seq与GWAS联合鉴定西瓜雄全同株调控基因
植物激素乙烯往往决定植物花朵的性别。但是植物雄全同株性状不止取决于乙烯，还包括了一个名为pa的隐性基因。Manuel Jamilena等人利用BSA-seq和全基因组关联分析（GWAS）联合鉴定了西瓜pa基因[6]。首先，作者构建了雌雄同株的混池与部分雄全同株的F2混池，并进行全基因组BSA-seq测序，比较两池间的Δ(SNP-index)如下图C。PA表示部分雄全同株，M表示雌雄同株。定位候选基因所在的区域大约0.5Mb。在此之上，作者利用Cucurbit Genomics Database数据库获取了122株西瓜SNP图谱进行GWAS分析。发现GWAS数据定位的区域与BSA-seq定位区域近似，进一步证实了功能位点的可靠性。在随后的基因精细定位过程中，作者利用了基因注释，外显子区域筛选等方式确定候选基因。
Figure2| BSA-seq用于鉴定西瓜性选择过程中，子代极端性状混池测序群体的构建及混池测序后两群体间SNP-index和Δ(SNP-index)情况，候选基因被定位于C中红色方框区域内。
目前，分离体分组混合分析法（BSA）结合测序技术目前已经被用于多种农业作物的基因定位，包括水稻，棉花，黄瓜，大豆和玉米[2]。其中RNA-seq与之结合的案例相对较多。其研究思路聚焦于RNA-seq检测到的差异表达基因和BSA-seq鉴定到的差异SNP指数区域，二者交集的位置可能包含潜在的功能基因。RNA-seq是缩小候选基因范围的有效工具。此外，在上述鉴定水稻耐冷基因的案例中[4]，indel-index与SNP-index同时作为筛选候选基因位点的指标，从而缩小筛选功能基因的筛选范围。
最后做一下总结
BSA-seq技术鉴定候选功能基因依赖于其混池群体的构建。目前混池群体的构建主要使用的是拥有一对极端性状的亲本杂交后的F2代，通常选取20-50株进行混池测序。相较于全基因组关联分析来说群体量更小。经过BSA-seq所定位到的候选基因区间在0.5Mb-4Mb之间，其中所包含的基因多者能达到数百个。精细度上相较于全基因组关联分析来说要更低。提高BSA-seq鉴定基因精确度，缩小功能基因所在的范围可以采用多个分子标记联合筛选的办法。开发具有多态性率高的不同种类标记，例如SSR、InDel和SNP，同时用于靶基因的精细定位[4]。
多组学的联合也为BSA-seq精细定位候选功能基因提供了强有力的支持。RNA-seq与BSA-seq的联合能够在转录组层面提供注释信息，从而更加精确地缩小候选功能基因的范围，以便于后继基因的精细定位。在基因组层面鉴定功能基因的主要方法还有全基因组关联分析（GWAS），二者的结合使用也可以一定程度上缩小功能基因的筛选范围，为功能基因筛选提供更多依据。
参考信息：
1. Michelmore RW, Paran I, Kesseli RV. Identification of markers linked to disease-resistance genes by bulked segregant analysis: a rapid method to detect markers in specific genomic regions by using segregating populations. Proc Natl Acad Sci U S A. 1991;88(21):9828-9832. doi:10.1073/pnas.88.21.9828
2. Zhang K, Li Y, Zhu W, et al. Fine Mapping and Transcriptome Analysis of Virescent Leaf Gene v-2 in Cucumber (Cucumis sativus L.). Front Plant Sci. 2020;11:570817. Published 2020 Sep 25. doi:10.3389/fpls.2020.570817
3. Lee SB, Kim JE, Kim HT, Lee GM, Kim BS, Lee JM. Genetic mapping of the c1 locus by GBS-based BSA-seq revealed Pseudo-Response Regulator 2 as a candidate gene controlling pepper fruit color. Theor Appl Genet. 2020;133(6):1897-1910. doi:10.1007/s00122-020-03565-5
4. Guo Z, Cai L, Chen Z, et al. Identification of candidate genes controlling chilling tolerance of rice in the cold region at the booting stage by BSA-Seq and RNA-Seq. R Soc Open Sci. 2020;7(11):201081. Published 2020 Nov 18. doi:10.1098/rsos.201081
5. Zhao Z, Sheng X, Yu H, Wang J, Shen Y, Gu H. Identification of Candidate Genes Involved in Curd Riceyness in Cauliflower. Int J Mol Sci. 2020;21(6):1999. Published 2020 Mar 15. doi:10.3390/ijms21061999
6. Aguado E, García A, Iglesias-Moya J, et al. Mapping a Partial Andromonoecy Locus in Citrullus lanatus Using BSA-Seq and GWAS Approaches. Front Plant Sci. 2020;11:1243. Published 2020 Aug 19. doi:10.3389/fpls.2020.01243
