CN105787294B - 确定探针集的方法、试剂盒及其用途 - Google Patents

Abstract

本发明公开一种确定探针集的方法，所述探针集包含多个探针，所述方法包括：获取DNA序列集，所述DNA序列集包含多条长度为K₀的DNA序列，获取所述DNA序列集包括，从长度为L的参考序列的一端的第i个核苷酸开始，沿另一端方向拷贝所述参考序列的K₀个连续核苷酸为一条DNA序列，i依次取{1，2，…，(L‑K₀)，(L‑K₀+1)}中的数值；筛选所述DNA序列集，以获得所述探针，其中包括，过滤掉与所述参考序列的匹配位置个数大于2的DNA序列；其中，200nt≥K₀≥50nt。本发明还公开一种试剂盒、该试剂盒的用途以及一种确定物种丰度的方法。

Description

确定探针集的方法、试剂盒及其用途

技术领域

本发明涉及生物技术领域，具体的，本发明涉及一种确定探针集的方法、一种试剂盒、试剂盒的用途以及一种确定物种丰度的方法。

背景技术

特定环境中的生物多样性状况是了解该生态环境健康程度及可持续发展能力的重要参考，也是解决管理过程中出现的问题的重要依据[1-3]。时至今日，迅速发展起来的第二代测序平台越来越多的应用到包括节肢动物、线虫等各种类群[4-9]及环境DNA中植物、真菌和蚯蚓的生物多样性研究中[10-12]。此类研究多利用基于PCR产物测序的技术，此技术有着难以克服的缺点，即物种偏向性[7，13]。

线粒体因其独特的系统发育历史常被用作系统发育研究的重要分子标记物。近期，许多研究为了进行物种划分和系统发育推演，都致力于利用全基因鸟枪法构建线粒体基因文库[14-17]。理论上，此方法与传统的叠瓦式PCR或大片段PCR相比，效率更高且对DNA质量要求更低[18]，但是离心技术难以回收线粒体基因导致大量测序数据的浪费(利用率仅0.5％)束缚了此技术的推广。因此，探索新理论、开发新技术以针对性地富集线粒体基因组可以加速推广以生物多样性为基础生物监测的应用。

人工合成的寡核苷酸探针对DNA进行杂交捕获的技术广泛应用于分子生物学、分子遗传学及临床医学等研究领域，进行基因诊断及核酸的定性定量检测。一般情况下，探针的杂交捕获对象与探针设计时使用的参考数据库物种是统一的，同时还可根据不同目的进行调整，如人外显子捕获[19]等，但也有研究者利用探针捕获非设计来源物种的DNA，有研究表明人类外显子捕获芯片对非人灵长目外显子的部落效率高达约95％[20]。也有报道称在高度分化的物种间可捕获保守的直系同源基因，但捕获效率较低[21，22]。杂交捕获在生物多样性研究中的应用受到限制主要由于探针效率低及缺少对异源序列捕获条件的摸索。

发明内容

本发明旨在至少一定程度解决上述问题之一或者至少提供一种商业选择。

依据本发明的一方面，提供一种确定探针集的方法，所述探针集包含多个探针，所述方法包括：获取DNA序列集，所述DNA序列集包含多条长度为K₀的DNA序列，获取所述DNA序列集包括，从长度为L的参考序列的一端的第i个核苷酸开始，沿另一端方向拷贝所述参考序列的K₀个连续核苷酸为一条DNA序列，i依次取{1，2，…，(L-K₀)，(L-K₀+1)}中的数值；筛选所述DNA序列集，以获得所述探针，其中包括，过滤掉与所述参考序列的匹配位置个数大于2的DNA序列；其中，L＞K₀，200nt≥K₀≥50nt。所说的过滤掉与参考序列匹配的位置个数大于1的DNA序列，是指弃去能够匹配到参考序列上两个不同位置或者两个以上不同位置的DNA序列，这里的匹配或者匹配到是指与参考序列上的某段序列一致，这里的一致可以是严格的一致，即一模一样，表现出与参考序列上的某段序列上的每个碱基和碱基顺序都相同，一致也可以是非严格绝对一致，在本发明的一个具体实施方式中，当与参考序列的碱基错配不超过3个也属于“一致”。参考序列指已知序列，可以是预先获得的目标个体所属生物类别中的任意的参考模板，例如，若目标个体是人类，参考序列可选择NCBI数据库提供的HG19，进一步地，也可以预先配置包含更多参考序列的资源库，在进行序列比对前，先依据目标个体的组成个数、地域等因素选择或是测定组装出更接近的序列来作为参考序列，有助于获得更全面更有针对性的探针集。

