CN109346130B - 一种直接从全基因组重测序数据中得到微单体型及其分型的方法 - Google Patents

Abstract

本发明公开了一种直接从全基因组重测序数据中获取微单体型及其分型的方法，包括以下步骤：获取待测个体的全部SNP位点信息，根据测序深度筛选SNP位点，获得全基因组范围的微单体型及其信息，针对待测个体上有两种及以上分型情况的二倍体类SNP，依据该SNP位点的分型情况的二项分布概率计算其属于测序误差的概率P，由此剔除潜在的重复序列，组装过程中若两个SNP之间有超过两种的连接方式，则同样依据上述方法计算是否为测序误差导致，最后根据信息熵从高到低的顺序依次对微单体型排列，即得到信息含量较高的微单体型标记。本方法可以直接从个体测序结果中获得单体型的情况，解决了非模式生物因测序误差及组装造成的影响，结果可靠。

Description

一种直接从全基因组重测序数据中得到微单体型及其分型的 方法

技术领域

本发明属于生物信息学技术领域，具体涉及一种直接从全基因组重测序数据中得到微单体型及其分型的方法。

背景技术

21世纪以前，分子标记主要是RAPD(Random Amplification of PolymorphicDNA)、AFLP(Amplified Fragment Length Polymorphism)、RFLP(Restriction FragmentLength Polymorphism)以及SSCP(Single Strand Conformation Po1ymorphism)。然而，上述四种标记已经不再广泛使用，现有的常用分子标记包括SSR(即Short Tandem Repeats，STR)以及SNP(Single Nucleotide Polymorphism)。由于SSR标记在实验过程会出现读条带不准确、筛选标记时人力物力耗费巨大等缺点，越来越多的群体分析，如群体结构分析及亲子鉴定等，选择使用SNP位点。SNP作为分子标记也有其缺点，如单个SNP位点的多态性不够高，要想达到与SSR标记同等效力，SNP标记的数量要比SSR标记多得多。现有寻找SNP标记的软件及算法主要是Stacks软件，然而它主要应用于寻找RAD-seq(Restriction siteAssociated DNA Sequencing)中的tagSNP，而RAD-seq技术比较依赖于限制酶酶切位点的数量，因此不一定能够得到全基因组上的tagSNP。

单体型(haplotype)作为连锁遗传的SNP组合，21世纪以来得到了较为迅速的发展。大量与人类疾病有关的单体型被鉴定出来，并且这也推动了GWAS(Genome-wideAssociation Study)技术的发展。现有的分析单体型的软件主要包括PHASE、HAPLOVIEW、WHATSHAP、IMPUTE2等。PHASE依据个体在每个SNP位点上的分型情况进行单体型的检测与分型，而不是直接从测序获得的短序列中获取，组装好的非模式生物基因组上的重复序列可能会对单体型的分析起到影响；HAPLOVIEW先从测序数据中分析重组率以获得SNP组合(SNPblock)，随后使用PL-EM算法(Partition-Ligation-Expectation-MaximizationAlgorithm)预测单体型，然而此方法获得的单体型是参数优化的结果；WHATSHAP使用动态规划算法(Dynamic Programming Algorithm)解决了直接从短序列获取单体型的过程中会出现的最小误差校正问题(Minimum Error Correction Problem,MEC)，然而这个算法比较适合长序列，且对测序深度有要求(不大于20X)；IMPUTE2、MACH、fastPHASE等软件则是依据HMM模型，使用已确定的单体型获取未知的单体型。上述与单体型获取有关的软件，除WHATSHAP，都与EM算法(Expectation-Maximization Algorithm)有关，此算法主要是对参数进行优化，以获取最确定的单体型分型结果，因此大部分与单体型分析有关的软件获得的结果都是优化的结果，而不是直接从测序结果中获取的确定的单体型结果，因此其结果可能还需要实验进行验证。

长度在10kbp以下的单体型称为小单体型(minihaplotype)，目前对小单体型的预测使用的计算机分析软件与单体型的一致。长度在200bp及以下的单体型称为微单体型(microhaplotype)，现今只有FLfinder一个软件专门用于分析这类单体型。FLfinder直接从测序获得的短序列获取该片段上存在的单体型及其分型情况，然而忽视了测序过程中可能存在的测序误差，以及非模式生物由二代测序组装的参考基因组上的重复片段的问题，因此获得的结果不准确。

