一种基于无参考基因组的群体进化分析方法

摘要

本发明公开了一种基于2d‑RAD测序后无参考基因组的群体进化分析方法，通过将样本进行数据拆分后过滤聚类得到群体SNP，并基于样本分组和群体SNP信息进行群体遗传参数分析后构建系统发育树，确定最佳K值后利用所述R自写脚本，根据群体SNP信息和指定的群体信息来寻找两大Group之间的共有和特有SNP信息进行无参考基因组的群体进化分析。本发明的整个数据分析比较自动化，节省了人力成本，提高了分析效率，避免了可能的人为失误且分析的数据图表更加美观。

说明书

技术领域

本发明涉及基因测序分析技术领域，具体涉及一种基于无参考基因组的群体进化分析方法。

背景技术

通过群体进化分析能更加深入的探究同物种内不同亚群之间的群体结构差异、基因交流情况，也能够研究不同物种之间的群体结构特征；但很多的物种还没有参考基因组发表，所以就要进行无参考基因组的群体进化分析。

因为无参建库方法有多种(RAD、GBS、2d-RAD、SLAF等)，不同的建库方法在无参分析的第一步数据拆分上会不同，而现有的基于2d-RAD建库的无参分析方法数据过滤流程复杂，效率较低，尤其是项目数量多且一个项目中包含的样本量大时，实际操作过程中一个项目可能会多次上机测序，这样就会得到不同批次的数据，现有无参分析方法无法智能地使用自动化流程合并不同批次的数据并且进行过滤，导致数据合并与过滤会耗费大量人工时间。

随着高通量测序的不断发展，已有的分析流程分析内容显得单薄，分析内容较少，新的无参分析内容种类更加多样化和个性化。以往的无参分析流程中有许多地方需要人工进行操作进行，现在新的无参分析方法更加自动化，该自动化流程提高服务器使用效率，减少分析人员的分析压力，便于控制分析内容。

发明内容

为了克服现有技术的上述缺陷，本发明的目的在于提供一种基于无参考基因组的群体进化分析自动化分析方法。

为了实现本发明的目的，所采用的技术方案是：

一种基于2d-RAD测序后无参考基因组的群体进化分析方法，包括如下步骤：

第一步：根据测序样本中的barcode、酶1和酶2的酶切位点信息利用拆分脚本进行数据拆分后将同个样本的多个下机测序数据进行合并后以fastq.gz的格式保存在第一文件夹中；

第二步：将第一步拆分合并后的数据通过滤脚本对数据进行fastQC的质控后按照碱基质量值：Q≥20和序列长度≥50bp的标准进行数据过滤得到过滤后的数据以fastq.gz的格式保存在第二文件夹中；

第三步：将单个样本内先进行序列聚类，在聚类前将单个样本的双端测序数据合并到一个文件中，然后利用软件Stacks中的ustacks命令进行聚类，得到每个样本的代表catalog序列，结果文件以tags.tsv.gz的格式保存在第三文件夹中；

第四步：将样本分组后，基于单个样本的catalog序列进行聚类得到所有样本的consensus序列，所述consensus序列为用于所有样本的类参考基因组序列；

第五步：读取所有文件指定每个样本的分组信息，同时指定缺失率参数利用软件Stacks中的cstacks命令检测群体SNP信息，并以VCF文件的格式保存群体SNP信息；

第六步：基于第五步的群体SNP信息，利用Stacks中的populations命令进行群体遗传参数的分析，计算得到群体分化指数Fst、群体核苷酸多样性π、群体期望杂合度和观测杂合度、单倍型多样性数据；

第七步：将第五步的群体SNP信息的VCF文件，使用软件vcftools和plink进行格式转换，然后使用软件GCTA对SNP进行降维分析，得到对群体影响较大的三个主成分并计算各主成分的贡献值，最后用R自写脚本绘制PCA分布图；

第八步：利用Python自写脚本将得到的群体SNP信息与单个样本的SNP信息转换格式连在一起，然后利用不同的模型构建系统发育树；

