本发明涉及生物学领域，特别地，涉及群体遗传学领域，更特别地，本发明涉及一种确定群体中受到选择作用的位点的方法和一种确定群体中受到选择作用的位点的装置。
背景技术：
：随着二代测序(nextgenerationsequencing，NGS)技术的成熟和成本的逐步降低，各项以此为基础、用于不同目的的研究技术层出不穷。RNA-Seq是一种基于NGS，通过对样品的转录组(transcriptome)进行测序，主要用于揭示样本中基因表达规律的技术，现已被广泛运用。同时，RNA-Seq的测序数据也可用于检测整个基因组转录区域的多态性位点，包括SNP位点。技术实现要素：依据本发明的一方面，本发明提供一种确定群体中受到选择作用的位点的方法，所述选择作用包括人工选择作用和自然选择作用的至少一种，该方法包括以下步骤：(1)获得群体样本的核酸测序数据，所述群体样本来自一个物种的多个个体，任选的，所述群体样本来自一个物种多个个体的相同组织或者一个物种的多个个体的相同部位，所述群体样本能够依据n对预定指标划分成2n个一级亚群体，n为自然数；(2)基于(1)中的核酸测序数据，检测以获得群体SNP数据，所述群体SNP数据包括多个一级亚群体SNP数据；(3)基于(2)中的群体SNP数据，比较不同一级亚群体的多态性的差异，以确定受到选择作用的SNP，所述受到选择作用的SNP为所述受到选择作用的位点。在本发明的一个实施例中，所述核酸测序数据是利用RNA-Seq技术得来的，为转录本测序数据。所称的预定指标可以是任意的两个个体样本的不一样的特征，在本发明的一个实施例中，预定指标是地理的和/或生物性状相关的，例如可以以不同地域来源、具有某个(些)不同性状等来作为初步划分群体的指标。在本发明的一个实施例中，在进行该方法的步骤(3)之前或者步骤(3)之后，进行群体结构分析，包括：基于(2)中的群体SNP数据，对所述群体样本进行群体结构分析，获得群体结构分析结果；任选的，进行所述群体结构分析包括构建系统发育树、主成分分析和STRUCTURE分析中的至少之一。而且，在本发明的另一个实施例中，进一步的，基于所述群体结构分析结果，对所述群体样本进行重新划分，以获得的划分结果即对群体的分类结果替代原先所述 一级亚群体，进而进行(3)来确定群体中受到选择作用的位点。依据本发明的另一方面，本发明提供一种基于群体转录本数据分析群体结构的方法，该方法包括：获得群体样本的核酸测序数据，所述群体样本来自一个物种的多个个体，任选的，所述群体样本来自一个物种多个个体的相同组织或者一个物种的多个个体的相同部位，所述群体样本能够依据n对预定指标分成2n个一级亚群体，n为自然数；基于所述核酸测序数据，检测以获得群体SNP数据，所述群体SNP数据包括多个一级亚群体SNP数据；基于所述群体SNP数据，比较不同一级亚群体的多态性的差异，确定受到选择作用的SNP，和/或，基于所述群体SNP数据，对所述群体进行群体结构分析。依据本发明的再一方面，本发明提供一种确定群体中受到选择作用的位点的装置，该装置用以实施上述本发明一方面的确定群体中受到选择作用的位点的方法，装置包括：数据输入单元，用于输入数据；数据输出单元，用于输出数据；处理器，用于执行机器可执行程序，执行所述机器可执行程序包括完成本发明一方面的或者任一实施例中的方法；存储单元，与所述数据输入单元、数据输出单元和处理器相连，用于存储数据，其中包括所述机器可执行程序。本领域技术人员能够理解，所说的机器可执行程序可以保存在存储介质中，所称存储介质可以包括：只读存储器、随机存储器、磁盘或光盘等。依据本发明的又一方面，本发明提供一种确定群体中受到选择作用的位点的系统，该系统能够用以实施上述本发明一方面的或者任一实施例中的方法的全部或部分步骤，该系统包括：测序数据获取装置，用以获取群体样本的核酸测序数据，所述群体样本来自一个物种的多个个体，任选的，所述群体样本来自一个物种多个个体的相同组织或者一个物种的多个个体的相同部位，所述群体样本能够依据n对预定指标分成2n个一级亚群体，n为自然数；SNP检测装置，与所述测序数据获取装置连接，用于基于所述核酸测序数据，检测以获得群体SNP数据，所述群体SNP数据包括多个一级亚群体SNP数据；目的位点确定装置，与所述SNP检测装置连接，用以基于所述群体SNP数据，比较不同一级亚群体的多态性的差异，以确定受到选择作用的SNP，所述受到选择作用的SNP为所述受到选择作用的位点。利用上述本发明的方法、装置和/或系统能够准确的判定出群体中受到选择作用的位点。本发明的方法和/或装置，集中于更具普遍重要性的基因组转录区域，能够基于获得的群体转录本数据，获得基因表达数据，揭示样本的基因表达规律，这将有利于揭示遗传背景差异条件下的基因表达规律，是对RAD、GBS等群体研究范围的进一步拓展。而且，又能够获得群体SNP数据，揭示群体结构和群体遗传进化规律。本发明方法、装置和/或系统能够用以规范群体转录组重测序分析流程，降低分析风险，能够高效率、高质量和高标准完成对群体项目的分析。附图说明本发明的上述和/或附加的方面和优点从结合下面附图对实施方式的描述中将变得明显和容易理解，其中：图1是本发明的一个实施例中的确定群体中受到选择作用的位点的方法的步骤流程图。图2是本发明的一个实施例中的确定群体中受到选择作用的位点的方法的步骤流程图。图3是本发明的一个实施例中的确定群体中受到选择作用的位点的方法的步骤流程图。图4是本发明的一个实施例中的确定群体中受到选择作用的位点的装置示意图。图5是本发明的一个实施例中的确定群体中受到选择作用的位点的系统示意图。图6是本发明的一个实施例中的Frappe基于群体SNP推测的群体遗传结构的示意图。图7是本发明的一个实施例中的基于群体SNPs采用邻接法推断的系统发生树的示意图。图8是本发明的一个实施例中的基于群体SNP的PCA分析结果示意图。图9是本发明的一个实施例中的Arlequin程序基于群体SNP检测受选择作用位点的结果示意图。图10是本发明的一个实施例中的GlobalFSTtest程序基于群体SNP检测受选择作用位点的结果示意图。图11是本发明的一个实施例中的BayeScan程序基于群体SNP检测受选择作用位点的结果示意图。具体实施方式下面详细描述本发明的实施例，所述实施例的示例在附图中示出，其中，自始至终相同或类似的标号表示相同或类似的元件或具有相同或类似功能的元件。下面通过参考附图描述的实施例是示例性的，仅用于解释本发明，而不能理解为对本发明的限制。需要说明的，本文中所使用的术语“一级”、“二级”等仅为方便描述，不能理解为指示或暗示相对重要性，也不能理解为之间有先后顺序关系。在本发明的描述中，除非另有说明，“多个”的含义是两个或两个以上。在本文中，除非另有明确的规定和限定，术语“相连”、“连接”等术语应做广义理解，例如，可以是固定连接，也可以是可拆卸连接，或一体地连接；可以是机械连接，也可以是电连接；可以是直接相连，也可以通过中间媒介间接相连，可以是两个元件内部的连通。根据本发明的一个实施例，如图1所示，本发明提供一种确定群体中受到选择作用的位点的方法，所述选择作用包括人工选择作用和自然选择作用的至少一种，该方法包括以下步骤：S10获得群体样本的核酸测序数据，所述群体样本来自一个物种的多个个体，任选的， 所述群体样本来自一个物种多个个体的相同组织或者一个物种的多个个体的相同部位，所述群体样本能够依据n对预定指标划分成2n个一级亚群体，n为自然数；S20基于S10中的核酸测序数据，检测以获得群体SNP数据，所述群体SNP数据包括多个一级亚群体SNP数据；S30基于S20中的群体SNP数据，比较不同一级亚群体的多态性的差异，以确定受到选择作用的SNP，所述受到选择作用的SNP为所述受到选择作用的位点。根据本发明的一个实施例，所述核酸测序数据是利用RNA-Seq技术得来的，为转录本测序数据。以同一物种、多个不同遗传背景的个体为研究对象，通过对转录组(transcriptome)样品进行高通量测序，一次性获得该特定物种群体水平的基因组转录区域多态性数据，包括群体SNP数据和全基因/转录本表达信息，可以用于揭示研究个体之间的进化关系和遗传组成差异、在特定选择作用下共同进化的基因簇、亚群体中受人工/自然选择作用的位点以及个体或亚群体之间的在表达上具有显著差异的功能模块和代谢通路等生物学问题。而且，相对于常规的少量样品的转录组重测序，相比于RAD、GBS等群体研究技术，本发明的研究区域相对集中于基因组转录区域，可以对基因表达进行定量，这将有利于揭示遗传背景差异条件下的基因表达规律，是对RAD、GBS等群体研究范围的进一步拓展。所称的预定指标可以是任意的两个个体样本的不一样的特征，根据本发明的一个实施例，预定指标是地理的和/或生物性状相关的，例如可以以不同地域来源、具有某个(些)不同性状等来作为初步划分群体的指标。根据本发明的一个实施例，如图2所示，在进行该方法的步骤S30之前，还包括进行S23群体结构分析，S23群体结构分析包括：基于S20中的群体SNP数据，对所述群体样本进行群体结构分析，获得群体结构分析结果；任选的，进行所述群体结构分析包括构建系统发育树、主成分分析(PCA)和GroupStructure分析中的至少之一。可以利用邻接法构建系统发育树，也可以利用MEGA软件构建关系，利用MEGA软件(http://www.megasoftware.net)，将每个样本所有SNP位点的基因型文件组成序列，一个个体样本对应一个序列，作为MEGA的输入文件，MEGA根据各个体样本序列上的差异，该软件有三种方法(Maximumlikelihood、LeastSquares和Maximumparsimony)来构建关系树。在统计学中，主成分分析(PrincipalComponentsAnalysis，PCA)是一种简化数据集的技术，是一个线性变换。这个变换把数据变换到一个新的坐标系统中，使得任何数据投影的第一大变量数在第一个坐标(称为第一主成分)上、第二大变量数在第二坐标(第二主成分)上，依次类推。主成分分析经常用于减少数据集的维数，同时保留对数据集贡献最大的特征变量。通过保留低阶主成分，忽略高阶主成分来实现的。这是由于低阶成分往往能够保留数 据集中最重要的方面。根据参考文献AtutorialonPrincipalComponentsAnalysis.LindsayISmith,2002-02和实施例中真实的SNP数据特点，首先将SNP数据转换成数字矩阵，例如设定与参考序列一致的为0、相反的为2、简并碱基为1，并作均一化。然后通过上述介绍的方法构建线性向量方程。其中i从1到k表示第i个样本。应用R语言软件包强大的解方程能力，解得矩阵a，根据各样本的数据特点提取前四个主成分向量，以向量作为坐标轴展示各个体聚类情况。GroupStructure分析可以利用Structure软件(http://pritch.bsd.uchicago.edu/software/structure2_1.html)进行，该软件基于SNP位点的基因分型数据，推断是否存在不同群体并判断每个个体所归属的群体。根据软件说明，将群体SNP的基因型文件转换格式，作为Structure输入文件并在混合模型中采用高达5万次模拟，假设多个群体存在时，计算每个个体归属各类(亚)群体的概率。经过以上，能够实现对个体的分类。在本发明的一个实施例中，在分类的基础上，还可以进一步筛选个体，例如根据上述群体结构分析结果，实现对个体的分类，提取每个个体样本信息，剔除存在异议的个体，比如分类不明确或明显离群样本。根据本发明的一个实施例，进一步的，基于所述群体结构分析结果，对所述群体样本进行重新划分，以获得的划分结果即获得的新的亚群体替代原先的一级亚群体，进而基于新的亚群体及其SNP数据进行步骤S30来确定群体中受到选择作用的位点，这样，以群体结构分析结果对群体/亚群体进行再分类或者重新分类，有利于准确判定受到选择作用的位点。根据本发明的一个实施例，如图3所示，在进行该方法的步骤S30之后，还包括进行S23群体结构分析，S23群体结构分析包括：基于S20中的群体SNP数据，对所述群体样本进行群体结构分析，获得群体结构分析结果；任选的，进行所述群体结构分析包括构建系统发育树、主成分分析(PCA)、GroupStructure分析和种群遗传结构Frappe检测中的至少之一。根据本发明的一个实施例，所述群体样本的核酸测序数据由组成群体样本的每个个体样本的核酸测序数据组成，要求每个个体样本的核酸测序数据不少于4G，以利于准确检测出SNP，进而有利于基于准确的群体SNP数据准确确定受到选择作用位点。根据本发明的一个实施例，群体样本来自同一物种的、具有不同遗传背景的个体。对于群体样本分析，建议群体样本中包含的个体样本数量不小于30个，而且，涉及的所有个体至少能够根据某种指标而被划分为两个及两个以上的亚群体，即所称的一级亚群体，以便于后续差异分析。根据本发明的一个实施例，较佳的，每个一级亚群体包括至少10个个体样本，以利于差异分析。根据本发明的一个实施例，将所有个体样本在相同的条件下进行培养，然 后在相同的组织或者部位取样，来获得群体样本，这样使得基于该群体样本数据进行群体分析包括进行基因差异表达分析有意义，原因在于，个体样本的遗传差异即变量已经存在，在相同条件下取样，能够使得到的差异表达基因能从遗传差异的角度去作解释，否则，多个变量的存在，会导致差异表达的原因模棱两可。例如，研究群体可以被分为抗盐碱和不抗盐碱两类，可以使用相同计量的盐水对生长在相同环境下的所有个体进行处理，然后对处理后特定时间(例如1小时)的根尖进行取样，这样，后续群体分析鉴定出来的差异表达基因可能能用于揭示此物种抗盐碱的机制，并且，能确定该差异表达是由于遗传背景的差异导致的。根据本发明的一个实施例，所述一级亚群体包括至少一个二级亚群体；任选的，一个所述二级亚群体包括至少10个个体。二级亚群体可以通过利用不同于划分群体的另一个(些)指标划分一级亚群体来获得。利用本发明的任一实施例中的方法能够对多次划分后的多级亚群体中的受到选择作用的位点进行准确判定。根据本发明的一个实施例，所述基于群体SNP数据，比较不同一级亚群体多态性的差异，以确定受到选择作用的SNP，包括：基于群体SNP数据，利用至少两种检验方法比较所述不同一级亚群体中的相同SNP位点的杂合度的差异，将得到至少两种检验方法支持的SNP位点确定为受到选择作用的SNP；任选的，所述检验方法包括F统计量、分子变异分析和多层贝叶斯方法。在本发明的一些实施例中，利用Arlequin程序、GlobalFSTtest程序和BayeScan程序中的两个或者全部三个，或者包括利用Arlequin，BayesScan和Datacal三种方法中的至少两个或者全部三种方法来判断比较位点的杂合度差异程度，当某SNP位点得到以上三种检验方法中的至少两种或者全部三种的支持，即至少其中的两种的检验结果都认定该SNP在不同亚群体中的杂合度的差异是显著的，则判定该SNP为受到选择作用的位点。这样，有利于准确判定。根据本发明的一个实施例，所述利用至少两种检验方法来比较所述不同一级亚群体中的相同SNP位点的杂合度的差异，将得到至少两种检验方法支持的SNP位点确定为受到选择作用的SNP，包括：计算所述SNP位点在不同一级亚群体中的杂合度差异值，将杂合度差异值不小于阈值的SNP位点确定为受到选择作用的位点。在本发明的一个实施例中，所称的杂合度差异值以FST(Fixationindex)表示。FST可以用来评价群体间的基因组距离和种群的差异，是度量种群间分化程度的一个指标，由SewallWright在1922年应用F-检验的一种特殊情况发展而来。FST的零假设是在群体没有分化时，多态性位点在(亚)群内和(亚)群间的次等位碱基的频率差别不具显著性。计算FST的方法很多，虽然具体计算方法不同，但基本理 论是一致的，即由Hudson(1992)给出的定义：其中，ΠBetween在这里表示从两个亚群体(Between)中分别抽取一个样本，组成一对，计算这对样本SNP基因型的差异，如此可以计算所有成对样本SNP基因型的差异，最后求平均值即为ΠBetween。ΠWithin表示从一个亚群体(Within)中分别抽取2个样本，组成一对，计算这对样本SNP基因型的差异，如此可以计算所有成对样本SNP基因型的差异，最后求平均值即为ΠWithin。如果有两个亚群体，可以两个亚群体分别先计算ΠWithin，然后累加。在该实施例中，结合已有亚群体SNP数据的结构，基于上述原理，推导公式如下：FST=ΠBetween-ΠWithinΠBetween=1-ΠWithinΠBetween=1-[Σj(2nj)Σj2nijnij-1(1-xij)]/Σj(2nj)Σj2nini-1xi(1-xi),]]>其中，xij是SNP位点i在亚群体j中的次等位碱基(第二碱基)的频率，而nij是SNP位点i在亚群体j中染色体上的物理位置，nj则是亚群体j中用于比较分析的SNP位点个数的总和。在本发明的一个实施例中，利用Arlequin，BayesScan和Datacal三种方法来比较检验SNP位点的次等位碱基频率在亚群体间的差异，各自设置的差异具有显著性的阈值分别为0.05，0.1和0.01。根据本发明的一个实施例，本发明提供一种基于群体转录本数据分析群体结构的方法，该方法包括：获得群体样本的核酸测序数据，所述群体样本来自一个物种的多个个体，任选的，所述群体样本来自一个物种多个个体的相同组织或者一个物种的多个个体的相同部位，所述群体样本能够依据n对预定指标分成2n个一级亚群体，n为自然数；基于所述核酸测序数据，检测以获得群体SNP数据，所述群体SNP数据包括多个一级亚群体SNP数据；基于所述群体SNP数据，比较不同一级亚群体的多态性的差异，确定受到选择作用的SNP，和/或，基于所述群体SNP数据，对所述群体进行群体结构分析。根据本发明的一个实施例，如图4所示，本发明提供一种确定群体中受到选择作用的位点的装置100，该装置100用以实施上述本发明一方面的确定群体中受到选择作用的位点的方法，装置100包括：数据输入单元110，用于输入数据；数据输出单元120，用于输出数据；处理器130，用于执行机器可执行程序，执行所述机器可执行程序包括完成本发明一方面的或者任一实施例中的方法；存储单元140，与所述数据输入单元110、数据输出单元120和处理器130相连，用于存储数据，其中包括所述机器可执行程序。本领域技术人员能够理解，所说的机器可执行程序可以保存在存储介质中，所称存储介质可以包括：只读存储器、随机 存储器、磁盘或光盘等。根据本发明的一个实施例，如图5所示，本发明提供一种确定群体中受到选择作用的位点的系统1000，该系统能够用以实施上述本发明一方面的或者任一实施例中的方法的全部或部分步骤，该系统1000包括：测序数据获取装置1100，用以获取群体样本的核酸测序数据，所述群体样本来自一个物种的多个个体，任选的，所述群体样本来自一个物种多个个体的相同组织或者一个物种的多个个体的相同部位，所述群体样本能够依据n对预定指标分成2n个一级亚群体，n为自然数；SNP检测装置1200，与所述测序数据获取装置1100连接，用于基于所述核酸测序数据，检测以获得群体SNP数据，所述群体SNP数据包括多个一级亚群体SNP数据；目的位点确定装置1300，与所述SNP检测装置1200连接，用以基于所述群体SNP数据，比较不同一级亚群体的多态性的差异，以确定受到选择作用的SNP，所述受到选择作用的SNP为所述受到选择作用的位点。利用上述本发明任一实施例中的方法、装置和/或系统能够准确的判定出群体中受到选择作用的位点。本发明的方法和/或装置，主要集中于更具普遍重要性的基因组转录区域，能够基于获得的群体转录本数据，获得基因表达数据，揭示样本的基因表达规律，这将有利于揭示遗传背景差异条件下的基因表达规律，是对RAD、GBS等群体研究范围的进一步拓展。而且，又能够获得群体SNP数据，揭示群体结构和群体遗传进化规律。本发明方法、装置和/或系统能够用以规范群体转录组重测序分析流程，降低分析风险，能够高效率、高质量和高标准完成对群体项目的分析。以下结合附图和具体样本数据实施例对本发明的确定受到选择作用的位点的方法、群体项目分析装置和/或系统进行详细的说明。通过参考附图描述的实施例是示例性的，仅用于解释本发明，而不能理解为对本发明的限制。除另有交待，以下实施例中涉及的未特别交待的试剂、序列(接头、标签和引物)、软件及仪器，都是常规市售产品或者开源的，例如购买Illumina的转录组文库构建试剂盒。实施例一参考序列、测序策略、样品要求及其他注意事项：i)参考序列：要求有用较高质量的基因组参考序列。ii)测序策略：采用PE91(双末端测序，获得多对paired-endreads，每条reads的长度都为91bp)的测序策略，单个样品达到过滤后数据量4G的标准。iii)样品应该来自同一物种、具有不同遗传背景的个体。iv)对于总的研究群体，推荐30个个体及以上的体量。同时，涉及的所有个体能够根据某种指标而被划分为两个及两个以上的亚群体(便于差异分析)，且每个亚群体最好多于10个个体。v)将所有样本在相同的条件下进行培养，然后在相同的组织、部位取样。原因在于，样品的遗传差异(变量)已经存在，只有在相同条件下取样，得到的差异表达基因才可能从遗传差异的角度去作解释。否则，多个变量的存在，会导致差异表达的原因模棱两可。例如研究群体可以被分为抗盐碱和不抗盐碱两类。可以使用相同计量的盐水对生长在相同环境下的所有个体进行处理，然后对处理后特定时间(例如1小时)的根尖进行取样。那么后续鉴定出来的差异表达基因则可能揭示出此物种抗盐碱的机制，因为，差异表达是由于遗传背景的差异导致的。为规范群体转录组重测序项目的分析流程，降低分析风险，以达到高效率、高质量、高标准完成项目的目的，在此提出一种群体转录组重测序分析方法，主要包括：一、实验流程提取样本总RNA并使用DNaseI消化DNA后，用带有Oligo(dT)的磁珠富集真核生物mRNA(若为原核生物，则用试剂盒去除rRNA后进入下一步骤)；加入打断试剂在Thermomixer中适温将mRNA打断成短片段，以打断后的mRNA为模板合成一链cDNA，然后配制二链合成反应体系合成二链cDNA，并使用试剂盒纯化回收粘性末端修复、cDNA的3'末端加上碱基"A"并连接接头，然后进行片段大小选择，最后进行PCR扩增；构建好的文库用Agilent2100Bioanalyzer和ABIStepOnePlusReal-TimePCRSystem质检合格后，使用IlluminaHiSeqTM2000或其他测序仪进行测序。二、信息分析内容1)标准RNA-Seq分析包括数据过滤，基因表达定量，组间差异基因鉴定及其GO、KEGGPathway富集分析，SNPcalling及注释等。2)基于群体SNP数据的分析基于标准RNA-Seq分析中对单个样品的一致性序列(consensus序列)的预测，即SNP识别(SNPcalling)的中间步骤，整理得到群体水平的SNP数据，用于下述多个方面的分析：a、群体结构分析：包括构建系统发育树、主成分(PCA)分析和STRUCTURE分析，三者都能够反映出群体的结构，但每个分析侧重点又有所不同。构建系统发育树侧重于揭示群体中个体之间的进化关系；主成分(PCA)分析侧重于揭示群体中个体之间遗传背景差异的主要因素；STRUCTURE分析侧重于对每个个体的遗传组成进行比较、量化，并以图示的 方式揭示个体之间遗传组成的异同。b、检测受到选择作用的位点：选择作用(来自于人工or自然)通常在种群的分化(亚群的形成)过程中起着非常重要的作用。从亚群的SNP数据出发，可以统计出所有位点在不同亚群之间多态性的差异(Fst)，并检验出Fst显著差异的位点。这些位点作为潜在的受到选择作用的位点，能够辅助研究者进一步认识针对于某些亚群的选择作用的过程。Fst(Fixationindex)主要用来评价群体间的基因组距离和种群的差异，是度量种群间分化程度的一个指标，由SewallWright在1922年应用F-检验的一种特殊情况发展而来。FST的零假设是在群体没有分化时，多态性位点在群内和群间的次等位碱基的频率差别是不显著的。计算FST的方法很多，虽然具体计算方法不同，但基本理论是一致的，即由Hudson(1992)给出的定义：FST=ΠBetween-ΠWithinΠBetween,]]>其中ΠBetween表示从两个群体(Between)中分别抽取一个样本，组成一对，计算这对样本SNP基因型的差异，如此可以计算所有成对样本SNP基因型的差异，最后求平均值即为ΠBetween。ΠWithin表示从一个群体(Within)中分别抽取2个样本，组成一对，计算这对样本SNP基因型的差异，如此可以计算所有成对样本SNP基因型的差异，最后求平均值即为ΠWithin。如果有两个群体，可以两个群体分别先计算ΠWithin，然后累加。3)基于基因表达数据的额外分析a、聚类分析、PCA分析：基于基因表达数据，可以对群体中的个体进行聚类、PCA分析，呈现个体与个体之间在基因表达层次上的差异。这一结果可与SNP数据构建出来的系统发育树和PCA分析结果相互印证、比较。b、共表达基因网络构建和组间比较：在各种生命活动中，多个基因(co-expressiongenes)通常在很多条件下协同地表达，以实现某些特定的功能。从多个不同个体的基因表达数据出发，可以构建出许多共表达基因的模块。以此为基础，研究者可以分析：i)在特定条件下，哪些共表达基因模块在发挥着作用(较高水平地表达)，这有利于认识这些特定条件背后的基因表达规律；ii)哪些共表达基因模块在哪个(哪些)特定的个体中发挥作用，这有利于解析部分共表达基因模块的生物学功能；iii)以上构建出的共表达基因模块还可以亚群体之间进行比较。从共表达基因模块这一更高的水平上去比较个体之间的差异，可以揭示出从常规 的基因差异表达数据(假定基因与基因之间相互独立，不考虑它们之间的相互作用)中无法体现出来的新内容。以上，以同一物种、多个不同遗传背景的个体为研究对象，通过对转录组(transcriptome)样品进行高通量测序，一次性获得该特定物种群体水平的基因组转录区域多态性数据(群体SNP)和全基因/转录本表达信息，进而可以揭示(i)研究个体之间的进化关系和遗传组成差异，(ii)在特定选择作用下共同进化的基因簇，(iii)亚群体中受人工/自然选择作用的位点，以及(iv)个体或亚群体之间的在表达上具有显著差异的功能模块和代谢通路等生物学问题。相对于常规的少量样品的转录组重测序，该方法还将给出群体SNP数据，该数据可用于揭示群体结构、群体进化历史、群体中每个个体的进化关系，以及潜在的受选择作用的位点等生物学问题。相比于RAD、GBS等群体研究技术，该方法的研究区域集中于更具普遍重要性的基因组转录区域。同时，本发明可以对基因表达进行定量，这将有利于揭示遗传背景差异条件下的基因表达规律，是对RAD、GBS等群体研究范围的进一步拓展。实施例二下面详细示例介绍分步骤操作过程：一、常规转录组重测序流程不同地域包括秦岭、岷山、梁山、邛崃和相岭的大熊猫，获取的大熊猫血液或组织样本数目总共34个，其中，来自梁山为2个——样本编号为GP37和GP52(均为血液样本)，来自岷山的有7个——样本编号为GP14-19和GP51(均为血液样本)，来自秦岭的有8个——样本编号分别为GP3-8(血液样本)、GP10(组织样本)和GP12(血液样本)，来自邛崃的有15个——样本编号分别为GP2、GP13、GP22-31、GP33和GP35-36(均为血液样本)，来自相岭的有2个——样本编号分别为GP38-39(均为血液样本)。样本转录组核酸提取、文库构建以及测序参照前面实施例进行，获得各样本测序数据。根据地域的不同，将34个样本分为5个一级亚群体。完成数据过滤、质控，将干净测序数据(cleandata)比对到基因组参考序列，比如利用SOAP或者BWA、按照其默认设置进行比对，对每个样品进行SNP识别(callsnp)，将cleandata比对到基因集参考序列上，计算每个基因的表达量并进行组间差异表达基因鉴定和GO、KEGGpathway富集分析。再次将cleandata比对到基因组参考序列，例如利用TopHat或者STAR进行比对，预测可变剪切及新的转录本，以及完成各种统计工作，包括原始、过滤后数据量统计、readsmapping信息统计、基因组覆盖度统计、生成文库随机性评估图等。二、识别(Call)群体SNP、以及基于群体SNP的群体进化分析从上一步获得的每个个体相对于基因组参考序列的consensus信息(即SOAPsnp输出的cns文件)出发，整合形成群体SNP数据，此为所有个体水平，即为取所有个体样本SNP的并集为群体SNP数据。以此群体SNP为基础，进行群体进化分析，群体进化分析包括进化树的构建、主成分分析、个体遗传组成分析等。此流程需要准备一些简单的配置文件，说明如下：individual.txt：样品(个体样本)信息文件，每一行是一个样品的信息，每行6列，如表1所示。表1snp.lst：群体SNP(genotype)文件列表，群体SNP文件格式如表2所示。表2第一列染色体编号第二列等位基因位置第三列对应参考序列位点的核苷酸第四列测序样本基因型，以空格隔开，顺序需与individual文件对应population.txt：进行位点选择分析的两个群体信息，第一列是亚群名称，可以与individual文件不同，第二列是样品缩写ID，需存在于individual文件第四列中。*.gff：基因组gff文件，进行位点选择分析时确定受选择位点所在基因，可以不提供。1)Call群体SNP利用SOAPsnp检测每个样本的SNP，整合所有单个样品的SNP数据获得群体SNP数据。具体包括：我们首先充分考虑并利用已公开的熊猫基因组信息(ZhaoS,etal.Whole-genomesequencingofgiantpandasprovidesinsightsintodemographichistoryandlocaladaptation.NatGenet.45(1):67-71(2013))，从NCBI网站下载熊猫基因组对应的dbsnp，作为SOAPsnp的先验概率，并依据目前确定的研究结果，设置杂合位点SNP的先验概率为0.0010，纯合位点SNP的先验概率为0.0005。在设置以上参数后，利用SOAPsnp软件将过滤后数据与熊猫参考基因组比对，得到比对结果为CNS文件。由于每个样本基因组存在一些低测序深度的区域，在此 综合所有样本基因型的可能性的文件，利用最大似然法整合所有样本的数据，产生包含所有样本每个位点的伪基因组(Pseudo-genome)。选择概率最大的基因型作为每个样本的一致基因型，通过基因型和测序深度等信息检测出高质量的SNPs。得到各个样本的一致性序列后，结果保存为群体SNPs格式，获得群体SNP数据。2)群体进化分析输入群体SNP结果，并以群体SNP为基础，整合调用多个软件进行群体进化分析，包括Tree，PCA，Structure和Frappe分析，具体如下。将软件命名PopuStruct.pl，相关参数说明如表3，需注意的是群体SNP文件必须与individual文件对应。Structure软件运行时间较长，如果时间紧急，建议首先用Frappe进行群体结构分析，得到初步分析结果。表3参数说明-indi<s>群体中每个个体信息，个体顺序与群体SNP文件一致，必须设置。-list<s>群体SNPsgenotype文件列表，必须设置。-OutDir<s>输出路径，默认当前路径。-prefix<s>输出脚本前缀信息，默认“Pop”。-Struct<y/n>是否用Structure软件进行群体结构分析，默认“y”-Tree<y/n>是否构建系统发生树，默认“y”-Frappe<y/n>是否用Frappe软件进行群体结构分析，默认“y”-PCA<y/n>是否进行主成分分析，默认“y”-queue<s>投递任务队列，默认bc.q-project<s>投递投任务-P参数值，默认rdtest-help帮助信息输出文件(结果)i)Frappe结果文件和Structure结果文件，可结合excel进行调整和作图。结果如图6所示，图6是Frappe基于群体SNP推测的群体遗传结构示意图，图中，分隔的每块代表一个群体，横坐标代表一个样本，不同分隔块代表K个不同或差异较大的祖先，分析每一个品系的遗传成分中，所具有的每一个假想祖先成分的比例。如果一个样品对应两个不同的分割块，则表示该样品可能是两个亚群之间的中间品种。当K值取得越大时，样品之间的差异性越被放大，分得越细，可根据实际结果来决定K值取到哪就可以完全体现出所有样品的结构关系。图中，K分别取2、3、4和5，可以看出K＝3即将群体分成3个亚群体基本可以完整体现出所有样本的结构关系。ii)tree结果文件利用mega软件进行调整，结果如图7所示。图7是基于群体SNP采用 邻接法推断的系统发生树的示意图，图中，分支距离越近，说明两分支间进化关系越近。对于同一亚群内的样本，应当显示能很好的分在一起或离得不远，通过该图可以说明品种之间的进化关系远近。从图7可看出，该群体可以分成3个亚群体。iii)PCA分析结果，需用excel进行作图，结果如图8所示。图8是基于群体SNP的PCA分析结果的示意图，图中不同形状的标记代表不同亚群的样本，一个标记点代表一个样品，点的横纵坐标分别是该样品对应的第一和第二特征向量中同一顺序元素的值，相应的特征值大小代表该主成分在整个关系中所占的比例，通过该图可以跟样品的实际分组进行对比，看出样品分组好坏。进而可以看要不要重新分类以获得新亚群。三、受选择作用位点的检测结合实施例一以及上述获得的群体SNP数据的结构，推导公式如下：FST=ΠBetween-ΠWithinΠBetween=1-ΠWithinΠBetween=1-[Σj(2nj)Σi2nijnij-1xij(1-xij)]/Σj(2nj)Σi2nini-1xi(1-xi)]]>上式中xij是SNP位点i在亚群体j中的次等位碱基(第二碱基)的频率；而nij是SNP位点i在亚群体j中染色体上的物理位置；nj则是亚群体j用于比较分析的SNP位点个数的总和。其中变量j依据上述群体结构分析结果，新取为3，变量i以最后判定的SNP位置代入。上述计算分析过程以群体SNP为基础，调用多个软件检测亚群体间可能存在的受到选择作用的位点，命名为SnpSelect.pl，使用的软件方法包括：Arlequin，BayesScan和Datacal三种，各软件对应参数说明，包括阈值的设置，详见表4。perlSnpSelect.pl<snp.list><individual><2population.txt>[options]；其中2population文件指的是参与位点选择分析的两个亚群信息，具体格式见说明。表4输出文件i)Arlequin分析结果，如图9所示。图9显示Arlequin程序基于群体SNP检测受选择作用位点的分析结果。横轴表示给定位点在群体水平的杂合度，纵轴表示亚群之间在给定位点上的杂合度差异值(Fst)。上部分圈起中的点表示受定向选择的位点(q<0.01或者q<0.05)，下部分圈起中的点表示受平衡选择的位点(q<0.01或者q<0.05)。ii)GlobalFSTtest分析结果，如图10所示。图10显示GlobalFSTtest程序基于群体SNP检测受选择作用位点的结果。横轴表示给定位点在群体水平的杂合度，纵轴表示亚群之间在给定位点上的杂合度差异值(Fst)。前1％Fst值所对应位点被认为是候选位点，即横线以上的点为检测出的受到选择作用的位点。iii)BayeScan分析结果，如图11所示。图11显示BayeScan程序基于群体SNP检测受选择作用位点的结果。横轴表示给定位点在群体水平的杂合度，纵轴表示将给定位点的检验q值(qvalue)取对数的值(以10为底数)。qvalue<0.1的位点被认为是候选受选择作用位点，即位于图上竖线右边的点为候选受选择作用位点。结合图9-图11，在位点选择分析时，获得有以上至少两种方法支持的判为最终的受选择作用位点。在本说明书的描述中，参考术语“一个实施例”、“一些实施例”、“示例”、“具体示例”、或“一些示例”等的描述意指结合该实施例或示例描述的具体特征、结构、材料或者特点包含于本发明的至少一个实施例或示例中。在本说明书中，对上述术语的示意性表述不一定指的是相同的实施例或示例。而且，描述的具体特征、结构、材料或者特点可以在任何的一个或多个实施例或示例中以合适的方式结合。尽管已经示出和描述了本发明的实施例，本领域的普通技术人员可以理解：在不脱离本发明的原理和宗旨的情况下可以对这些实施例进行多种变化、修改、替换和变型，本发明的范围由权利要求及其等同物限定。当前第1页1 2 3