一种批量合并eQTL分析结果中存在连锁不平衡SNP的方法

摘要

本发明公开了一种批量合并eQTL分析结果中存在连锁不平衡SNP的方法。本发明所提供的合并eQTL分析结果中存在连锁不平衡SNP的方法基于SNP与靶基因的位置信息将其分为顺式和反式SNP，基于SNP位置信息合并相邻的SNP为一个SNP cluster，基于SNP cluster间的连锁不平衡程度进一步合并结果。本发明的脚本由python3语言写成，速度快，灵活性高，可靠性强，实现了批量化、自动化和流程化计算。本发明将在eQTL分析结果的简化上发挥重要作用。

说明书

技术领域

本发明属于生物技术领域，涉及一种批量合并eQTL分析结果中存在连锁不平衡SNP的方法。

背景技术

关联分析(association analysis)是一种以位点间的连锁不平衡(linkagedisequilibrium)为基础，采用统计方法检测遗传多态性与性状之间关联的分析方法。全基因组关联分析(genome-wide association study，GWAS)最初应用于人类疾病相关的研究中，这些研究对人类理解相关疾病的遗传基础和分子机制具有显著的贡献。近年来，随着高密度SNP基因分型芯片和全基因组测序等技术的发展，GWAS广泛应用于作物复杂性状遗传结构的解析。与传统的连锁分析相比，GWAS有多个优势：可以利用自然群体，无需针对特定性状构建作图群体，花费时间少；能同时检测群体同一基因的多个等位基因，有利于优良等位基因的挖掘；当群体连锁不平衡程度低并且标记覆盖度高时，定位精度高，可以达到单基因水平。

eQTL(expression quantitative trait loci)分析是将每个基因的表达量作为表型，将它们与基因组变异位点进行关联分析，从而发掘调控基因表达的遗传变异的分析方法。通过eQTL分析，不但可以获取基因的表达调控位点，还可以构建基因之间的调控网络，对于解析基因的调控机制具有很好的指导作用。目前，eQTL分析已经成功地应用到多个物种中，包括拟南芥、玉米和水稻等。根据eQTL位点与调控的目标基因的位置关系可以将eQTL分成顺式eQTL(cis eQTL)和反式eQTL(trans eQTL)。顺式eQTL与目标基因物理距离较近，表明是该基因本身的差别引起mRNA水平的差别，反式eQTL与目标基因物理距离较近或位于不同的染色体，表明是其他基因的差别控制该基因mRNA水平的差异。

做eQTL分析前首先要对基因的表达谱进行正态化，其次利用正态化的表达谱计算隐性因素，利用基因型数据计算群体结构，最终以群体结构和隐性因素作为协变量，对群体的基因表达谱和基因型数据进行关联分析。由于得到的eQTL结果非常多，常常数以十万、百万计，而相邻的SNP间往往连锁不平衡(linkage disequilibrium，LD)程度高，因此需要合并冗余SNP简化结果。手动合并存在LD的SNP非常耗时耗力，但目前没有软件能提供批量合并的功能，这成为一个亟待解决的问题。

发明内容

本发明提供的是一种批量合并eQTL分析结果中存在连锁不平衡SNP的方法，具体包括如下步骤：

(1)在Windows操作系统下创建工作目录eqtl_analysis及其子文件夹gene_info，将待分析植物基因注释信息文件***.gff3和脚本abstract_gene_info.py放在gene_info文件夹下，运行“python abstract_gene_info.py***.gff3”命令，得到每条染色体各自的基因信息文件，记为G数据集。

所述“***.gff3”为研究物种的基因注释文件，油菜Darmor-bzh参考基因组对应的文件为Brassica_napus.annotation_v5.gff3，油菜中双11参考基因组对应的文件为ZS11.annotation.gff3。

G数据集文件命名方式为染色体名+“_gene_info.txt”，不保留标题行，以油菜为例，染色体A01基因型文件为A01_gene_info.txt。文件包括5列，分别为基因名、染色体、基因起始位置、基因中止位置和正负链信息(图1)。

(2)将待分析植物的eQTL结果文件记为A数据集，格式为eQTL分析常用软件MatrixeQTL的结果文件格式；脚本separate_cis_trans.py参考G数据集提供的基因物理位置，分析SNP与基因的染色体和物理距离，将所有SNP划分为两类，顺式SNP和反式SNP。A数据集和脚本eqtl_cis_trans.py均放在工作目录eqtl_analysis下，运行命令“pythonseparate_cis_trans.py XXX1.txt cis_dis”，得到“XXX1_cis.txt”和“XXX1_trans.txt”两个文件。