在本发明的一个具体实施方式中，所述参考序列为线粒体参考序列，并且包含多种昆虫的线粒体参考序列，在获得覆盖所有多种线粒体参考序列的、数目很大的DNA序列后，在筛选DNA序列集时，将所述DNA序列集比对到所说的参考序列上，设置允许的错配碱基不大于3，即比对上(匹配到)参考序列的每条DNA序列最多允许有3个碱基错配(mismatch)，若一条DNA序列中有超过3个碱基发生错配，则视为该条DNA序列无法比对到参考序列。比对可以利用SOAP(Short OligonucleotideAnalysis Package)，bwa等软件进行，本发明方法对此不作限制。在获得能够全面覆盖目标区域的DNA序列集后，筛选所述DNA序列集从而获得的探针具有高特异性。

在本发明的一个具体实施方式中，筛选所述DNA序列集还包括：计算过滤后的各条DNA序列的长度为K₁的单元序列的频数，过滤掉频数超出预定范围的DNA序列。所述预定范围是依据参考序列的单元序列的频数中位数来确定的，所述参考序列的单元序列的频数中位数通过以下确定：首先，将所述参考序列拆分成多个长度为K₁的单元序列，所述单元序列通过以下获得——从长度为L的参考序列的一端的第j个核苷酸开始，沿另一端方向拷贝所述参考序列的K₁个连续核苷酸为一条单元序列，j依次取{1，2，…，(L-K₁)，(L-K₁+1)}中的数值，接着，依据序列相同的单元序列分到相同的组，对所述单元序列进行分组，获得q个单元序列组，将单元序列组r包含的单元序列的个数表示为A_r，A_r为所述参考序列对单元序列组r中的单元序列的频数，取{A₁，A₂，…，A_q}中的中位数为所述参考序列的单元序列的频数中位数，其中，4*K₁≤K₀，r为单元序列组的编号，1≤r≤q，q≤L-K₁+1，A_r≥1。在本发明的一个具体实施方式中，所述DNA序列的单元序列的频数＝匹配到该DNA序列的单元序列的数目/包含匹配到该DNA序列的单元序列的单元序列组的个数。频数也称“次数”，为对总数据按某种标准进行分组，统计出各个组内含个体的个数。预定范围可以依据经验或者预使所得的探针集的使用环境和目的，包括捕获的混合物种的种类数目、混合物种目的核酸的序列复杂度等来调整，例如，总目标区域大(例如物种多和/或每个物种的目标区大)和/或不同物种的目标区域序列相似程度低，可以设置预定范围为上下浮动较大的比例的频数中位数，比如0.7～1.3倍的参考序列的单元序列的频数中位数，预定范围越大，会有较多的DNA序列经筛选后得以保留，在本发明的一个具体实施方式中，所述预定范围为0.8-1.2倍的所述参考序列的单元序列的频数中位数，使得筛选得的DNA序列集数目明显减少，即使探针集中的探针数目减少但其覆盖度没有降低且特异性提高。通过该方式的利用单元序列频数筛选得的DNA序列集中的DNA序列的组成碱基比例更为相近，有利于确定的探针集在同一反应条件下进行捕获和/或洗脱。