微单体型作为一种新型的分子标记，目前还没有得到十分广泛的使用，用于预测和分析此类分子标记的软件也没有很多，本发明旨在提供一种直接从全基因组重测序数据中得到微单体型及其分型的方法，能够在全基因组范围内精确获取微单体型，以便后续使用此类分子标记进行其他生物学研究。

发明内容

本发明的目的在于提供一种直接从全基因组重测序数据中得到微单体型及其分型的方法，其优点在于考虑到了全基因组重测序中可能出现的测序误差以及基因组组装过程中合并的重复序列导致的影响。

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

S1)对待测个体进行全基因组重测序；

S2)使用比对软件(如BWA等)，将测序结果与相同物种或近似物种(依据研究对象进行选取)的参考基因组进行比对；

S3)对得到的比对文件使用检测SNP的软件(如GATK等)进行SNP检测；

S4)对各SNP进行测序深度的筛选：分别统计各SNP位点测序深度，并依据所得的测序深度作出的箱型图，选取全部个体的集中区域的最小值作为测序深度的阈值，并筛选深度大于该阈值的SNP位点；统计各SNP位点在全部个体中的分型情况，若完全一致，则舍弃该SNP位点；

进一步地，为了获得由信息量较高的SNP位点组成的微单体型标记，还可以各SNP位点的MAF值进行筛选，MAF的阈值一般设为0.2；

S5)获得全基因组范围内的微单体型：全基因组搜索两个SNP之间距离小于设定阈值的全部SNP组合，所述的阈值设为0-999bp(依据不同测序平台而进行设置，如使用Illumina测序，获得的序列片段长度为150bp，因此该值可根据需要设为不低于150bp的任一数值)，若连续若干SNP，其两两之间距离均小于该阈值，则将这些SNP作为整体输出；若无其他SNP与该SNP距离小于该阈值，则该SNP单独输出，上述“距离”全部指物理距离；

S6)获得各单体型的具体SNP位点信息：通过S5)获得的各SNP组合，筛选其中SNP数量大于等于2的SNP组合作为潜在的微单体型，对每一个微单体型，依据S3)中获得的SNP信息文件(VCF文件)获得覆盖该微单体型片段的全部SNP位点信息，并按从小到大的顺序(碱基位置的顺序)将各SNP位点一一列举，加上所在Scaffold的ID、SNP个数以及该片段长度等三个信息组成该微单体型的信息；

S7)使用下述原理和方法进行微单体型的检测：

S7.1)获得全部个体覆盖某标记的全部测序短序列(以下简称“短序列”)信息；

S7.2)依据该标记的SNP信息获取待测个体在这些SNP位点上的碱基情况，并将这些SNP位点分成“一致的”与“二倍体的”两类，“一致的”指该SNP在待测个体上只有一种分型情况，“二倍体的”指该SNP在待测个体上有两种及以上分型情况(图2)；测序误差筛选方法如下：对每个“二倍体”类SNP依据各分型情况的二项分布概率计算其属于测序误差的概率P：

其中n为覆盖该位点的短序列总数，m为上述短序列中某一分型的数量，概率p为测序误差概率，该值设为0.001-0.02(测序误差概率p值依据不同测序平台的测序误差设置，如测序平台为Illumina，则依据文献报告将此值设置为0.008或0.01，其默认值为0.01)，若计算所得概率大于0.01，则认为该分型是由测序误差导致，并将该分型剔除，并接受该SNP位点，以用于后续微单体型的组装，如对于某一个体的某一SNP位点，存在两种分型即“A”和“T”,其中“A”有一条短序列支持，“T”有30条，则通过(1)式可以将“A”剔除，保留“T”，该SNP位点在此个体中以“T”分型存在；若计算概率后该位点依旧含有两种以上分型情况，则认为该位点所在微单体型片段可能是重复序列，为降低分析误差从而将该片段从全部微单体型中去除；

