基于基因拷贝数变异位点对林木群体基因型分型的方法与流程

本发明涉及分子生物学技术领域，尤其涉及一种基于基因拷贝数变异位点的基因型分型方法，更具体涉及一种基于基因拷贝数变异位点对林木自然群体的基因型的分型方法。

背景技术：

拷贝数变异(Copy Number Variations,CNVs)是指不同个体基因组中长度范围从1kb至数Mb的DNA区段，与参考基因组比较而发生拷贝数变异的亚微观染色体结构变异。在人类中，CNV已被证实与许多复杂疾病相关。而在植物基因组中，由于CNV的存在影响了植株的表型性状、代谢过程、生理过程和适应性进化等。因此，对群体中基因拷贝数变异进行全面研究显得尤为重要。尽管在过去的几年里，高通量测序平台的问世，使得大量检测CNV的方法和技术也陆续出现和不断完善，但对于精确的发现和检测CNV，尤其是对CNV位点基因型检测的准确计算仍具有强烈的挑战。

目前，利用高通量测序结果对CNV位点基因型进行分型大多是基于CNV位点附近的SNP位点的基因型进行的分型，该方法操作复杂、结果分析困难，不能精确的检测CNV位点的基因型，使其应用受到很大的局限。此外，在植物领域，尤其是在林木自然群体中对CNV位点基因型进行检测的方法研究更是空白。因此，研究植物的CNV分型方法对培育高产、优质、抗病性强等优良品种具有深远的意义。

技术实现要素：

鉴于此，本发明的目的是提供一种适用于不同群体大小的林木(或植物)CNV位点的基因型分型方法，操作简便易行，并且较精确地发现和检测CNV位点的基因型。

为了实现上述发明目的，本发明提供以下技术方案：

本发明提供了基于基因拷贝数变异位点对林木群体基因型的分型方法，包括以下步骤：

1)将林木物种中已经公布的基因组序列作为参考序列，操作界面的每个窗口上显示连续的、非重叠的、大小相等区域，每个窗口能空间上左右调整；

2)将所述林木物种的不同个体测序，将所述测序得到的reads比对到所述步骤1)的参考序列上，将每个窗口内比对到所述窗口上的reads数目作为读取深度信号；

3)根据每个窗口内reads的GC含量和偏差，空间上左右调整得到所述窗口的读取深度信号值，根据调整后的读取深度信号值得到不同个体在每个窗口内的读取深度信号值形成的数据集，再根据所述数据集计算得到中值；

4)根据步骤3)得到的中值对所述步骤2)得到的读取深度信号值进行数据标准化处理，得到校正后的读取深度信号值；

5)根据步骤4)得到的校正后的读取深度信号值计算不同个体在所述窗口中拷贝数，以拷贝数为2的拷贝数设为正常拷贝，正常拷贝的基因型为(1；1)；

6)当所述步骤5)中得到的拷贝数发生变异，且所述窗口内读取深度信号与无拷贝数变异的窗口相比有显著地增高时，属于基因拷贝数发生重复的结构变异；对重复拷贝数CN为3或4的基因位点进行分型；

7)所述步骤6)的分型方法具体为：将标准化校正后的读取深度信号值在1.25<读取深度信号值<1.75范围内的基因拷贝数变异位点设为杂合重复，杂合重复的基因型为(1；2)；将标准化校正后的读取深度信号值在1.75<读取深度信号值<2.25范围内的基因拷贝数变异位点设为纯合重复，纯合重复的基因型为(2；2)；

8)当所述步骤5)中得到的拷贝数发生变异，且所述窗口内读取深度信号与无拷贝数变异的窗口相比有显著地降低时，属于基因拷贝数发生缺失的结构变异；

9)将所述步骤8)中缺失的结构变异进行分型，所述分型的方法具体为：将标准化校正后的读取深度信号值<0.10的基因拷贝数变异位点设为纯合缺失，设定纯合缺失基因型为(0；0)；将标准化校正后的读取深度信号值为0.10<读取深度信号值<0.75的位点设为杂合缺失，设定杂合缺失的基因型为(0；1)。

优选的，所述步骤1)中窗口的大小固定，所述窗口的大小为500bp。

优选的，所述步骤2)中读取深度信号值由CNVnator软件计算得到。

优选的，所述步骤3)调整所述窗口的读取深度信号值的具体方法是：调整读取深度信号值使GC含量达到48～52％和偏差小于5％。

优选的，所述步骤3)中值的计算方法为：将每个个体得到的读取深度信号值的数据集按从小到大的顺序排列，位于该数列中间位置的数值为该个体的读取深度信号值的中值。

优选的，所述步骤4)中数据标准化处理按照式I计算；所述式I为x′＝x_i-Me/Std(x)，其中x′为得到的新数据，x_i为原始数据，Me为该列数据集中的中值，Std(x)为标准差。

