大豆非同义SNP相邻碱基组分分析

2022-07-12王娟，常玮

南阳师范学院学报 2022年4期

王娟，常玮

(南阳师范学院生命科学与农业工程学院，河南南阳 473061)

单核苷酸多态性(single nucleotide polymorphism，SNP)是指在基因组上单个核苷酸的变异，包括转换、颠换、缺失和插入.据报道，在人类基因组中，每1000～2000个碱基就有一个SNP，在一些高频的区域，甚至能达到每300个碱基一个SNP[1-2].由于SNP在基因组中的高密度，使其成为一种理想的分子标记，并已被用于构建人类遗传分析图谱和人类复杂疾病性状的研究[3-5].SNP在植物中也广泛存在，例如在水稻基因组中的发生率接近1/268[6].根据SNP在基因中的位置，可将其分为基因编码区SNP(coding SNP，cSNP)、基因内含子区SNP和基因调控区SNPs(regulatory SNPs，rSNP).其中cSNP根据是否改变编码的氨基酸类型又可分为同义cSNP(synonymous SNP，sSNP)和非同义cSNP(non-synonymous SNP，nsSNP).其中nsSNP由于能引起蛋白质编码序列的改变，因而以其作为研究基因功能的一种功能标记.这使得nsSNP成为遗传研究中应用前景广阔的理想分子标记.

目前，开发SNP标记(SNPs)的主要方法有：EST序列挖掘法[7]、基因芯片检测法[8]、基因组比较法[9]、扩增子重测序技术[10]以及基于下一代测序技术(Next-Generation Sequencing，NGS)的SNPs挖掘法[11].随着测序技术的提高与成本的下降，NGS技术已经成为SNPs挖掘的主要方法.例如，基于NGS技术的RNA-seq研究，由于可以直接获得基因转录本的序列，也为挖掘nsSNP标记(nsSNPs)提供了基础[12].尽管基于NGS的RNA-seq方法可以获得存在于细胞中所有的mRNA分子，但由于RNA剪接的复杂性会导致假阳性结果的出现[13]，因此如何获得更为可靠的nsSNPs，已成为功能基因分子标记开发、分子辅助育种的关键.随着测序数据的不断增长和计算机技术的发展，目前已经出现了许多降低SNP假阳性率的算法与程序，如基于决策树的机器学习法[14]、ISMU(Integrated SNP Mining and Utilization)和GBS-SNP-CROP等程序包[15-16]以及具有交互式图形用户界面的软件PLANET-SNP[17].这些算法及程序充分利用了SNP位点自身的信息进行判断，并能够显著降低SNP假阳性率，但在计算过程中没有分析SNP位点侧翼序列信息.赵辉等研究发现，基因组的碱基组分在不同程度上与碱基替换突变的发生相关[18].该研究同时表明水稻A&T碱基影响突变的范围局限在突变位点上下游2个碱基内，而拟南芥A&T碱基的影响范围不超过上下游4个碱基.因此对SNP位点侧翼序列成分进行分析统计，可以进一步降低SNP假阳性率[18].

大豆是重要的粮食作物，相比于玉米等作物，近年来产量提升缓慢，分子标记辅助技术将是加快育种进程的重要手段，而开发大豆nsSNPs是实现大豆分子育种的基础.本研究拟采用生物信息的方法对大豆ESTs数据中所包含的候选SNPs进行挖掘，结合数据库中公布的大豆SNPs，依据其编码蛋白的情况挖掘其中的nsSNPs，并对这些nsSNPs分布规律及其附近的序列特征进行统计分析，以期为降低nsSNPs假阳性率提供理论依据.

1 材料与方法

1.1 材料

1.1.1 序列数据