S7.3)组装过程主要针对于“二倍体”类SNP进行，首先对“二倍体”类SNP位点的位置进行筛选，若两个SNP位点之间的距离大于设定的阈值(此值由不同测序平台指定，如Illumina测序平台获得的短序列长度为150bp，则该阈值上限为150，此阈值范围可设置为50-900bp)，则认为肯定无短序列覆盖，遂将该片段在此处切断，分成两个子标记；若不存在上述情况，则继续依赖于该“二倍体”类SNP与其后一个“二倍体”类SNP的关系，若二者没有短序列覆盖，则无法进行组装，于是将该标记在此处打断，分成两个子标记，若打断后依旧出现相邻两个SNP位点没有短序列覆盖，则在这两个SNP位点处打断，以此类推；若有短序列将二者覆盖，则依据其覆盖的短序列个数进行归类(统计短序列的覆盖数量)(图3)，若这两个SNP之间有两种以上连接方式，则首先依据各组合方式的二项分布概率计算是否为测序误差导致，使用的参数与计算方法和(1)一致(由于是两个SNP位点连接情况的筛选，因此此值与测序平台无关，此值越高，对SNP连接方式的过滤性能越强，建议设置范围为0-0.02，默认为0.01)，若经计算依旧含有两种以上连接方式，则认为该片段位于重复片段上，从而将该片段去除；由于某一条链可能在某两个SNP位点上没有短序列覆盖，最终组装后只得到一条链，则另一条链依据之前的“二倍体”类SNP的结果，完成组装；

S7.4)上述组装过程最终会得到两条确定的只含有“二倍体”类SNP组成的链，随后将“一致”类SNP与已组装好的两条链按照参考基因组上的顺序进行组合(图2)，即可得到该个体在该标记上的二倍体分型情况；

S7.5)采用上述方法获取全部个体的微单体型分型结果，并对每个标记(即微单体型)在全部个体的分型情况进行信息含量的计算，计算方式为：

其中，N表示该SNP组合中等位基因的个数，p_i表示第i个等位基因在群体中的基因频率，H表示该组合在群体中的信息熵；

S7.6)依据计算结果，按信息熵从高到低的顺序依次排列，即可得到信息含量较高的微单体型标记。

步骤S7)使用了二项分布概率对测序结果中可能存在的测序误差进行了计算与排除，也对某两个SNP之间的多种连接方式进行了测序误差方面的排除，若经过上述排除后某两个SNP位点之间依旧含有两种以上连接方式，则认为它们位于基因组上的重复片段，但由于使用二代测序进行组装时合并到一起的区域，从而剔除掉该标记。

本发明的有益效果：可以直接从个体测序结果中获得单体型的情况，解决了非模式生物使用二代测序技术组装而成的参考基因组中重复序列的影响以及测序误差可能造成的影响，可以获得可靠的分子标记，由于是直接从测序结果中获得的微单体型，而不是由参数优化得到，因此，不需要进行实验验证。同时，本发明还可以获得信息含量较高的分子标记，可用于亲子鉴定、亲缘关系鉴定、群体结构鉴定等方面。

附图说明

图1为微单体型标记筛选流程图。

图2为SNP分类依据原理图(举例)，其中左上角第一行为参考基因组，其中加粗的碱基为在该个体中为纯合子的SNP位点，其下方为覆盖该片段的短序列，加粗和下划线的碱基为该短序列中该SNP位点的碱基情况，右上角为各短序列SNP碱基序列情况，下划线的碱基为该个体中纯合子的SNP位点，未覆盖的SNP位点使用“N”代替。

图3为某微单体型组装的原理图(举例)，其中圆圈表示某一个SNP位点上的碱基分型情况，箭头表示依据短序列中匹配数量得到的前后两SNP位点连接方式，下方的碱基顺序为通过上述连接方式得到的“二倍体的”SNP的连接结果。

图4为待测个体全部SNP位点测序深度分布箱型图。横坐标为个体ID，由于R软件作图限制，F1表示个体1jing，M2表示个体2jing，F13表示个体m_F13；纵坐标表示各SNP测序深度，图中粗实线为测序深度为18的阈值。

具体实施方式

下面结合图1对本发明所述的直接从全基因组重测序数据中获取微单体型及其分型的方法进行说明。

S1)选取本实验室2017年用于繁殖的草鱼个体(雌性：m_F13、F40、2jing，雄性：M10、1jing)作为实验样品，并使用全基因组重测序，测序平台为Illumina X10，测序深度为30X。