优选的，所述步骤5)中计算不同个体在窗口中拷贝数的方法为：以2个拷贝为中心，将标准化后的读取深度信号值经过四舍五入后，取最接近的整数作为该个体在该窗口中的拷贝数。

优选的，所述步骤5)，步骤7)和步骤9)中整个步骤是基于自然群体中的平衡选择定律。

本发明提供的基于基因拷贝数变异位点对林木群体基因型的分型方法，是一种适用于不同群体大小的林木CNV位点的基因型分型方法，利用高通量测序结果的读取深度信号确定CNV位点的基因型，其算法复杂度较低，操作简便易行，可以较精确地发现和检测CNV位点的基因型。

另外，本发明提供的分型方法，利用读取深度信号值数据集的中值对读深信号值进行标准化处理，可以减少测序过程中产生的测序深度、技术和人员操作等方面的误差，降低假阳性，使整体数值更为均一化，分型结果更为精确。

说明书附图

图1为实施例1中利用读深信号值的方法在4个样本中检测到的缺失示意图；

图2为实施例1中利用读深信号值的方法在4个样本中检测到的重复示意图。

具体实施方式

本发明提供了基于基因拷贝数变异位点对林木群体基因型的分型方法，包括以下步骤：

1)将林木物种中已经公布的基因组序列作为参考序列，操作界面的每个窗口上显示连续的、非重叠的、大小相等区域，每个窗口能空间上左右调整；

2)将所述林木物种的不同个体测序，将所述测序得到的reads比对到所述步骤1)的参考序列上，将每个窗口内比对到所述窗口上的reads数目作为读取深度信号；

3)根据每个窗口内reads的GC含量和偏差，空间上左右调整得到所述窗口的读取深度信号值，根据调整后的读取深度信号值得到不同个体在每个窗口内的读取深度信号值形成的数据集，再根据所述数据集计算得到中值；

4)根据步骤3)得到的中值对所述步骤2)得到的读取深度信号值进行数据标准化处理，得到校正后的读取深度信号值；

5)根据步骤4)得到的校正后的读取深度信号值计算不同个体在所述窗口中拷贝数，以拷贝数为2的拷贝数设为正常拷贝，即无拷贝数变异；所述无拷贝数变异的基因型为(1；1)；

6)当所述步骤5)中得到的拷贝数发生变异，且所述窗口内读取深度信号与无拷贝数变异的窗口相比有显著地增高时，属于基因拷贝数发生重复的结构变异；对重复拷贝数CN为3或4的基因位点进行分型；

7)所述步骤6)的分型方法具体为：将标准化校正后的1.25<读取深度信号值<1.75范围内的基因拷贝数变异位点设为杂合重复，其基因型为(1；2)；将标准化校正后的1.75<读取深度信号值<2.25的基因拷贝数变异位点设为纯合重复，其基因型为(2；2)；

8)当所述步骤5)中得到的拷贝数发生变异，且所述窗口内读取深度信号与无拷贝数变异的窗口相比有显著地降低时，属于基因拷贝数发生缺失的结构变异；

9)将所述步骤8)中缺失的结构变异进行分型，所述分型的方法具体为：将标准化校正后的读取深度信号值<0.10的基因拷贝数变异位点设为纯合缺失，设定纯合缺失的基因型为(0；0)；将标准化校正后的读取深度信号值为0.10<读取深度信号值<0.75的位点设为杂合缺失，设定杂合缺失的基因型为(0；1)。

本发明将林木物种中已经公布的基因组序列作为参考序列，操作界面的每个窗口上显示连续的、非重叠的、大小相等区域，每个窗口能空间上左右调整。

本发明中，林木物种中已经公布的基因组序列优选为从NCBI网站中下载得到。

本发明中，提供所述操作界面的软件优选为Linux系统的Shell窗口。所述Linux系统的Shell窗口是基于大型计算机服务器。

本发明中，所述操作界面的每个窗口上显示的连续的、非重叠的、大小相等区域优选利用CNVnator算法划分。所述CNVnator算法优选利用均值漂移技术(mean-shift technique)将读取深度信号分成有潜在CNV的小片段中。所述窗口的大小固定，所述窗口的大小优选为500bp。

得到参考序列后，本发明将所述林木物种的不同个体测序，将所述测序得到的reads比对到所述参考序列上，将每个窗口内比对到所述窗口上的reads数目作为读取深度信号。所述林木物种优选为杨树，更优选为毛白杨。

本发明中，利用了林木基因杂合度高、DNA序列多态性丰富的特点，结合Illumina的pair-end和454的mate-pair两种测序方法对物种的不同个体进行测序，这种测序方案很好地结合了短序列插入片段和长序列插入片段的各自突出优势，更好地提高了测序结果的准确性，为后续的精确分型提供了保障。