从GenBank/dbEST(http：//www.ncbi.nlm.nih.gov/dbEST/)下载共1 453 823条大豆EST序列,大豆基因组序列及大豆cDNA序列来自phytozome数据库(http://www.phytozome.net).

1.1.2 软件工具

用于序列拼接的程序：CAP3[19]；用于本地化的序列比对程序：BLAST[20]；用于在Windows平台上执行Perl程序的工具软件：ActivePerl；在Windows平台上运行的unix模拟环境：Cygwin；用于SNPs是评估的程序：MAVIANT[21]；用于非同义SNPs挖掘的工具：UltraCompare；用于序列分析的Perl程序工具包：BioPerl[22].

1.2 方法

1.2.1 SNPs挖掘及筛选

将全部1 453 823条EST序列用BLAST进行聚类，E值小于10-10；将每一类的序列用CAP3进行拼接，每一个重叠群必须有大于95%的相似度；将得到的拼接结果与相应的质量文件输入MAVIANT，去除低质量的结果，利用MAVIANT的SNP评估功能进行候选SNPs的筛选.根据报道的主要作物SNP的出现频率设置窗口大小为150bp，进一步对所得的SNPs进行筛选.

1.2.2 nsSNPs挖掘

将得到的包含候选SNPs的一致性序列与大豆cDNA序列进行blastn分析，从而将SNPs定位在大豆cDNA序列上，记录每个位点的位置；将包含cSNPs的编码序列，按照野生型与突变型用BioPerl中的make_mrna_protein.pl脚本进行翻译，并将生成的两蛋白序列文件用UltraCompare比较，找出其中的nsSNPs.

1.2.3 nsSNPs相邻碱基组分分析

根据赵辉等对水稻和拟南芥基因组相邻碱基组分与产生SNP的转换或颠换的研究方法[18]，对大豆SNPs相邻碱基组分进行统计分析：首先统计基因组中对应两种碱基替换类型(转换和颠换)SNP位点上下游各12个位置共24个位点组成的区域(称为背景区域)碱基A&T的使用频率，记为f0=(f01,f02)，其中f01对应于碱基A的使用频率，f02对应于碱基T的使用频率.然后统计SNP位点上下游各12个位置上A&T碱基的使用频率，记为fi=(fi1,fi2)，i=-12,-11,…,-2,-1,1,2,…,11,12,其中fi1表示第i位置上碱基A的使用频率，fi2表示第i位置上碱基T的使用频率，SNP位点的位置记为0位，将5’端为负号位，3’端记为正号位.最后计算fd=fi-f0，并根据fd的值的大小，估计该位置A&T碱基对SNP发生概率的影响程度和影响范围.

根据SNP位点上下游紧邻的两个位点上的A&T个数将其分为三类，即没有A或T，只有一个A或T(另一个是G或C)和两个都是A或T，分别用AT0，AT1和AT2记之.首先，统计这三类内6种类型SNP的个数，然后再统计每一类(AT0，AT1和AT2)转换(A/G和C/T型)和颠换(A/C，G/T，A/T和C/G型)的个数，从而得出每类的转换和颠换的比值(Ts/Tv).

2 结果与分析

2.1 大豆EST序列中SNPs的挖掘

全部1 453 823条EST序列共拼接得到50 867个contigs和183 847个single ESTs，利用CAP3产生的50 867个contigs以及与之相对应的质量文件，经筛选共得到16 093个候选SNPs.如表1所示：其中：[A/T]、[G/T]、[A/C]、[G/C]、[A/G]和[C/T]突变类型分别检测到1 903、1 516、1 396、854、4 875和5 576次.其中发生转换的次数为：[A/G]+[C/T]= 4 875+5 576=10 451次；发生颠换的次数为：[A/T]+[G/T]+[A/C]+[G/C]=1 903+1 516+1 369+854=5 642次.从结果中可以看出转换发生的次数明显高于颠换，在转换中[C/T]突变类型高于[A/G]；而在颠换中，[A/T]突变类型最多，[G/C] 突变类型最少.

进一步分析每种突变类型中两种突变事件,如 [A/T]突变类型中的两种突变事件：[A→T]与[T→A]的发生频率，结果见表1.从表1中我们可以看出，在[A/G]转换中，[G→A]发生的次数明显高于[A→G]发生的次数；在[C/T]转换中，[C→T]发生的次数明显高于[T→C]发生的次数.而其他四种颠换类型的两种突变事件的发生次数之间没有明显的差别.

表1 大豆SNP类型统计表

2.2 编码序列中nsSNPs的挖掘

将包含候选SNPs位点的16 093条一致性序列与大豆cDNA序列进行blastn分析，共有14 297条能够匹配相应的cDNA序列.将这14 297条与植物蛋白数据库比对，共获得13 888个匹配结果，保留Score≥150、Identities≥65%的结果共7 085个.搜寻这些cDNA序列上的ORF的起始位置，根据SNPs位点的位置，共获得5’UTR-SNPs：63个、3’UTR-SNPs：1 508个以及cSNPs：5 514个.根据UltraCompare的比较结果，其中sSNPs：3 748个；nsSNPs：1 766个(图1).位于编码区的SNPs占总数的77%，在这些SNPs中，sSNPs所占的比例约为nsSNPs的2倍，而在剩下的22%位于非编码区内的SNPs，约95%的位点位于3’UTR.分别统计分布于不同区域的SNPs中转换Ts与颠换Tv的比值，5’UTR-SNPs、3’UTR-SNPs、sSNPs与nsSNPs的Ts/Tv分别为2、1.5、2.25、1.7，四者的关系为：sSNPs>5’UTR-SNPs>nsSNPs>3’UTR-SNPs.

图1 大豆编码序列SNP类型统计表

2.3 相邻碱基组分与产生nsSNPs的转换或颠换的分析

分析结果表明显示，无论总SNPs还是不同分布的SNPs，均表现为Ts/Tv大于1，进一步的分析结果还表明：大豆nsSNPs邻近碱基中，A&T碱基在远离突变位点的其他位置的碱基组成无差异，而在与突变位点紧密相邻的两个碱基则表现出了一定的规律，即无论是转换(图2A)还是颠换(图2B)，突变点5’端的A&T 碱基使用频率偏差接近0；突变点3’端的A碱基使用频率偏差约为-0.05；T碱基使用频率偏差约为0.02，表明A&T 碱基的影响范围局限在下游第一个碱基.大豆nsSNPs的AT0，AT1和AT2三类中 Ts/Tv 的比值依次为1.682、1.685、1.77.结果表明：只有在两侧碱基都为A或者T时，Ts/Tv 的比值才会升高.由于突变位点上游第一个碱基无明显的频率偏差，故推测下游第一个碱基的频率偏差会导致Ts/Tv 的比值的升高，同时Ts中T的出现频率(-0.052)小于Tv中T的出现频率(-0.048)，说明在突变位点下游第一个碱基T频率偏差与Ts/Tv的比值成反比.

3 讨论

位点的可靠性是SNP挖掘的难点，一般来说可以从以下四方面进行评价：第一，具有同一突变类型的片段越多，结果就越可靠，这也是本研究采用EST数据进行SNP挖掘的主要原因；第二，由于测序的原因，相较于测序片段的两端，读序中段假阳性率更低；第三，由于连锁不平衡现象的存在，相邻两个SNP位点应紧密连锁；第四，由于编码区受选择压力更大，SNP的比例低于全基因组的水平[23].

图2 nsSNPs相邻碱基组分分布图

本研究获得的SNPs位点中，转换变异10 451个，约为颠换变异的2倍.这是由于胞嘧啶易经甲基化和脱氨作用转换为胸腺嘧啶，致使置换发生率高于颠换[23].本研究发现[G→A]发生的次数明显高于[A→G]发生的次数；[C→T]发生的次数明显高于[T→C]发生的次数，与上述结论一致.此外，由于5’UTR中包含一些调控元件，因此相比较于3’UTR受到的选择压力更大一些，所以位于3’UTR的SNPs数目远远大于5’UTR；而由于编码区受选择压影响，编码区内不改变氨基酸序列的候选SNP比改变氨基酸序列的候选SNP具有更高的存在概率，因此位于编码区的SNPs中nsSNPs数目少于sSNPs的数目[24]，本研究也得到了相同的结果.

DNA复制酶的严谨性、各种化学物理诱变因素等都会影响基因组的碱基替换变异.我们的研究发现，大豆nsSNPs位点影响范围局限在下游第一个碱基，且T的含量越低，越有利于转换的发生.这一结果与赵辉等的研究一致[18].Radman和Wagner的研究对此现象的解释为：当某个位点周围只有A或T碱基时，形成CG二联码的概率为0，因而不易发生转换突变[25].