S2)首先使用BWA软件或SAMtools软件对上述五个个体的全基因组重测序数据进行比对，所使用的参考基因组为草鱼基因组，使用默认参数。

S3)使用GATK(https://software.broadinstitute.org/gatk/)进行SNP检测，共检测到4,955,011个SNP(SNP信息保存在VCF文件中)，筛选条件为QUAL<30.0，SOR>10.0，QD<2.0，FS>200.0，InbreedingCoeff<-0.8。

S4)使用SAMtools(http://samtools.sourceforge.net/)depth统计检测到的SNP测序深度，使用R软件(https://www.r-project.org/)作出箱型图获得测序深度阈值(依据统计学中箱型图相关概念，其盒子所在范围表示数据较为集中的区域，因此获取全部个体的集中区域的最小值)。从图4中可以看出，其大部分SNP位点测序深度都在18X以上，因此测序深度阈值设置为18X。

为了获得由信息量较高的SNP位点组成的微单体型标记，还可以进一步根据各SNP位点的MAF值进行筛选，MAF值的阈值可依据文献进行设置，一般设为0.2，如若获得的标记将用于亲子鉴定，则需要存在个体在每个SNP位点的分型情况与其他个体不同即可，此时阈值可设置为0。本实施例为了尽可能获得能将上述5尾个体分开的微单体型，MAF值设置为0，随后通过VCF文件中各SNP位点在个体中分型情况的统计，过滤掉在全部个体中分型一致(如都为“1/1”)的位点，共获得2,915,134个SNP。从比对得到的VCF文件中判断各SNP位点在每个个体上的支持的短序列数量，若总数量与分型对应的数量之和的差值大于2，则认为可能这个位点所在的微单体型是重复序列，因此记录该SNP位点(本实施例中共2,736个)，以用于后续对微单体型的筛选；其他过滤后的SNP信息记录在文件中。

S5)设置同一微单体型中两个SNP位点间距为100bp以内(依据不同测序方式而进行设置，如使用Illumina测序，获得的序列片段长度为150bp，该值可根据需要设为低于150bp的任一数值)，由于结果文件中含有只有一个SNP位点的“微单体型”，因此对获得的“微单体型”进行筛选，共得到572,101个微单体型。

S6)获取得到的微单体型信息：为了获得信息含量较高的微单体型位点，我们对获得的微单体型依据其所含有SNP位点的个数由多到少排序，取排名前100的微单体型位点，储存于文件中进行分型分析。