在本发明的一个具体实施方式中，筛选所述DNA序列集还包括：将所述参考序列划分成不连接的大小为m的窗口，相邻窗口间距为n，将过滤后的DNA序列与所述参考序列比对，获得各个窗口的比对结果，当一个窗口的比对结果包含的DNA序列数目大于等于1且包含的DNA序列中至少有一条的GC含量为40％～60％时，选择其中的GC含量为40％～60％的一条DNA序列，剔除该窗口比对结果包含的其它DNA序列，而当一个窗口的比对结果包含的DNA序列数目大于等于1且不包含GC含量为40％～60％的DNA序列时，选择其中的GC含量最接近60％的一条DNA序列并且拷贝该条DNA序列一次，剔除该窗口比对结果包含的其它DNA序列，其中，K₀+10nt≥m+n≥K₀，m≥n。GC富含区容易发生非特异性反应，经过该筛选，使留下来的DNA序列在不降低覆盖度的同时进一步减少数目及提高整体特异性，而且筛选得的DNA序列的整体的GC含量处于相同的40％-60％这一水平，利于获得的探针集在同一反应条件同一反应系统中一起使用。在本发明的一个具体实施方式中，当一个窗口的比对结果包含的DNA序列数目为多条且其中的GC含量为40％～60％的DNA序列也是多条，在选择GC含量为40％～60％的其中一条DNA序列时，还进一步考虑DNA序列的熔解温度，优先选择GC含量符合要求且熔解温度最靠近平均熔解温度的那条DNA序列。平均熔解温度指过滤后的所有DNA序列的熔解温度温度的平均值，将这里的DNA序列当成双链DNA序列来计算其熔解温度(Tm)，Tm的计算可利用已知公式或软件，例如，使用Santa Lucia或者Bolton and McCarthy公式来计算。这样，进一步使筛选得的能代表/覆盖各个窗口的DNA序列的GC含量和GC分布更加相近，使所确定的探针集包含尽量少的探针且整体具高特异性，且适于在同一条件体系下使用，使能在不降低对目标区域的覆盖能力的同时具有高特异性和高有效性。

在本发明的一个具体实施方式中，所述方法还包括：反转录筛选后的DNA序列，获得RNA序列，以使获得的探针集包含RNA探针。在本发明的一个具体实施方式中，在反转录之前，先扩增筛选后的DNA序列以使获得足够的DNA探针数。本领域技术人员可以理解，通过逆转录本发明上述的或者任一具体实施方式方法获得的DNA序列而获得的RNA探针也具有上述DNA探针集的优点和特点。在本发明的一个具体实施方式中，使所述RNA探针游离于溶液中，即获得液相探针/芯片。根据目前已知的探针设计和制备成本计算，液相捕获芯片的主要成本在探针合成这一步，一张96K芯片的成本约为19000元，每个芯片产生的探针可以做2,000此捕获实验，即每个杂交反应和样品的探针成本不高于10元，即利用本发明方法获得的探针集仅花费20元对混合核酸进行捕获之后再测序，在相同测序数据量的前提下得到的线粒体基因组数据量是传统方法得到线粒体基因组数据量的近一百倍，对于测序成本有极大的降低。

依据本发明的另一方面，本发明提供一种试剂盒，其包含探针集，所述探针集能够特异性识别线粒体DNA，所述探针集利用上述本发明的一方面或者任一具体实施方式的方法确定。前述对通过本发明一方面的方法确定的探针集的优点和技术特征的描述，也适用于本发明这一方面的试剂盒，在此不再赘述。

依据本发明的再一方面，本发明提供上述试剂盒在捕获线粒体基因组序列中的用途。前述的对本发明一方面的试剂盒的优点和技术特征的描述，也适用于本发明这一方面的试剂盒的用途，在此不再赘述。

依据本发明的又一方面，本发明提供一种确定物种丰度的方法，所说的物种丰度是指该物种在一个群体中的相对丰度，所述方法包括：(1)获得多个物种的混合线粒体DNA的测序数据，其中包括，利用前述本发明的试剂盒捕获多个物种的混合核酸中的混合线粒体DNA，所述测序数据包括多个读段；(2)将所述测序数据与各个所述物种的线粒体参考序列比对，获得各个物种的比对结果；(3)基于所述物种的比对结果中的读段数目，确定该物种的丰度。在本发明的一个具体实施方式中，所说的线粒体参考序列包含多个物种的已知线粒体序列，利用测序数据中比对到某物种线粒体参考序列的读段的数目来代表该物种在这个群体中的相对丰度。由于利用前述本发明一方面的试剂盒来获得测序数据，利用试剂盒中的探针集高特异性、高效率捕获线粒体DNA序列，使得本发明这一方法中的测序数据质量高、利用率高，使得利用本发明的这一方面方法确定的物种丰度能够准确反映物种的实际丰度，提供一种监测生物多样性、检测生物量的简便且低成本的方法。