“XXX1.txt”代表所述A数据集的文件名，格式与eQTL分析常用软件MatrixeQTL的结果文件格式相同，包含6列“SNP”，“gene”，“beta”，“t-stat”，“p-value”和“FDR”(图2)，脚本利用其中“SNP”，“gene”和“p-value”三列信息。文件需按“gene”和“SNP”两列信息进行排序。SNP的命名方式用染色体名+物理位置，染色体为3位或者10位，如A03和C03_random；物理位置为8位数，不足的位数用0补全。以油菜为例，SNP的染色体为A03，物理位置为41437，SNP命名为A0300041437；SNP的染色体为C03_random，物理位置为44754，SNP命名为A01_random00044754。

“cis_dis”为划分SNP为顺式SNP和反式SNP的距离阈值，默认设定为24,000bp。

所述“XXX1_cis.txt”为包含所有顺式SNP的文件名，记为B1数据集；“XXX1_trans.txt”为包含所有反式SNP的文件名，记为B2数据集。

(3)脚本combine_near_snp.py合并相邻的显著SNP，得到SNP cluster，并用其中最显著、物理位置小的SNP作为代表。将脚本combine_near_snp.py放在工作目录eqtl_analysis下，针对B1、B2数据集分别运行命令“python combine_near_snp.py XXX1_cis.txt part_dis”和“python combine_near_snp.py XXX1_trans.txt part_dis”，得到“XXX1_cis_median.txt”和“XXX1_trans_median.txt”两个文件。

“part_dis”为合并相邻SNP的距离阈值，默认设定为10,000bp。所述“XXX1_cis_median.txt”为合并相邻的顺式SNP后得到的结果文件，记为C1数据集；“XXX1_trans_median.txt”为合并相邻的反式SNP后得到的结果文件，记为C2数据集。

(4)为计算SNP cluster间的LD系数r

M数据集各染色体基因型文件命名方式为染色体名+“_snp_info.txt”，不保留标题行，以油菜为例，染色体A03基因型文件为A03_snp_info.txt。文件格式为SNP+基因型信息，SNP只能包含两个等位基因，分别用0和2表示，杂合、缺失用NA表示(图3)。

所述“XXX1_cis_final.txt”为合并所有相邻、存在LD的顺式SNP的最终结果文件，记为D1数据集；所述“XXX1_trans_final.txt”为合并所有相邻、存在LD的反式SNP的最终结果文件，记为D2数据集。

在上述方法步骤(1)中，所述脚本abstract_gene_info.py关于提取基因信息的内容，是基于如下原理进行编程的：提取基因注释文件中第三列为“mRNA”的所有行，提取第九列“Parent”后内容作为基因名，第一列内容作为染色体信息，第四、五和七列作为起始、中止和正负链信息。

在上述方法步骤(2)中，所述脚本separate_cis_trans.py关于划分SNP为顺式和反式的内容，是基于如下原理进行编程的：根据G数据集，提取基因染色体和物理位置；若SNP与基因位于不同染色体，划分为反式SNP；SNP与基因位于相同染色体时，若SNP位于基因上下游24kb内，划分为顺式SNP，否则划分为反式SNP。

在上述方法步骤(3)中，所述脚本combine_near_snp.py关于合并相邻SNP为SNPcluster的内容，是基于如下原理进行编程的：判断同一染色体物理位置最大和最小的SNP距离是否超过10kb；距离小于10kb时，合并所有SNP为一个cluster，用cluster中最显著、物理位置最小的SNP作为代表；距离大于10kb时，以10kb为间距，将SNP分为多组，每组分别进行合并，得到多个cluster。

在上述方法步骤(4)中，所述脚本combine_ld_snp.py关于合并SNP cluster的内容，是基于如下原理进行编程的：参考M数据集，计算SNP cluster间的LD系数r

在本发明的一个实施例中，所述脚本abstract_gene_info.py具体为：

在本发明的一个实施例中，所述脚本separate_cis_trans.py具体为：

在本发明的一个实施例中，所述脚本combine_near_snp.py具体为：