本发明中，所述林木物种的不同个体测序具体是测定林木物种中不同个体的基因组。所述林木物种个体的个数优选为400～500个，更优选为435个。

本发明中，所述比对优选包括以下步骤：

Ⅰ.利用软件BWA-0.7.8中的aln算法将不同个体得到的大量reads分别比对到参考基因组序列上；

Ⅱ.以picard软件包中的Markduplicate工具标记可能的PCR重复；

Ⅲ.利用软件Samtools文件包中的flagstat工具统计得到所有样品的比对信息文件，并以bam格式进行保存。

比对结束后，本发明将同一个个体中比对到参考序列上gap位置的reads删除。

本发明中，所述窗口的大小优选为固定值，所述窗口的大小优选为500bp。

本发明中，所述读取深度信号值优选由CNVnator软件划分的连续的、非重叠的、大小相等区域中映射的reads数目计算得到。利用读取深度信号值分型的思路为：假定读取深度是一个泊松分布，然后利用读深信号的随机分布检测目标样本中的重复和缺失类型的CNV，并进行分型。

得到每个窗口的读取深度信号值后，本发明根据每个窗口内的GC含量和偏差，空间上左右调整得到所述窗口的读取深度信号值，根据调整后的读取深度信号值得到不同个体在每个窗口内的读取深度信号值形成的数据集，再根据所述数据集计算得到中值。

本发明中，所述调整所述窗口的读取深度信号值的具体方法优选是：调整读取深度值使GC含量达到48～52％或偏差小于5％。

本发明中，所述步骤3)中值的计算方法优选为：将每个个体得到的读取深度信号值的数据集按从小到大的顺序排列，位于该数列中间位置的数值为该个体的读取深度信号值的中值，也称中位数。

得到中值后，本发明根据中值对所述读取深度信号值进行数据标准化处理，得到校正后的读取深度信号值。

本发明中，所述数据标准化处理按照式I计算；

x′＝x_i-Me/Std(x) 式I

其中为x′得到的新数据；

x_i为原始数据；

Me为该列数据集中的中值；

Std(x)为标准差。

得到校正后的读取深度信号值x’后，本发明根据校正后的读取深度信号值计算不同个体在所述窗口中的拷贝数，以CN为2的拷贝数设为正常拷贝，即无拷贝数变异，所述无拷贝数变异的基因型为(1；1)。

本发明中，所述计算不同个体在窗口中拷贝数的方法优选为：以2个拷贝(即二倍体)为中心，将标准化后的读取深度信号值经过四舍五入后，取最接近的整数作为该个体在该窗口中的拷贝数。

本发明中，当个体某个基因存在结构变异时，个体样本所测得到的reads映射到窗口中的读取深度信号较无拷贝数变异的区域会有显著地增高或降低，说明在该区域发生了基因拷贝数的重复(duplication)或缺失(deletion)。

当所述窗口内读取深度信号与无拷贝数变异的窗口相比有显著地增高时，属于基因拷贝数发生重复的结构变异；本发明仅对重复拷贝数CN为3或4的基因位点进行分型。所述分型方法具体为：将标准化校正后的1.25<读取深度信号值<1.75范围内的基因拷贝数变异位点设为杂合重复，杂合重复的基因型为(1；2)；将标准化校正后的1.75<RD-value<2.25的CNV位点设为纯合重复，纯合重复的基因型为(2；2)。本发明中，所述显著地增高是指发生结构变异的窗口中映射匹配的reads数目显著的比参考基因组上对应的该窗口中的reads数目较多。

当所述窗口内读取深度信号与无拷贝数变异的窗口相比有显著地降低时，属于基因拷贝数发生缺失的结构变异；将所述缺失的结构变异进行分型，所述分型的方法具体为：将标准化校正后的读取深度信号值<0.10的基因拷贝数变异位点设为纯合缺失，设定纯合缺失的基因型为(0；0)；将标准化校正后的读取深度信号值为0.10<读取深度信号值<0.75的位点设为杂合缺失，设定杂合缺失的基因型为(0；1)。本发明中，所述显著地降低是指发生结构变异的窗口中映射匹配的reads数目显著的比参考基因组上对应的该窗口中的reads数目较少。

本发明中，所述显著地增高或降低具体是将读取深度信号值采用数学上的显著性统计算法进行计算。所述算法为单样本t-检验(one-sample t-test)(p<0.05)。

本发明中，所述检测标准及分型方法均基于自然群体中的平衡选择定律。

下面结合实施例对本发明提供的一种基于基因拷贝数变异位点对林木群体基因型分型方法进行详细的说明，但是不能把它们理解为对本发明保护范围的限定。

实施例1

