基于亲本基因型与子代表型的全基因组关联分析算法
阅读说明：本技术 基于亲本基因型与子代表型的全基因组关联分析算法 (Whole genome association analysis algorithm based on parental genotype and progeny phenotype ) 是由 夏晓勤 石米娟 张婉婷 程莹寅 于 2021-09-06 设计创作，主要内容包括：本发明公开了一种基于亲本基因型与子代表型的全基因组关联分析算法,包括以下步骤：获取待分析群体的亲本基因型、子代表型以及子代父母本信息,根据亲本对的基因型建立子代的组合基因型矩阵,获得每个SNP位点不同子代的可能基因型组合信息,以及对应的子代群体分组表型数据,构建子代组合基因型与对应子代表型的统计模型进行关联分析,获取表型与各SNP位点的关联P值；之后区分SNP位点的关联类型,计算是否符合加性或完全显性效应,对候选位点进行筛选；再根据全基因组SNP标记的LD衰减获得强相关标记的平均距离,用于最终标记集合的筛选。(The invention discloses a whole genome association analysis algorithm based on parental genotypes and offspring phenotypes, which comprises the following steps: acquiring parent genotypes, offspring phenotypes and offspring parental information of a population to be analyzed, establishing an offspring combined genotype matrix according to the genotypes of parent pairs, acquiring possible genotype combined information of different offspring of each SNP locus and corresponding offspring population grouping phenotype data, constructing a statistic model of the offspring combined genotypes and the corresponding offspring phenotypes for association analysis, and acquiring association P values of the phenotypes and the SNP loci; then distinguishing the association type of the SNP locus, calculating whether the additive or complete dominant effect is met, and screening candidate loci; and then obtaining the average distance of the strong related markers according to the LD attenuation of the whole genome SNP marker, and using the average distance for screening a final marker set.)
技术领域
本发明属于生物信息学技术领域，具体为一种基于亲本基因型与子代表型的全基因组关 联分析算法。
背景技术
基因序列往往包含着不同数量级的序列变异(如人类基因包含着百万种序列变异)，即 单核苷酸多态性(SNP)，这些变异对于疾病的形成、生长的进程、药物的反应等性状有着 直接或间接的影响。全基因组关联分析(Genome-Wide Association Studies，简称GWAS)是 指应用基因组中的SNP为分子遗传标记，对大量个体在全基因组范围的遗传标记进行检测， 获得基因型，进而与可观测的性状(表型)进行相关性分析，筛选出与性状相关的SNP， 挖掘出影响性状的基因变异。
现有的GWAS方法分析的是同一个体的基因型与表型构成的对应关系，为了获取群体 的基因型信息，需要对大量的单样本进行全基因组重测序以获取各样本的标记信息，同时采 集各样本的表型数据。通过在全基因组范围内分析大量个体的基因型－表型对应关系，确定 某些基因座与表型之间的关联关系。通常样本量越大，获得的关联基因信息准确率越高，一 般实验样本数量从数百到数万不等，需耗费非常大的实验与测序成本。
实际上，性状虽然由个体本身的基因型决定，但其基因型与亲本的基因型是存在关系的。 对于同一基因座，在全同胞子代群体中只可能存在有限的几种基因型类型，且这些基因型的 占比符合孟德尔分离比，即一群全同胞群体的单一基因座的基因型组合是可以由对应的亲本 的基因型推算出来的。原理上，当子代数量较多，即群体抽样数目足够时，这种组合和子代 群体表型的关系是可以用来进行关联分析且获得可靠结果的。鱼类通常能产生数目巨大的子 代群体，且体外受精模式也使得大规模泛交群体的构建十分容易，能恰好满足这种关联分析 所需的条件。因此本发明根据鱼类数据的特有繁殖体系特点，利用亲本的基因型与子代表型 数据，基于广义线性模型(GLM)，检测单个位点的加性显性效应，并结合标记的连锁不平 衡等信息，构建了新型全基因组关联分析算法。本发明的方法仅需对亲本进行重测序而不需 要对大量子代进行测序，极大地减少了一般关联分析中的实验及测序成本，且关联的标记位 点原理上不存在家系效应。同时本方法还可以获得亲本的全基因组SNP、InDel标记库作为 副产物。
发明内容
本发明的目的在于提供一种基于亲本基因型和子代表型的全基因组的关联算法，其优点 在于仅需对亲本进行重测序而不需要对大量子代进行测序，极大地减少了一般关联分析中的 测序成本，且关联的标记位点原理上不存在家系效应，同时本方法还可以获得亲本的全基因 组SNP、InDel标记库作为副产物。
为了实现上述目的，本发明采用以下的技术方案：
基于亲本基因型与子代表型的全基因组关联分析算法，包括以下步骤：
S1)对待分析亲本进行全基因组重测序；
S2)获取子代的表型数据(如体重、体长等性状数据)；
S3)使用BWA、Bowtie2等比对软件对S1)中亲本的全基因组重测序数据与其物种基因组进行比对；
S4)对得到的比对结果使用SNP检测工具(如GATK、Samtools、BCFtools等)进行SNP检测(SNP信息保存为VCF文件格式)，获得全基因组范围的SNP信息；
S5)从VCF文件中提取亲本SNP位点的基因型信息，得到SNP位点在亲本中的分型矩阵；
S6)接着使用下述原理和方法进行亲本基因型和子代表型的关联分析：
S6.1)过滤子代表型异常值；
S6.2)依据子代的亲本信息，结合SNP位点在亲本中的分型矩阵，获得每个SNP位点不同子代的可能基因型组合信息，以及对应的子代群体分组表型数据；
S6.3)建立子代组合基因型与表型的统计模型(如广义线性模型(GLM)、混合线性模 型(MLM)等)，进行基因型与表型的关联分析，获取表型与各SNP位点的关联P值，以 此初步确定是否为显著关联位点；此时，如果某位点的单个组合基因型中含有无法分型的情 况(即“NN”)，则剔除该处无法分型的值后再进行关联分析，同时，去除所有子代基因型 组合相同(子代基因型没有产生分离)；
S6.4)依据加性显性理论及统计检验原理进行SNP位点的加性效应和显性效应分析；
进一步地，位点加性效应和显性效应分析方法具体如下：
首先统计SNP位点的子代基因型组合信息，具体包括各SNP位点子代基因型组合的种 类数、子代基因型组合、各基因型组合对应的子代数目、各基因型组合对应的子代表型平均 值，然后，结合子代基因型组合及其对应子代表型信息，运用下述方法检验位点是否具有加 性或显性效应：
位点加性效应是指各基因型组合对应子代表型具有以下特征：
W(1/1×1/1)>W(0/1×1/1)>W(0/0×1/1)≈W(0/1×0/1)>W(0/0×0/1)>W(0/0×0/0)
位点显性效应是指各基因型组合对应子代表型具有以下特征：
W(1/1×1/1)≈W(0/1×1/1)≈W(0/0×1/1)>W(0/1×0/1)>W(0/0×0/1)>W(0/0×0/0)
式中的“0”表示某一条染色体上，某一SNP位点与参考基因组一致的分型，“1”表示与参考基因组不一致的分型，因此二倍体“0/0”表示两条染色体与参考基因组均一致的纯合子分型，“1/1”表示两条染色体均不一致的纯合子分型，“0/1”表示一条分型一致、另一条不一致的杂合子分型；
“≈”表示符号左右两组子代表型经统计检验没有显著性差异，“＞”表示符号左边的子 代表型均值大于右边的子代表型均值，且符号左右两边的子代表型数据经统计检验具有显著 性差异；
同一个位点可能同时具有加性效应和显性效应，例如，子代在某一个位点的基因型组合 只有0/0×0/1和0/0×0/0两种类型，0/0×0/1基因型组合的子代表型均值大于0/0×0/0基因型 组合的子代表型均值，且经过统计检验验证0/0×0/1类型的子代表型与0/0×0/0类型的子代 表型具有显著性差异，则认为其满足W(0/0×0/1)>W(0/0×0/0)的特征，依据上述公 式及说明知此位点既具有加性效应又具有显性效应；
S6.5)利用软件(如PopLDdecay)检测全基因组SNP位点的连锁不平衡衰减(LDdecay) 效应，确定超强相关标记(如r2>0.9)的平均距离DLD，若最大的相关系数r2对应的平均距离 小于10，则将DLD确定为10；
S6.6)使用snpEFF软件分析各位点在基因组上的位置，获得位点的基因组注释信息， 对所有SNP按照所属基因进行分组；
S6.7)筛选显著关联基因，具体分为以下步骤：
①设定距离阈值K，K设置为小于5000的正数，获取基因上以及距离基因上下游Kbp 以内的SNP位点；
②筛选出S6.4)中标记出的具有加性或显性效应的位点；
③确定关联P值的阈值p0，依据S6.3)中计算出的各位点关联P值，筛选出以基因为分组单元的各组区域内连续2个及以上位点关联P值小于p0的区域块；
④计算上述区域块上起始位点与终止位点之间的距离D(即多标记间的距离)，设置距 离阈值d＝max{350，DLD}，筛选出D大于d的区域块，将这些区域块确定为候选关联标记；
⑤计算候选关联标记上各位点关联P值的几何均值Pm，设置几何均值的阈值p1，从候选关联标记中筛选出Pm小于p1的标记，获得关联标记，根据关联标记的注释信息，获 得关联基因。
与现有技术相比，本发明具有以下优点及有益效果：由于不似现有技术获取的都是同一 个体的基因型与表型，而是对亲本的基因型和子代的表型建立对应关系，所以不需对数百乃 至数万的个体进行重测序，而只需对少量亲本进行重测序，与子代的表型进行关联，极大的 减少了人力物力及时间成本。并且由于直接利用亲本的基因型，在获取的结果中囊括了亲本 的所有基因型，能构建非常完善且有效的分子标记库。在方法上考虑了加性与显性效应，采 用广义线性模型、T-检验等相结合的统计方法，可以直接获得性状高关联的单体型及其基因。
附图说明
图1为1729尾子代体重箱形图，纵轴表示体重(单位：g)。
图2为本发明的主要算法流程图。
图3为亲本基因型与子代表型结合的方法示意图。
具体实施方式
接下来以草鱼样本数据对本发明所描述的基于亲本基因型和子代表型的关联分析方法 进行说明。
S1)对本实验室2014年用于繁殖的草鱼亲本(雌鱼、雄鱼各15尾)进行全基因组重测序，测序平台为Illumina Xten，测序深度为20×。
S2)获取S1)中描述亲本的子代体重数据，共采集到190个亲本对(部分亲本对没有采集到有效的子代表型数据)的1729尾子代表型数据(示例：表1)：
表1：1729尾子代体重数据(示例)
样本编号
母本
父本
体重(g)
A1681
F14
M12
109.16
B0039
F5
x9
119.01
A1693
F14
2
99.95
A2479
F8
M9
134.84
…
…
…
…
A1664
F14
M9
87.39
S3)使用BWA软件对30尾亲本的全基因组重测序数据进行比对(默认参数)，参考基因组为已发表的草鱼基因组。
S4)使用GATK工具(https://gatk.broadinstitute.org/)对得到的亲本比对结果进行SNP 检测，共检测到6,933,331个SNP位点(筛选指标为QUAL<30.0，SOR>10.0，QD<2.0，FS>200.0，InbreedingCoeff<-0.8)。
S5)从VCF文件中提取亲本SNP位点的基因型信息(例如提取到编号F14的雌鱼在CI01000000.1296位点基因型为GG、编号M12的雄鱼在其基因型为AG)，得到SNP位点 在亲本中的6933331×30阶分型矩阵(示例：表2)：
表2：30尾亲本的基因型(示例)
位点
F14
F15
…
M12
CI01000000.1296
GG
GG
AG
CI01000000.4999
TA
TT
…
TT
CI01000000.5103
GG
GG
…
GG
…
…
…
…
…
CI01000000.5253
GG
GG
…
GA
...
...
...
...
...
S6)接着使用下述原理和方法进行亲本基因型和子代表型的关联分析：
S6.1)子代体重异常值过滤：用R软件绘制的箱形图(如图1所示)，确定 [Q1-3IQR,Q3+3IQR]区间外的极端异常值14个，保留剩余1715个子代样本数据进行后续 分析。其中：Q1、Q3分别表示子代体重数据的下四分位和上四分位，IQR表示四分位距。
S6.2)依据子代的亲本信息，结合SNP位点在亲本中的分型矩阵，获取子代在各SNP位点的基因型组合信息(例如编号A1681的样本为F14与M12的子代，其在CI01000000.1296位点的组合基因型即为GG×AG)，依此得到SNP位点在子代中的6933331×1715阶分型矩 阵(示例：表3)：
表3：子代基因型组合(示例)
位点
A1681
B0039
…
A1693
CI01000000.1296
GG×AG
GG×GG
...
GG×GG
CI01000000.4999
TA×TT
TA×TT
…
TT×TT
CI01000000.5103
GG×GG
GG×GG
…
GG×GG
…
…
…
…
…
CI01000000.5253
GG×GA
GA×GG
…
GG×GG
...
...
...
...
...
结合各子代体重数据，获得每个SNP位点不同子代的可能基因型组合信息，以及对应 的子代群体分组表型数据。
S6.3)基于广义线性模型，在R环境下进行子代的组合基因型与其体重的关联分析，获 取各SNP位点的关联P值。其中组合基因型均无法获取有效信息(即“NN”)的位点标记为“ALLcontainN”(572个位点)，所有子代基因型组合相同的位点标记为“NULL”(32728 个位点)，使用剔除这两种标记类型位点后的数据进行后续分析。
S6.4)统计SNP位点的子代基因型组合信息，按不同的基因型组合对子代进行分组，获 取每组的子代体重数据并计算每组的体重均值，判断各组子代体重均值的大小关系，并使用 T-检验方法检验各组子代体重是否具有显著性差异，依据加性显性效应的原理与方法确定各 位点是否具有加性或显性效应，加性用“additivity”标记，显性用“dominance”标记。
以CI01000030.197128位点对上述各步骤产生的结果进行说明：
G C 1.07728069309273e-14 2 GC×GG|GG×GG 426|1289 27.5840610328639|23.1349495733126 1；15；15|14；15；174
2.458345e-13 additivity|dominance
F14:426；19:19,2:15,7:1,M12:48,M13:15,M14:26,M15:25,M19:25,M3:67,M30:21,M4:14,M5:39, M7:49,M9:57,x9:5；F14×19:19,F14×2:15,F14×7:1,F14×M12:48,F14×M13:15,F14×M14:26,F14 ×M15:25,F14×M19:25,F14×M3:67,F14×M30:21,F14×M4:14,F14×M5:39,F14×M7:49,F14×M 9:57,F14×x9:5|17:49,F10:33,F13:12,F15:206,F17:312,F28:43,F32:117,F34:200,F5:10,F8:140,x2: 29,x3:36,x4:32,x5:70；19:60,2:46,7:8,M12:133,M13:20,M14:113,M15:86,M19:79,M3:124,M30:7 6,M4:91,M5:131,M7:128,M9:162,x9:32；17×19:2,17×M12:4,17×M13:1,17×M14:6,17×M15:5,1 7×M19:3,17×M3:5,17×M30:1,17×M4:3,17×M5:5,17×M7:4,17×M9:7,17×x9:3,F10×19:2,F10 ×2:3,F10×M12:3,F10×M14:6,F10×M15:2,F10×M19:2,F10×M3:2,F10×M30:3,F10×M4:2,F10 ×M5:4,F10×M7:2,F10×M9:2,F13×2:1,F13×M14:1,F13×M15:3,F13×M19:1,F13×M3:1,F13× M30:1,F13×M7:2,F13×M9:1,F13×x9:1,F15×19:8,F15×2:5,F15×M12:33,F15×M13:7,F15×M1 4:14,F15×M15:8,F15×M19:12,F15×M3:16,F15×M30:11,F15×M4:15,F15×M5:23,F15×M7:21, F15×M9:24,F15×x9:9,F17×19:15,F17×2:7,F17×7:2,F17×M12:26,F17×M13:3,F17×M14:30,F1 7×M15:18,F17×M19:24,F17×M3:32,F17×M30:16,F17×M4:29,F17×M5:41,F17×M7:37,F17× M9:25,F17×x9:7,F28×19:4,F28×2:2,F28×7:1,F28×M12:4,F28×M13:1,F28×M14:5,F28×M15: 6,F28×M30:1,F28×M4:3,F28×M5:5,F28×M7:3,F28×M9:8,F32×19:1,F32×2:12,F32×7:2,F32× M12:17,F32×M13:2,F32×M14:9,F32×M15:5,F32×M19:6,F32×M3:3,F32×M30:7,F32×M4:4,F 32×M5:6,F32×M7:10,F32×M9:29,F32×x9:4,F34×19:11,F34×2:5,F34×7:1,F34×M12:12,F34× M13:3,F34×M14:14,F34×M15:11,F34×M19:9,F34×M3:40,F34×M30:8,F34×M4:14,F34×M5:2 1,F34×M7:23,F34×M9:27,F34×x9:1,F5×19:1,F5×2:1,F5×M12:2,F5×M19:1,F5×M5:2,F5×M7: 1,F5×M9:2,F8×19:6,F8×2:5,F8×M12:12,F8×M14:13,F8×M15:15,F8×M19:14,F8×M3:14,F8× M30:14,F8×M4:8,F8×M5:11,F8×M7:9,F8×M9:14,F8×x9:5,x2×19:3,x2×M12:4,x2×M14:2,x2 ×M15:3,x2×M19:1,x2×M3:2,x2×M4:1,x2×M7:3,x2×M9:8,x2×x9:2,x3×7:1,x3×M12:4,x3×M1 3:1,x3×M14:3,x3×M15:3,x3×M19:2,x3×M3:3,x3×M30:3,x3×M4:6,x3×M5:4,x3×M7:2,x3×M 9:4,x4×19:4,x4×2:1,x4×M12:5,x4×M13:1,x4×M14:1,x4×M15:3,x4×M19:2,x4×M3:1,x4×M30 :2,x4×M4:1,x4×M5:3,x4×M7:6,x4×M9:2,x5×19:3,x5×2:4,x5×7:1,x5×M12:7,x5×M13:1,x5× M14:9,x5×M15:4,x5×M19:2,x5×M3:5,x5×M30:9,x5×M4:5,x5×M5:6,x5×M7:5,x5×M9:9 intron_variant MODIFIERGENE_CI01000030_00200402_00425242.ExON
其中“CI01000030.197128”表示位点信息；
“G”表示CI01000030.197128位点在参考基因组上的碱基为G；
“C”表示CI01000030.197128位点的碱基突变为C；
“1.07728069309273e-14”表示由GLM模型进行统计检验计算得到的P值；
“2”表示此位点子代基因型组合种数为2；
“GC×GG|GG×GG”表示子代基因型组合有GC×GG、GG×GG两种，并且已经按照 其对应子代体重均值由大到小排序；
“426|1289”表示基因型组合为GC×GG的子代有426个，基因型组合为GG×GG的子代有1289个；
“27.5840610328639|23.1349495733126”表示基因型组合为GC×GG的子代体重均值 为27.5840610328639，基因型组合为GG×GG的子代体重均值为23.1349495733126；
“1；15；15|14；15；174”表示每种基因型组合下的母本、父本、亲本对个数，即基因型组 合为GC×GG的子代母本有1个，父本有15个，亲本对有15个，基因型组合为GG×GG 的子代母本有14个，父本有15个，亲本对有174个；
“2.458345e-13”表示为分析加性显性效应获取的T-检验P值。此处的P值远远小于我 们常设置的0.05或者0.01，故认为两种基因型的子代体重具有显著差异。而两种基因型组 合分别为GC×GG和GG×GG，此位点原始碱基为G，变异碱基为C，其满足W(0/0×0/1) >W(0/0×0/0)特征，根据加性显性效应原理故得到此位点既具有加性效应又具有显性效 应的特征；
“additivity|dominance”表示此位点的基因型组合即具有加性效应又具有显性效应；
“F14:426；19:19,2:15,7:1,M12:48,M13:15,M14:26,M15:25,M19:25,M3:67,M30:21,M4:14, M5:39,M7:49,M9:57,x9:5；F14×19:19,F14×2:15,F14×7:1,F14×M12:48,F14×M13:15,F14×M14: 26,F14×M15:25,F14×M19:25,F14×M3:67,F14×M30:21,F14×M4:14,F14×M5:39,F14×M7:49, F14×M9:57,F14×x9:5|...”表示各基因型组合下的具体母本、父本、亲本对ID、分别有几个 子代。
“|”分割各基因型组合，如上即是GC×GG|GG×GG的对应信息；“；”分割母本、父本、亲本对，即“母本；父本；亲本对”；“,”分隔单个母本/父本/亲本对，"19:19,2:15,7:1,M12:48,...；" 即表示基因型组合为GC×GG子代群体其15个父本分别为19、2、7、M12、...，子代个数 对应分别为19、15、1、48、...。
S6.5)使用PopLDdecay软件检测全基因组SNP位点的连锁不平衡衰减效应，计算强相 关标记的平均距离DLD，其值低于最低值10bp，故DLD＝10。
S6.6)使用snpEFF软件分析各位点在基因组上的位置，获得位点的基因组注释信息， 对所有SNP按照所属基因进行分组。
S6.7)使用以下步骤筛选显著关联基因：
①设定距离阈值K＝2000，获取基因上以及距离基因上下游2000bp以内的SNP位点；
②筛选出S6.4)中标记出的具有加性或显性效应的位点；
③确定关联P值的阈值p0＝10^(-12)，依据S6.3)中计算出的各位点关联P值，筛选出以基因为分组单元的各组区域内连续2个及以上位点关联P值小于p0的区域块；
④计算上述区域块上起始位点与终止位点之间的距离D(即多标记间的距离)，距离阈 值d＝max{350，DLD}＝350，筛选出D大于350的区域块，将这些区域块确定为候选关联标 记；
⑤计算候选关联标记上各位点P值的几何均值Pm，确定几何均值的阈值p1＝1e-17， 从候选关联标记中筛选出Pm小于1e-17的标记，获得关联标记。根据这些关联标记注释到 的基因信息，获得28个基因，部分展示如下：
表4：基于新型关联分析方法的高支持组别关联基因列表(示例)
本实施例使用的群体是长期处于饥饿状态的草鱼子代群体，其生长的快慢取决于个体获 取食物的能力。关联结果的基因主要和个体的视觉、听觉、小脑的发育有关，这些都是保障 个体能更有效找食与摄食的先决条件。另外，过强的氧化应激和免疫应答都会抑制个体的生 长，在表4中也出现了上述相关基因。综上所述，本方法关联到的基因是符合群体的生存环 境特征的。本方法适用于二倍体且子代数目较大的物种，比如鱼类与植物等，目标性状可以 是数量性状，也可以是质量性状。
- 上一篇：一种医用注射器针头装配设备
- 下一篇：一种同源重组修复基因变异的解读方法