在本发明的一个实施例中，所述脚本combine_ld_snp.py具体为：

在所述方法中，所述基因注释文件可以通过下载已公开的物种基因注释文件获得。

具体地，本发明所述物种具体为甘蓝型油菜，所述基因注释文件记录于法国GenoScope网站(https://www.genoscope.cns.fr/brassicanapus/)。

本发明的优点在于：1.可批量合并eQTL分析结果中存在连锁不平衡SNP，策略为先合并相邻位点再合并存在LD的位点，速度快，效率高；2.脚本可以自定义合并相邻SNP的区间大小和划分顺式/反式SNP的距离阈值，灵活性高；3.当输入文件很大时，脚本需要按分段读取文件，会造成部分跨段且存在LD的SNP无法合并，导致结果的冗余。本方法考虑到该问题并在脚本中添加了相应代码予以解决，确保分析结果的代表性。

附图说明

图1为G数据集文件格式，包括5列，分别为基因名、染色体、基因起始位置、基因中止位置和正负链信息。

图2为A数据集文件格式，包含6列“SNP”，“gene”，“beta”，“t-stat”，“p-value”和“FDR”。

图3为M数据集文件格式，包含SNP名和基因型信息。

图4为本发明批量合并eQTL分析结果中存在连锁不平衡SNP的方法的流程图。

具体实施方式

下述实施例中所使用的实验方法如无特殊说明，均为常规方法。

下述实施例中所用的材料，如无特殊说明，均可从商业途径得到。

实施例1、批量合并eQTL分析结果中存在连锁不平衡SNP方法的建立

本发明所提供的批量合并eQTL分析结果中存在连锁不平衡SNP方法的流程图见图4，具体包括如下步骤：

(1)在Windows操作系统下创建工作目录eqtl_analysis及其子文件夹gene_info，将基因注释信息文件***.gff3和脚本abstract_gene_info.py放在gene_info文件夹下，运行“python abstract_gene_info.py***.gff3”命令，得到每条染色体各自的基因信息文件，记为G数据集。

所述“***.gff3”为研究物种的基因注释文件，油菜Darmor-bzh参考基因组对应的文件为Brassica_napus.annotation_v5.gff3，油菜中双11参考基因组对应的文件为ZS11.annotation.gff3。

G数据集文件命名方式为染色体名+“_gene_info.txt”，不保留标题行，以油菜为例，染色体A01基因型文件为A01_gene_info.txt。文件包括5列，分别为基因名、染色体、基因起始位置、基因中止位置和正负链信息。其中，所属脚本abstract_gene_info.py具有如下特点：提取基因注释文件中第三列为“mRNA”的所有行，提取第九列“Parent”后内容作为基因名，第一列内容作为染色体信息，第四、五和七列作为起始、中止和正负链信息。

(2)将待分析的eQTL结果文件记为A数据集，格式为eQTL分析常用软件MatrixeQTL的结果文件格式；脚本separate_cis_trans.py参考G数据集提供的基因物理位置，分析SNP与基因的染色体和物理距离，将所有SNP划分为两类，顺式SNP和反式SNP。A数据集和脚本eqtl_cis_trans.py均放在工作目录eqtl_analysis下，运行命令“python separate_cis_trans.py XXX1.txt cis_dis”，得到“XXX1_cis.txt”和“XXX1_trans.txt”两个文件。所述“XXX1.txt”代表所述A数据集的文件名，格式与eQTL分析常用软件MatrixeQTL的结果文件格式相同，包含6列“SNP”，“gene”，“beta”，“t-stat”，“p-value”和“FDR”，脚本利用其中“SNP”，“gene”和“p-value”三列信息。文件需按“gene”和“SNP”两列信息进行排序。SNP的命名方式用染色体名+物理位置，染色体为3位或者10位，如A03和C03_random；物理位置为8位数，不足的位数用0补全。以油菜为例，SNP的染色体为A03，物理位置为41437，SNP命名为A0300041437；SNP的染色体为C03_random，物理位置为44754，SNP命名为A01_random00044754。所述“cis_dis”为划分SNP为顺式SNP和反式SNP的距离阈值，默认设定为24,000bp。

所述“XXX1_cis.txt”为包含所有顺式SNP的文件名，记为B1数据集；“XXX1_trans.txt”为包含所有反式SNP的文件名，记为B2数据集。其中，所述脚本separate_cis_trans.py具有如下特点：根据G数据集，提取基因染色体和物理位置；若SNP与基因位于不同染色体，划分为反式SNP；SNP与基因位于相同染色体时，若SNP位于基因上下游24kb内，划分为顺式SNP，否则划分为反式SNP。

(3)脚本combine_near_snp.py合并相邻的显著SNP，得到SNP cluster，并用其中最显著、物理位置小的SNP作为代表。将脚本combine_near_snp.py放在工作目录eqtl_analysis下，针对B1、B2数据集分别运行命令“python combine_near_snp.py XXX1_cis.txt part_dis”和“python combine_near_snp.py XXX1_trans.txt part_dis”，得到“XXX1_cis_median.txt”和“XXX1_trans_median.txt”两个文件。

所述“part_dis”为合并相邻SNP的距离阈值，默认设定为10,000bp。

所述“XXX1_cis_median.txt”为合并相邻的顺式SNP后得到的结果文件，记为C1数据集；“XXX1_trans_median.txt”为合并相邻的反式SNP后得到的结果文件，记为C2数据集。

其中，所述脚本combine_near_snp.py具有如下特点：判断同一染色体物理位置最大和最小的SNP距离是否超过10kb；距离小于10kb时，合并所有SNP为一个cluster，用cluster中最显著、物理位置最小的SNP作为代表；距离大于10kb时，以10kb为间距，将SNP分为多组，每组分别进行合并，得到多个cluster。

(4)为计算SNP cluster间的LD系数r

M数据集各染色体基因型文件命名方式为染色体名+“_snp_info.txt”，不保留标题行，以油菜为例，染色体A03基因型文件为A03_snp_info.txt。文件格式为SNP+基因型信息，SNP只能包含两个等位基因，分别用0和2表示，杂合、缺失用NA表示。

所述“XXX1_cis_final.txt”为合并所有相邻、存在LD的顺式SNP的最终结果文件，记为D1数据集；所述“XXX1_trans_final.txt”为合并所有相邻、存在LD的反式SNP的最终结果文件，记为D2数据集。

其中，所述脚本combine_ld_snp.py具有如下特点：参考M数据集，计算SNPcluster间的LD系数r

实施例2、利用实施例1建立的方法批量合并甘蓝型油菜eQTL分析结果中存在连锁不平衡SNP

进入GenoScope网站下载甘蓝型油菜Darmor-bzh的基因注释文件Brassica_napus.annotation_v5.gff3。我们之前对29个甘蓝型油菜的分枝进行转录组测序，获得基因的表达谱数据，结合已获得的基因型数据，利用MatrixeQTL软件进行eQTL分析，结果文件包含8,489个基因及其5，616，239个显著SNP，命名为“29_all_snp.txt”(A数据集)。利用本发明中的脚本对该文件进行SNP的合并。分。具体包括如下步骤：

1)参照实施例1的步骤(1)进行。