附图说明

本发明的上述和/或附加的方面和优点，结合下面附图对实施方式的描述将变得明显和容易理解，其中：

图1是本发明的一个具体实施方式中的一次性获取多个物种线粒体基因组的步骤流程图；

图2是本发明的一个具体实施方式中的蛋白编码基因在分类阶元里的平均遗传距离和AT含量示意图；

图3是本发明的一个具体实施方式中的个体相对丰度的预估值的示意图。

具体实施方式

结合本发明的方法和/或试剂盒，基于混合样品的先验物种分类知识，结合第二代高通量测序以及大数据的信息分析，能够一次性得到多个物种的线粒体基因组，包含的步骤流程一般如图1所示，详述如下：

1.捕获探针的设计

为了保证探针的覆盖度，我们选取了来自1KITE项目的线粒体数据。1KITE是千种昆虫转录组进化(1K Insect Transcriptome Evolution)项目的简称，此项目旨在研究覆盖昆虫纲所有目的1000种昆虫的转录组(www.1kite.org)。该示例中的线粒体参考序列利用第一批379个物种的1KITE数据获得的，线粒体参考序列共包含2,553条平均长度1,902bp的组装序列，利用全部这些组装序列进行随后的探针设计和分析。探针设计整个过程可分为以下四步：

(1)探针长度选择，该步是可选择的：首先，把线粒体组装序列分为长度为65、73和93的三个DNA序列集(在该示例中也称为探针集)并计算每个组的理论熔解温度。对65、73和93的理论熔解温度做线性回归，长度73的熔解温度整体靠近线性回归性(73长度那组的熔解温度的线性相关系数最高)，因此认为73长度的整体性能更加接近。

选择长度为73bp的DNA序列集用于后续分析因为此长度有更为均衡的理论熔解温度和相对较好的耐突变性。

(2)探针序列的筛选或性能评估：通过BWA(Burrows-Wheeler Aligner)[23]将上述生成的每一条DNA序列(在该示例中也称为探针)比对回由全部线粒体组装序列构成的一个线粒体参考序列库。为使最终的探针能够特异性捕获目标区域，比对次数超过2次及以上的DNA序列将被剔除。将全部的参考序列打断成15-mer，构建15-mer的频数表，通过对参考序列的15-mer的频数表进行分析，定位频数表的中位数，以该数值上下浮动20作为频数的筛选阈值，之后对所有剩余的DNA序列计算其15-mer频数表，若该DNA序列的频数值超出筛选阈值(预定范围)，则过滤掉该DNA序列。

在进行上述筛选后，还可进一步考虑一些参数来对剩余的DNA序列进行排序用以进行下一步，这些参数包括：1)GC含量，将GC含量处于[40％,60％]区间的DNA序列排在前面，2)理论熔解温度，将满足1)的且熔解温度更接近平均熔解温度的DNA序列排在前面。

(3)获取探针集：将线粒体组装序列(参考序列)分成45bp长，间隔距离20bp的窗口。随后同一窗口内的DNA序列根据排名，将排名最高/最前的作为此窗口的探针。此外，对GC含量高于60％的探针进行一次拷贝。

(4)探针合成：将探针集导入到CustomArray B3TM合成仪中，根据标准流程合成DNA寡核苷酸后用浓氨水洗脱收集。经过PCR和反转录后，最终得到187674条RNA探针。

此设计过程中，发明人利用千种昆虫转录组(1KITE)项目中379个物种线粒体上的编码蛋白基因设计了一套线粒体捕获芯片，这些物种涵盖了昆虫纲的所有目及颚足纲、软壳纲、多足纲、桨足纲和甲壳纲的代表性物种，在昆虫系统进化树上分布范围更广并补充了蜻蜓目、襀翅目、革翅目、螳螂目和膜翅目的数据。随后发明人利用49个已知物种的混合DNA检测了这套芯片的线粒体捕获效率。

2.用于杂交的混合DNA