第九步：

利用Perl自写脚本将群体SNP格式转换为软件structure要求的格式，然后指定分析中使用的SNP数目、群体数目，计算每个样本属于指定的祖先百分比；

然后确定最佳K值(祖先个数)，通过这个结果可以得到样本的分群信息和最初指定的是否一致；

第十步：

利用Perl自写脚本，根据群体SNP信息和指定的群体信息来寻找两大Group之间的共有和特有SNP信息。

在本发明的一个优选实施例中，所述过滤脚本为filter_batch_v2.pl。

在本发明的一个优选实施例中，所述构建系统发育树的模型包括maximumparsimony(MP)、neighbor-joining(NJ)、Maximum Likelihood(ML)或Bayesian method(BI)中的任意一种或多种。

在本发明的一个优选实施例中，所述最佳K值为ln likelihood进入平台期后拐点所对应的K值。

本发明的有益效果在于：

整个数据分析比较自动化，节省了人力成本，提高了分析效率，避免了可能的人为失误且分析的数据图表更加美观。

附图说明

图1为本发明的流程图。

图2为本发明的PCA分布图。

图3为本发明基于NJ模型的进化树分布图。

图4为本发明最佳K＝3时群体遗传结构分布图。

具体实施方式

本发明的原理在于：

基于2d-RAD无参简化的自动化过滤流程，可进行批量数据拆分和过滤，数据过滤后续的各项分析也都可以自动化完成，提高数据处理效率和服务器使用效率，节约人工时间，同时降低人为错误，最终缩短了整个项目分分析周期，实现了丰富分析内容的无参分析高效自动化。

结合图1，本发明的一种基于2d-RAD测序后无参考基因组的群体进化分析方法，包括如下步骤：

(1)数据拆分步骤

根据测序样本的barcode、酶1和酶2的酶切位点信息用自写的脚本进行数据自动化拆分，格式大致为一行表示一个样本的信息，每列的元素分别为样本名，barcode碱基，酶1的酶切位点，酶2的酶切位点，其中间隔符设置为制表符；若一个样本有多次下机测序，分析流程会自动匹配进行合并，合并后的数据以fastq.gz的格式被统一存放在1_RawData的文件夹中。

拆分脚本具体是：

一个文库内会包含多个样本，将样本名、barcode、酶1和酶2的酶切位点序列这四列作为输入文件1，将文库下机原始双端数据fastq.gz作为输入文件2和3；

若一条序列的R1的5’端前7bp与barcode一致，接下来的5个碱基与酶1酶切位点一致，且该reads对应的R2的5’端前4bp与酶2的酶切位点序列一致，则会把这条reads拆分到该样本中，多次循环并输出最终每个样本的拆分数据结果。

(2)数据质控和过滤步骤

利用自己写的自动化过滤脚本filter_batch_v2.pl对样本进行fastQC的质控，同时按照碱基质量值(Q≥20)和序列长度(≥50bp)的标准进行数据过滤。运行结束后所有高质量数据以fastq.gz的格式被存放在2_HQData中。

过滤脚本为filter_batch_v2.pl：

该脚本首先读取位于1_RawData中的样本下机原始数据的双端序列文件${name}_R1.fastq.gz和${name}_R1.fastq.gz作为输入文件，然后对文件进行重命名，通过软件fastqc对输入文件进行质控，得到原始数据的碱基质量等信息的fastq文件；

再利用软件AdapterRemoval以原始数据的fastq.gz文件作为输入文件，去除测序接头，同时将新产生的结果文件以fastq格式保存在2_HQData中，然后再将上一步新产生的fastq文件作为序列质量过滤程序的输入文件，采用滑动窗口法进行质量过滤，窗口大小设置为5bp，步长设置为1bp；

每一次往前移动一个碱基，取5个碱基计算窗口的平均Q值，若窗口的平均Q值≤20，则仅保留该窗口倒数第二个碱基及之前的碱基；