采用abstract_gene_info.py脚本从甘蓝型油菜基因注释文件Brassica_napus.annotation_v5.gff3中提取每条染色体的基因信息，分别生成文件(G数据集)保存在gene_info子文件夹中。

2)参照实施例1的步骤(2)进行。

利用separate_cis_trans.py脚本对A数据集进行分析，参考G数据集提供的基因物理位置，分析SNP与基因的染色体和物理距离，将所有SNP划分为两类，顺式SNP(B1数据集)和反式SNP(B2数据集)。最终检测到203,868顺式SNP和5,412,371反式SNP。

3)参照实施例1的步骤(3)进行。

采用combine_near_snp.py脚本分别对B1和B2数据集进行相邻SNP(<10kb)的合并，得到6,324顺式(C1数据集)和490,971反式SNP cluster(C2数据集)。

4)参照实施例1的步骤(4)进行。

提取C1和C2数据集所有SNP的名称，从29份材料的基因型文件(vcf格式)中提取这些SNP的基因型信息，按说明书部分的描述整理成相应的格式，存放在snp_info子文件夹中(M数据集)。脚本combine_near_snp.py根据M数据集信息，分别计算C1和C2数据集中同一靶基因不同SNP cluster间的LD系数r

随机挑选10个基因SNP的eQTL结果进行验证，与脚本分析得到的结果完全相同，这证实了本方法的可靠性。