混合49个物种DNA的混合DNA样品来自唐敏等人[18]的研究工作，这49个物种涵盖了47个属、42个科，几乎每个物种都代表了一个科，详情见参考文献。每个物种的基因组DNA按照Ivanova等人[24]发布的流程单独提取。每个物种抽取100ng DNA混合建库，插入片段长度为200bp。探针通过生物素与磁珠连接后与DNA文库杂交72小时。最后利用IlluminaHiseq2000的100bp双末端法对捕获的DNA测序。

3.生物信息学分析

对测序后的数据进行分析，评估探针捕获率及对各个物种进行丰度和覆盖度的预估，具体步骤如下：

(1)探针遗传距离和AT含量评估，该步为可选步骤：遗传距离表示两条序列上核苷酸位点的差异性。通过ClustalW[25](断开及断开延伸的参数为50)进行全长比对，之后计算遗传距离和AT含量。图2呈现的是每个蛋白编码基因在分类阶元(大部分是目一级)里的平均遗传距离和AT含量。

(2)对49个物种的线粒体基因进行分析，该步为可选择步骤：49个线粒体基因组数据来自唐敏等人的报道[18]。大致上，现有方式对照结果的获得可以通过通过：混合DNA不经过捕获直接在Hiseq2000上进行测序，总共产生35G数据，随后利用软件TGICL[26]结合之前利用SOAPdenovo，SOAPDenov trans和IDBA_UD得到的组装结果进行分析，最后根据公共数据库中最相近的分类信息将组装结果匹配到每个物种上。除了Aposthonia borneensis(婆罗洲丝蚁)缺少CYTB及Opiliones缺少ND4，ND4L，ND6和CYTB之外[18]，49个线粒体基因组包含其他全部的蛋白编码基因，缺失的这些基因可作为后续分析的参考，线粒体蛋白分布情况见图2。通过以下步骤寻找与每个蛋白编码基因最相近的探针：1)每个编码蛋白基因与探针数据集比对寻找同源性最近的同源基因，2)利用ClustalW(断开及断开延伸的参数为50)将与每个蛋白编码基因最近似的5条探针做全长比对，3)拥有最小遗传距离的探针作为最近似的同源基因而相应的遗传距离会用于随后的分析。

(3)丰度和覆盖度评估：利用BWA将原始测序短序列与49个物种的线粒体基因组进行比对，错配两次以上的测序短序列会被过滤掉。之后通过锚定(比对上/匹配)的测序短序列(读段，reads)数计算区域丰度。考虑到线粒体片段即使没有通过捕获富集也会被随机测序，超过3个测序短序列比对上的区域是在计算覆盖度的时候算有效区域。

基于上述技术流程，能够大幅度提高DNA提取物中线粒体的捕获效率，并可以通过生物信息学分析手段进行物种丰度和覆盖度的预估。利用直接测序和芯片捕获后测序的数据计算测序短序列与49个线粒体基因组的匹配度。结果表明测序短序列与参考数据的匹配度从捕获之前的0.47％大幅度提高到了42.52％。

通过后续的生物信息分析结果表明，捕获芯片得到的每个物种的丰度和他们的初始丰度一致，说明该芯片可用于生物量的评估。发明人惊奇的发现，比对不同物种之间丰度差异时发现，遗传距离较小的物种与遗传距离较大的物种相比更倾向于拥有更高的物种丰度，这就会在某种程度上影响最终的丰度评估结果。对于那些研究特定类群或某些环境指示生物的研究来说，定制化探针更有利于获得更准确的物种丰度信息。这也进一步扩大了未来捕获芯片的应用范围。

以下结合具体实施例对依据本发明的方法和/或试剂盒获得的结果进行详细的描述。下面示例，仅用于解释本发明，而不能理解为对本发明的限制。除另有交待，以下实施例中涉及的未特别交待的试剂、序列(接头、标签和引物)、软件及仪器，都是常规市售产品或者开源的，比如购自Illumina公司的hiseq2000测序平台建库相关试剂盒来进行文库构建等。

实施例