(1)原材料的获得：从毛白杨自然分布区收集了435株个体作为研究对象。利用CTAB法提取每株个体的基因组DNA后，送至上海伯豪生物技术有限公司进行测序。测序选用Illumina的pair-end和454的mate-pair两种测序方法对毛白杨个体进行测序。

(2)比对：利用上述比对的工具、软件和算法，将每株个体得到的测序片段即reads与参考基因组序列进行比对，去除PCR重复、冗余和测序过程中引入的接头序列。

(3)统计：利用CNVnator软件和其算法统计每个具有潜在CNV位点区域的相关信息，如CNV的起始-结束坐标、CNV的长度、类型(deletion或duplication)、读取深度信号值(RD-value)、可进行数学上显著性统计的P值、确定候选CNV假阳性的q₀值等。为了提高分型结果的精确性，避免假阳性，本发明中选择p<0.01和q₀<0.5的CNV位点进行后续分析。

结合上一步骤得到的读取深度信号值，根据当前窗口内的GC含量和偏差，基于不同个体在该窗口内的读取深度信号值的数据集的中值，对该窗口的读深信号值进行数据标准化处理。

(4)校正：利用不同个体在某个窗口内的读深信号值数据集的中值对读深信号进行标准化校正，减少误差，降低假阳性。

(5)合并：由于所采用的材料毛白杨本身具有基因杂合度较高，DNA序列多态性丰富的特点，再加上测序过程中不可避免的因为测序技术、实验人员的操作等造成一定的试验误差。因此，即使同一个区域中的CNV位点也会因上述存在的问题而出现每株个体与每株个体之间的CNV起始和结束坐标并非完全一致。为了解决这一问题，本发明采取小于等于5个bin窗口(每个bin＝500bp，5个bin的长度＝500*5＝2,500bp)的算法进行合并，合并后再利用上述步骤5)、6)和7)中的检测标准和分型方法进行CNV位点基因型的检测。该步骤具体如下(以本发明中的实际例子说明)：

如，第1号染色体上发现了一个deletion类型的CNV位点，发生该CNV位点的毛白杨株数为144株，经过上述所有步骤检测到该CNV位点发生在染色体上的位置为(1260001-1265000，即1260001为起始坐标，1265000为结束坐标)，其中有30株个体以1260001为起始坐标，3株个体以1260501为起始坐标，111株个体以1261001为起始坐标。基于上述小于等于5个bin窗口的算法对该deletion类型的CNV位点的染色体位置进行合并如下：chr01：1261001-1265000。再如：同样在第1号染色体上发现了一个duplication类型的CNV，发生该CNV位点的毛白杨株数为114株，经过上述所有步骤和算法检测到该CNV发生在染色体上的位置为(1292001-1327500)，其中12株个体的结束坐标为1325000，6株个体的结束坐标为1326000，15株个体的结束坐标为1326500，66株个体的结束坐标为1327000，15株个体的结束坐标为1327500。基于上述小于等于5个bin窗口的算法对该duplication类型的CNV在染色体上的位置进行合并如下，chr01：1292001-1325000。即，如果某CNV位点存在不一致的起始坐标或结束坐标时，其起始坐标总是以发生该CNV位点最大的起始坐标作为该CNV位点的起始坐标，其结束坐标总是以发生该CNV位点的最小结束坐标作为该CNV位点的结束坐标。基于此原则，本发明对毛白杨基因组中的19条染色体上发生的CNV位点进行合并，然后进行分型。

(6)分型：利用标准化后的读深信号值对完成合并的每个CNV位点进行基因型分型。其具体操作如下：

i.对于deletion类型的CNV，将标准化校正后的RD-value<0.10的CNV位点设为纯合缺失，其基因型为(0；0)；将标准化校正后的0.10<RD-value<0.75的位点设为杂合缺失，其基因型为(0；1)；

ii.对于没有发生结构变异的基因位点，将其设为正常拷贝(CN＝2)，其基因型为(1；1)；

iii.对于duplication类型的CNV，将标准化校正后的1.25<RD-value<1.75的CNV位点设为杂合重复，其基因型为(1；2)；将标准化校正后的1.75<RD-value<2.25的CNV位点设为纯合重复，其基因型为(2；2)。

统计结果显示：本发明在毛白杨基因组的前4条染色体上共成功分型到1,628个CNV位点，其中第1号染色体上成功分型了743个CNV位点，包括523个deletion和220个duplication；第2号染色体成功分型了321个CNV位点，包括199个deletion和122个duplication；第3号成功分型了284个CNV位点，包括175个deletion和109个duplication；第4号染色体上分别成功分型了280个CNV位点，包括195个deletion和85个duplication。

由以上实施例可知，本发明提供的基于基因拷贝数变异位点对林木群体基因型的分型方法，利用读深信号进行林木自然群体中的CNV位点的基因型检测，对CNV与数量性状的关联分析具有重要意义。

以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。