S7)对于上述100个微单体型中的每一个微单体型，将其所在scaffold的ID与其起始位点和终止位点作为该微单体型的具体区域，结果储存于文件中；获取SNP碱基分型情况方法有三种，可以依据VCF文件中已确定的分型情况确定该位点可能的碱基，也可以依据GATK软件(https://software.broadinstitute.org/gatk/)获取SNP信息的方法获得该位点的碱基，还可以完全依据短序列上在该SNP位点上各碱基分型的数量确定该位点可能的碱基，本实施例获取各SNP位点碱基分型情况的方法为第三种，即依据每个样本的短序列中各SNP位点上不同碱基出现的次数，并进行测序误差的筛选，所使用的测序过程导致的误差参数p值为0.01，用于筛选的测序误差概率为0.01，共得到116个微单体型。若不进行测序误差筛选，共可获得188个微单体型。现列举三个微单体型片段进行实施案例分析，其中“：”之前的为所在Scaffold，“-”之前为起始位点，其后为终止位点，列举碱基数量时，一个“[]”为某一个SNP位点，其中引号“‘’”内表示碱基分型，“-”后为该分型对应的数量：

①片段CI01000354:1612160-1612353:

该片段所含有的SNP所在位点如下：1612160,1612168,1612172,1612271,1612313,1612345,1612353。以个体2jing为例，其短序列中各位点上碱基数量如下：[‘A’-31，‘G’-1],[‘C’-31],[‘C’-30],[‘A’-33],[‘G’-33],[‘C’-36],[‘A’-37]。第一个位点上经过(1)式计算，“G”碱基分型由测序误差导致，因此去除该分型，第一个位点的最终分型结果为[‘A’-31]。综合上述7个位点，均为“一致类”SNP位点，因此直接进行组装得到个体2jing在该片段下的微单体型为“ACCAGCA”。

②片段CI01000149:246126-247776及CI01000149:247661-247723：

该片段中含有的SNP位点如下：246126,246144,246159,246186,246270,246464,246524,246552,246554,246565,246568,246608,246635,246638,246639,246717,246752,246757,246863,246951,246968,247034,247047,247107,247108,247134,247171,247226,247227,247242,247293,247309,247353,247387,247447,247453,247481,247505,247590,247661,247703,247705,247723,247742,247762,247776。以个体F40为例，各SNP位点碱基分型如下：[‘A’-11,‘C’-2],[‘C’-15],[‘T’-17],[‘T’-16],[‘C’-29],[‘G’-11,‘T’-8],[‘C’-13,‘T’-7],[‘C’-13,‘G’-6],[‘A’-5,‘G’-13],[‘C’-5,‘T’-13],[‘C’-13,‘T’-5],[‘A’-13,‘G’-7],[‘C’-18,‘T’-1],[‘C’-19],[‘A’-12,‘G’-7],[‘C’-1,‘T’-22],[‘T’-20],[‘T’-21],[‘A’-19,‘G’-11],[‘C’-27],[‘C’-25],[‘C’-23],[‘A’-20],[‘A’-9,‘T’-10],[‘C’-9,‘T’-10],[‘T’-14,‘G’-11],[‘C’-28],[‘C’-14,‘T’-19],[‘A’-33,‘C’-1],[‘C’-21,‘G’-16],[‘T’-33],[‘T’-35],[‘G’-40],[‘C’-40],[‘A’-24,‘G’-19],[‘C’-19,‘T’-25],[‘C’-23,‘T’-18],[‘T’-42],[‘A’-46,‘C’-1],[‘G’-44],[‘G’-34],[‘T’-34],[‘T’-35],[‘G’-19,‘C’-13],[‘A’-32,‘G’-1],[‘A’-38,‘G’-1]。

对上述SNP位点进行测序误差过滤，可得最终各SNP位点碱基分型如下：

[‘A’-11,‘C’-2],[‘C’-15],[‘T’-17],[‘T’-16],[‘C’-29],[‘G’-11,‘T’-8],[‘C’-13,‘T’-7],[‘C’-13,‘G’-6],[‘A’-5,‘G’-13],[‘C’-5,‘T’-13],[‘C’-13,‘T’-5],[‘A’-13,‘G’-7],[‘C’-18],[‘C’-19],[‘A’-12,‘G’-7],[‘T’-22],[‘T’-20],[‘T’-21],[‘A’-19,‘G’-11],[‘C’-27],[‘C’-25],[‘C’-23],[‘A’-20],[‘A’-9,‘T’-10],[‘C’-9,‘T’-10],[‘T’-14,‘G’-11],[‘C’-28],[‘C’-14,‘T’-19],[‘A’-33],[‘C’-21,‘G’-16],[‘T’-33],[‘T’-35],[‘G’-40],[‘C’-40],[‘A’-24,‘G’-19],[‘C’-19,‘T’-25],[‘C’-23,‘T’-18],[‘T’-42],[‘A’-46],[‘G’-44],[‘G’-34],[‘T’-34],[‘T’-35],[‘G’-19,‘C’-13],[‘A’-32],[‘A’-38]。

以“[]”表示某一个SNP位点，其中“,”之前的内容与上同，“,”后表示该SNP位点所在位置。“一致类”的SNP位点为：[‘C’-15,246144],[‘T’-17,246159],[‘T’-16,246186],[‘C’-29,246270],[‘C’-18,246635],[‘C’-19,246638],[‘T’-22,246717],[‘T’-20,246752],[‘T’-21,246757],[‘C’-27,246951],[‘C’-25,246968],[‘C’-23,247034],[‘A’-20,247047],[‘C’-28,247171],[‘A’-33,247227],[‘T’-33,247293],[‘T’-35,247309],[‘G’-40,247353],[‘C’-40,247387],[‘T’-42,247505],[‘A’-46,247590],[‘G’-44,247661],[‘G’-34,247703],[‘T’-34,247705],[‘T’-35,247723],[‘A’-32,247762],[‘A’-38,247776]。

“二倍体类”的SNP位点为：[‘A’-11,‘C’-2,246126],[‘G’-11,‘T’-8,246464],[‘C’-13,‘T’-7,246524],[‘C’-13,‘G’-6,246552],[‘A’-5,‘G’-13,246554],[‘C’-5,‘T’-13,246565],[‘C’-13,‘T’-5,246568],[‘A’-13,‘G’-7,246608],[‘A’-12,‘G’-7,246639],[‘A’-19,‘G’-11,246863],[‘A’-9,‘T’-10,247107],[‘C’-9,‘T’-10,247108],[‘T’-14,‘G’-11,247134],[‘C’-14,‘T’-19,247226],[‘C’-21,‘G’-16,247242],[‘A’-24,‘G’-19,247447],[‘C’-19,‘T’-25,247453],[‘C’-23,‘T’-18,247481],[‘G’-19,‘C’-13,247742]。

依据S7.3)，对上述“二倍体”类SNP位点进行筛选，位于246126的SNP与位于246464的SNP间距大于150bp，位于246639的SNP与位于246863的SNP间距大于150bp，位于247107的SNP与位于246863的SNP间距大于150bp，位于247242的SNP与位于247226的SNP间距大于150bp，位于247447的SNP与位于247242的SNP位点间距大于150bp，位于247742的SNP与位于247481的SNP间距大于150bp。遂将该微单体型区域依据其含有的全部SNP位点进行打断，分为246126-246270，246464-246757，246863-247047，247107-247387，247447-247590，247661-247723，247742-247776。