为了保证探针的覆盖度，选取了来自1KITE项目的线粒体数据。线粒体基因组利用第一批379个物种的1KITE数据获得，之后通过基因注释得到线粒体上的蛋白编码基因。最终得到2,553条平均长度1,902bp的组装序列用于随后的探针设计。在探针设计的过程中，首先将线粒体组装序列分长度为64、73和93的序列集，综合考虑理论熔解温度和相对较好的耐突变型最终选择长度为73bp的序列集，即将1KITE的2553条组装序列分割成73bp长的探针。随后根据探针的特异性区域情况、理论熔解温度和GC含量将所有的探针进行综合排名。将所有的线粒体组装序列分割成45bp的窗口，滑动距离为20bp，在每个滑动窗口中将此窗口里所有的探针根据之前的综合排名，选取排名最前/高的一条探针作为此窗口的探针，将所有窗口的探针集合起来的到最后的探针集并合成。

验证探针所用的49个物种混合DNA样品来自唐敏等人的研究工作，这49个物种涵盖了47个属、42个科，几乎每个物种都代表了一个科，详见参考文献[18]。每个物种的基因组DNA按照Ivanova等人发布的流程单独提取。每个物种抽取100ng DNA混合建库，插入片段长度为200bp。探针通过生物素与磁珠连接后与DNA文库杂交72小时。最后利用IlluminaHiseq2000的100bp双末端法对捕获的DNA测序。

随后的生物信息学分析对测序后的数据进行分析，并评估探针捕获率及各个物种丰度和覆盖度的预估。首先评估探针遗传距离和AT含量，遗传距离表示两条序列上核苷酸位点的差异性。通过ClustalW(断开及断开延伸的参数为50)进行全长比对，之后计算遗传距离和AT含量。图2呈现的是每个蛋白编码基因在分类阶元(大部分是目一级)里的平均遗传距离和AT含量。

验证用的49个线粒体基因组数据来自唐敏等人的报道。大致上，混合DNA不经过捕获直接在Hiseq2000上进行测序，总共产生35G数据。随后利用TGICL软件结合之前利用SOAPdenovo，SOAPDenov trans和IDBA_UD得到的组装结果进行分析。最后根据公共数据库中最相近的分类信息将组装结果匹配到每个物种上。除了Aposthonia borneensis缺少CYTB及Opiliones缺少ND4，ND4L，ND6和CYTB之外，49个线粒体基因组包含其他全部的蛋白编码基因，缺失的这些基因可作为后续分析的参考。之后通过以下步骤寻找与每个蛋白编码基因最相近的探针：1)每个编码蛋白基因与探针数据集比对寻找同源性最近的同源基因，2)利用ClustalW(断开及断开延伸的参数为50)将与每个蛋白编码基因最近似的5条探针做全长比对，3)拥有最小遗传距离的探针作为最近似的同源基因而相应的遗传距离可用于随后的分析。

最后利用BWA将原始测序短序列与49个物种的线粒体基因组进行比对进行丰度和覆盖度评估，错配两次以上的测序短序列会被过滤掉。之后通过锚定的测序短序列数计算区域丰度。考虑到线粒体片段即使没有通过捕获富集也会被随机测序，超过3个测序短序列比对上的区域是在计算覆盖度的时候算有效区域，每个个体相对丰度预估值见图3。

上述结果表明不仅线粒体丰度数据利用率由最初的0.47％提高到了42.52％，还可对每个个体的相对丰度进行预估，与初始丰度的相关性较高(r²＝0.81)。

结合图2和物种丰度评估结果，对比不同物种之间丰度差异时发明人发现，遗传距离较小的物种与遗传距离较大的物种相比更倾向于拥有更高的物种丰度，这就会在某种程度上影响最终的丰度评估结果。对于那些研究特定类群或某些环境指示生物的研究来说，定制化探针更有利于获得更准确的物种丰度信息。这也进一步扩大了未来捕获芯片的应用范围。

[参考文献]

1.Board,M.A.,Millennium ecosystem assessment.Washington,DC:NewIsland,2005.

2.Naidoo,R.,et al.,Effect of biodiversity on economic benefits fromcommunal lands in Namibia.Journal of Applied Ecology,2011.48(2):p.310-316.

3.Keesing,F.,et al.,Impacts of biodiversity on the emergence andtransmission of infectious diseases.Nature,2010.468(7324):p.647-52.