之后双末端中任意一条reads的长度≤50bp，则去除该双末端reads。最终结果输出为${name}_HQ_R1.fq和${name}_HQ_R2.fq。

(3)单个样本内序列聚类步骤

因为无参分析没有参考基因组，所以单个样本内先进行序列聚类，在聚类前将单个样本的双端测序数据合并到一个文件中，然后利用软件Stacks中的ustacks命令进行聚类，得到每个样本的代表catalog序列，结果文件以tags.tsv.gz的格式存放在3_Stacks文件夹中。

(4)所有样本的catalog序列聚类步骤

指定样本的分组信息，同时基于单个样本的catalog序列进行聚类得到所有样本的consensus序列，这个consensus序列被当做所有样本的类参考基因组序列。

(5)检测群体SNP的步骤

读取所有文件指定每个样本的分组信息，同时指定缺失率参数利用软件Stacks中的cstacks命令检测群体SNP信息，并以VCF文件的格式保存。

(6)群体遗传参数(Fst、π、杂合度、单倍型多样性)分析

根据群体SNP信息，利用Stacks中的populations命令进行群体遗传参数的分析，计算得到群体分化指数Fst、群体核苷酸多样性π、群体期望杂合度和观测杂合度、单倍型多样性。

(7)群体PCA分析的步骤

根据群体SNP的VCF文件使用软件vcftools和plink进行格式转换，然后使用软件GCTA对SNP进行降维分析，得到对群体影响较大的三个主成分并计算各主成分的贡献值，最后用R自写脚本绘制PCA分布图。

R自写脚本首先读取GCTA软件输出的PC1、PC2向量信息做为输入文件，计算PC1和PC2的贡献率，然后利用R中的ggplot2包画散点图。

(8)群体系统发育树分析的步骤

自写脚本将得到的群体SNP信息将每个样本的SNP信息转换格式连在一起，然后利用不用的模型构建系统发育树。

构建进化树的常用模型包括maximum parsimony(MP)、neighbor-joining(NJ)、Maximum Likelihood(ML)、Bayesian method(BI)；

其中，MP模型适用于位点不存在回复突变和平行突变，序列相似度较高，核苷酸或氨基酸数目大、替代速率稳定的长序列。NJ模型适用于进化距离不大，信息位点少的短序列。在进化模型确定的情况下，ML法是与进化事实吻合最好的建树方法。BI模型保留了最大似然法的基本原理，又引进了马尔科夫链的蒙特卡洛方法，适用于推导系统树、评估系统树的不确定性、检测选择作用、比较系统树、参考化石记录计算分歧时间和检测分子钟。

(9)群体遗传结构分析的步骤

自写脚本将群体SNP格式转换为软件structure要求的格式，然后指定分析中使用的SNP数目、群体数目，计算每个样本属于指定的祖先百分比。然后确定最佳K值(祖先个数)，通过这个结果可以得到样本的分群信息和最初指定的是否一致。

每个K值基于贝叶斯模型的计算方法模拟的结果，都会产生对应最大似然值(likelihood)，它是取了自然对数后输出的(ln likelihood)。ln likelihood越大，说明K值越接近于真实情况，但一般随着K值升高，ln likelihood值也会进入平台期。最优K值就是进入平台期的那个拐点对应的K值)。

(10)群体特有SNP分析的步骤

自写脚本根据群体SNP信息和指定的群体信息来寻找两大Group之间的共有和特有SNP信息。

首先根据基因型缺失的情况和SNP位点的测序深度对原始SNP进行过滤，群体SNP的特异性通过两个阈值(A和B)来定义，一是SNP在目标群体中出现的频率高于某阈值(A)，二是SNP在非目标群体中出现的频率低于某阈值(B)，阈值一般设置为0.8。

基于上述步骤本发明的优点在于：

(1)整个数据分析比较自动化，节省了人力成本，提高了分析效率，避免了可能的人为失误。

(2)分析内容更加丰富，分析结果的图形更加美观(如图2-4所示)。