现以247661-247723片段(SNP位点为247661,247703,247705,247723)与个体1jing为例，覆盖该片段的短序列经过SNP位点的选择后，得到的各短序列对应的微单体型序列如下：TAAT,GGTT,GGTT,GGTT,GGTT,GGTT,GGTT,GGTT,GGTT,TAAT,GGTT,GGTT,NAAT,NGTT,NGTT,NAAT,NAAT,GGTN,NGTT,TNNN,GNNN,GNNN,TNNN,GNNN,GNNN,TNNN,GNNN,NNNT,TNNN,NNNT,GNNN。统计上述片段中各SNP位点碱基分型情况如下：[‘G’-17,‘T’-6],[‘G’-14,‘A’-5],[‘A’-5,‘T’-14],[‘T’-20]。依据S7.2)方法进行测序误差筛选后，将它们分成两类，其中“一致”类：[‘T’-20，247723]；“二倍体”类：[‘G’-17,‘T’-6，247661],[‘G’-14,‘A’-5,247703],[‘A’-5,‘T’-14,247705]。

依据S7.3)中的方法，对“二倍体”类SNP进行连接，其中以“[]”表示某两个SNP连接，以“‘X-Y’”表示碱基X与Y连接，“：”后为各连接方法的数量，其连接结果为[‘T-A’:2,‘G-G’:11],[‘G-T’:14,‘A-A’:5]。则“二倍体”类SNP组装后的结果为：[‘TAAN’,‘GGTN’]。依据S7.4)的方法，将“二倍体”类组装结果与“一致”类组装结果结合，得到最终组装结果为[‘TAAT’,‘GGTT’]。

③片段CI01000358:825019-825060及个体1jing：

该片段SNP位点为825019,825025,825027,825060，短片段序列为CAAT,CGCT,CAAT,CGCT,CGCT,CGCT,CAAT,CGCT,CGCT,CGCT,CGCT,CAAT,CAAT,CAAT,NNNT,NNNT,NNNT,CGCN,CAAN,CGCN,CAAN,TACN,CAAN,TACN,CGCN,TACN,CGCN,CGNN。使用S7.2)中的方法，统计各SNP位点碱基分型情况如下(表示方式与②相同)：[‘C’-22,‘T’-3],[‘G’-13,‘A’-12],[‘C’-15,‘A’-9],[‘T’-17]。经过测序误差筛选后，进行分类，“二倍体”类SNP为：[‘C’-22,‘T’-3，825019],[‘G’-13,‘A’-12,825025],[‘C’-15,‘A’-9,825027]；“一致”类SNP为[‘T’-17，825060]。

使用S7.3)中的方法进行连接，通过短序列可以得到(表示方法与②相同)：[‘C-G’:13,‘C-A’:9,‘T-A’:3]。对于第一个“二倍体”类SNP与第二个“二倍体”类SNP，其连接方式经过测序误差筛选后依旧含有三种，因此认为该片段可能位于重复序列中，遂舍弃该片段。