4.Yu,D.W.,et al.,Biodiversity Soup:Metabarcoding of arthropods forrapid biodiversity assessment and biomonitoring.Methods in Ecology andEvolution,2012.3(4):p.613-623.

5.Hajibabaei,M.,et al.,Environmental barcoding:a next-generationsequencing approach for biomonitoring applications using river benthos.PloSone,2011.6(4):p.e17497.

6.Zhou,X.,et al.,Ultra-deep sequencing enables high-fidelity recoveryof biodiversity for bulk arthropod samples without PCRamplification.GigaScience,2013.2(1):p.4.

7.Liu,S.,et al.,SOAPBarcode:revealing arthropod biodiversity throughassembly of Illumina shotgun sequences of PCR amplicons.Methods in Ecologyand Evolution,2013.4(12):p.1142-1150.

8.Porazinska,D.L.,et al.,Evaluating high throughput sequencing as amethod for metagenomic analysis of nematode diversity.Molecular EcologyResources,2009.9(6):p.1439-1450.

9.Porazinska,D.L.,et al.,Reproducibility of read numbers in high-throughput sequencing analysis of nematode community composition andstructure.Molecular Ecology Resources,2010.10(4):p.666-676.

10.Baird,D.J.and M.Hajibabaei,Biomonitoring 2.0:a new paradigm inecosystem assessment made possible by next-generation DNAsequencing..Molecular Ecology,2012.21(8):p.2039-2044.

11.Bienert,F.,et al.,Tracking earthworm communities from soilDNA.Molecular Ecology,2012.21(8):p.2017-2030.

12.Coissac,E.,T.Riaz,and N.Puillandre,Bioinformatic challenges forDNA metabarcoding of plants and animals.Molecular Ecology,2012.21:p.1834–1847.

13.Taberlet,P.,et al.,Towards next‐generation biodiversity assessmentusing DNA metabarcoding.Molecular Ecology,2012.21(8):p.2045-2050.

14.Hahn,C.,L.Bachmann,and B.Chevreux,Reconstructing mitochondrialgenomes directly from genomic next-generation sequencing reads—a baiting anditerative mapping approach.Nucleic Acids Research,2013:p.gkt371.

15.Williams,S.,P.Foster,and D.Littlewood,The complete mitochondrialgenome of a turbinid vetigastropod from MiSeq Illumina sequencing of genomicDNA and steps towards a resolved gastropod phylogeny.Gene,2014.533(1):p.38-47.

16.Groenenberg,D.S.,et al.,The complete mitogenome of Cylindrusobtusus(Helicidae,Ariantinae)using Illumina next generation sequencing.BMCgenomics,2012.13(1):p.114.

17.Timmermans,M.J.,et al.,Why barcode？High-throughput multiplexsequencing of mitochondrial genomes for molecular systematics.Nucleic acidsresearch,2010.38(21):p.e197-e197.

18.Min,T.,et al.,Multiplex sequencing of pooled mitochondrialgenomes—a crucial step toward biodiversity analysis using mito-metagenomics.Nucleic Acids Research,2014:p gku917.

19.Bamshad,M.J.,et al.,Exome sequencing as a tool for Mendeliandisease gene discovery.Nature Reviews Genetics,2011.12(11):p.745-755.

20.Vallender,E.J.,Expanding whole exome resequencing into non-humanprimates.Genome biology,2011.12(9):p.R87.

21.Lemmon,A.R.,S.A.Emme,and E.M.Lemmon,Anchored hybrid enrichment formassively high-throughput phylogenomics.Systematic biology,2012:p.sys049.

22.Li,C.,et al.,Capturing protein-coding genes across highlydivergent species.BioTechniques,2013.54(6):p.321-326.

23.Li,H.and R.Durbin,Fast and accurate short read alignment withBurrows–Wheeler transform.Bioinformatics,2009.25(14):p.1754-1760.

24.Ivanova,N.V.,J.R.Dewaard,and P.D.Hebert,An inexpensive,automation‐friendly protocol for recovering high‐quality DNA.Molecular ecology notes,2006.6(4):p.998-1002.

25.Larkin,M.A.,et al.,Clustal W and Clustal X version2.0.Bioinformatics,2007.23(21):p.2947-2948.

26.Pertea,G.,et al.,TIGR Gene Indices clustering tools(TGICL):asoftware system for fast clustering of large EST datasets.Bioinformatics,2003.19(5):p.651-652.

Claims (10)

1.一种确定探针集的方法，所述探针集包含多个探针，其特征在于，所述方法包括，

获取DNA序列集，所述DNA序列集包含多条长度为K₀的DNA序列，获取所述DNA序列集包括，

从长度为L的参考序列的一端的第i个核苷酸开始，沿另一端方向拷贝所述参考序列的K₀个连续核苷酸为一条DNA序列，i依次取{1，2，…，(L-K₀)，(L-K₀+1)}中的数值；

筛选所述DNA序列集，以获得所述探针，其中包括，过滤掉与所述参考序列的匹配位置个数大于1的DNA序列；其中，

L＞K₀，200nt≥K₀≥50nt；

其中筛选所述DNA序列集还包括，计算过滤后的各条DNA序列的长度为K₁的单元序列的频数，过滤掉频数超出预定范围的DNA序列，

所述预定范围是依据参考序列的单元序列的频数的中位数来确定的，所述参考序列的单元序列的频数中位数通过以下确定，

将所述参考序列拆分成多个长度为K₁的单元序列，所述单元序列通过以下获得，

从长度为L的参考序列的一端的第j个核苷酸开始，沿另一端方向拷贝所述参考序列的K₁个连续核苷酸为一条单元序列，j依次取{1，2，…，(L-K₁)，(L-K₁+1)}中的数值，

依据序列相同的单元序列分到相同的组，对所述单元序列进行分组，获得q个单元序列组，将单元序列组r包含的单元序列的个数表示为A_r，A_r为所述参考序列对单元序列组r中的单元序列的频数，

取{A₁，A₂，…，A_q}中的中位数为所述参考序列的单元序列的频数中位数，其中，

4*K₁≤K₀，r为单元序列组的编号，1≤r≤q，q≤L-K₁+1，A_r≥1。

2.权利要求1的方法，其特征在于，所述参考序列为线粒体参考序列。

3.权利要求1的方法，其特征在于，所述参考序列为多个物种的线粒体参考序列，所述物种为昆虫。

4.权利要求1的方法，其特征在于，所述DNA序列的单元序列的频数＝匹配到该DNA序列的单元序列的数目/包含匹配到该DNA序列的单元序列的单元序列组的个数。

5.权利要求1的方法，其特征在于，所述预定范围为0.8-1.2倍的所述参考序列的单元序列的频数中位数。

6.权利要求1的方法，其特征在于，筛选所述DNA序列集还包括，

将所述参考序列划分成不连接的大小为m的窗口，相邻窗口间距为n，

将过滤后的DNA序列与所述参考序列比对，获得各个窗口的比对结果，

当一个窗口的比对结果包含的DNA序列数目大于等于1且包含的DNA序列中至少有一条的GC含量为40％～60％时，选择其中的GC含量为40％～60％的一条DNA序列，剔除该窗口比对结果包含的其它DNA序列，

当一个窗口的比对结果包含的DNA序列数目大于等于1且不包含GC含量为40％～60％的DNA序列时，选择其中的GC含量最接近60％的一条DNA序列并且拷贝该条DNA序列一次，剔除该窗口比对结果包含的其它DNA序列，其中，K₀+10≥m+n≥K₀，m≥n。

7.权利要求6的方法，其特征在于，还包括，

反转录筛选后的DNA序列，获得RNA序列。

8.一种试剂盒，其包括探针集，所述探针集能够特异性识别线粒体DNA，所述探针集利用权利要求1-7任一方法确定。

9.权利要求8的试剂盒在获取线粒体DNA中的用途。

10.一种确定物种丰度的方法，其特征在于，包括，

(1)获得多个物种的混合线粒体DNA的测序数据，其中包括，利用权利要求8的试剂盒捕获多个物种混合核酸中的混合线粒体DNA，所述测序数据包括多个读段；

(2)将所述测序数据与各个所述物种的线粒体参考序列比对，获得各个物种的比对结果；

(3)基于所述物种的比对结果中的读段数目，确定该物种的丰度。