Claims (2)

1.一种直接从全基因组重测序数据中获取微单体型及其分型的方法，包括以下步骤：

S1)对待测个体进行全基因组重测序；

S2)使用比对软件将测序结果与相同物种或近似物种的参考基因组进行比对；

S3)使用检测SNP的软件进行SNP位点检测，获得SNP信息文件；

S4)根据测序深度筛选SNP位点：分别统计待测个体的SNP位点测序深度，根据所得的测序深度作出箱型图，选取全部个体的集中区域的最小值作为测序深度的阈值，并筛选深度大于该阈值的SNP位点；统计各SNP位点在全部个体中的分型情况，若完全一致则舍弃该SNP位点；

S5)获得全基因组范围内的微单体型：全基因组搜索两个SNP之间距离小于设定阈值的全部SNP组合，所述的阈值为0-999bp；

S6)获得潜在微单体型的具体SNP位点信息：通过S5)获得的各SNP组合，筛选其中SNP数量大于或等于2的SNP组合作为潜在的微单体型，并依据S3)中获得的SNP信息文件获得覆盖该微单体型片段的全部SNP相关信息，包括全部SNP具体位置、所在Scaffold的ID、该片段的SNP个数以及该片段长度；

S7)按照下述方法进行微单体型的检测：

S7.1)获得全部个体覆盖某标记的测序短序列信息；

S7.2)根据待测个体在SNP位点上的碱基情况，将SNP位点分成一致与二倍体两类，一致类指该SNP在待测个体上只有一种分型情况，二倍体类指该SNP在待测个体上有两种及以上分型情况，筛选有两种及以上分型情况的SNP位点；并依据SNP位点的分型情况的二项分布概率计算其属于测序误差的概率P，计算公式如下：

其中n为覆盖该位点的短序列总数，m为上述短序列中某一分型的数量，概率p为测序误差概率，该值设为0.001-0.02，若计算后概率P大于0.01则认为该分型是由测序误差导致，将该分型剔除并接受该SNP位点，以用于后续微单体型的组装，若计算概率后该位点依旧含有两种以上分型情况，则认为该位点所在微单体型片段是潜在的重复序列，为降低分析误差从而将该片段从全部微单体型中去除；

S7.3)首先对二倍体类SNP位点的位置进行筛选，若两个SNP位点之间的距离大于设定的阈值，所述的阈值范围为50-900bp，则将该片段在此处切断，分成两个子标记；若不存在上述情况，则继续依赖于该二倍体类SNP与其后二倍体类SNP的关系，若二者没有短序列覆盖，则无法进行组装，于是将该标记在此处打断，分成两个子标记，若打断后依旧出现相邻两个SNP位点依旧没有短序列覆盖，则在这两个SNP位点处打断，以此类推；若有短序列将二者覆盖，则统计各连接方式的短序列个数，若这两个SNP之间有超过两种的连接方式，则首先依据各组合方式的二项分布概率计算是否为测序误差导致，使用的参数与计算方法和公式(1)一致，若经计算依旧含有超过两种的连接方式，则认为该片段位于重复片段上，从而将该片段去除；若某一条链在某两个SNP位点上没有短序列覆盖，最终组装后只得到一条链，则另一条链依据之前的二倍体类SNP的结果，完成组装；

S7.4)上述组装过程最终会得到两条确定的只含有二倍体类SNP组成的链，随后将一致类SNP与已组装好的两条链按照参考基因组上的顺序进行组合，即可得到待测个体在该标记上的二倍体分型情况；

S7.5)采用上述方法获取全部个体的微单体型分型结果，并对每个微单体型在全部个体的分型情况进行信息含量的计算，计算方式为：

其中，N表示该SNP组合中等位基因的个数，p_i表示第i个等位基因在群体中的基因频率，H表示该组合在群体中的信息熵；

S7.6)按信息熵从高到低的顺序依次对微单体型排列，即可得到信息含量较高的微单体型标记。

2.根据权利要求1所述的直接从全基因组重测序数据中获取微单体型及其分型的方法，其特征在于，所述的S4)中进一步根据各SNP位点的MAF值进行筛选，MAF的阈值设为0.2。

Images