分子标记是遗传学研究的重要内容。单核苷酸多态性标记因其数量多、分布广、易于准确鉴定等特点，被广泛地应用于人类和动植物各项遗传学研究，包括遗传图谱构建、连锁分析、关联分析、基因组预测等。Taqman探针、Sanger测序、Sequenom质谱检测、Fluidigm微流体芯片等技术提供了中低通量下准确的单核苷酸多态标记分型信息，但是其通量难以达到全基因组水平，所以并不适合基因组水平的研究。如何高效准确地对全基因组范围内的SNPs标记进行分型是遗传学重要的研究内容，目前主流的全基因组标记获取方法包括SNP分型芯片、简化基因组测序以及全基因组低深度重测序等。
另一方面，在获得了全基因组标记信息后，如何使用全基因组标记进行基因组预测也是遗传学研究的重要内容。对于数量性状来说，由于通常由微效多基因控制，使用少量的分子标记往往只能解释部分遗传方差，有必要使用全基因组标记进行个体表型预测。然而，一个个体的全基因组能够找到数以百万计的标记，而通常情况下的研究样本群体中具有某特定表型的个体数仅有千数量级[1]。传统的预测方法使用这种高维小样本的数据，容易引发“维度灾难” (dimensionality curse)和过拟合(over fitting)问题。为了解决这个问题，目前主流的基因组预测方法包括：正则化线性回归、基因组选择方法和机器学习方法这几类。本文分别从全基因组基因分型方法与基因组预测方法两个重要方面展开阐述。
1 全基因组基因型分型策略 1.1 全基因组SNP芯片基因分型方法基因组SNP芯片是第一种广泛应用的基因组水平遗传标记基因分型方法。目前主流的基因组SNP芯片技术为Illumina Infinium技术和Axiom (原Affymetrix公司，被Thermo公司收购)原位光刻技术。Illumina Infinium技术采用全基因组扩增的方式，不需要进行PCR反应；该方法采用50 mer寡核苷酸探针退火，利用特异荧光基团判定基因型。Axiom基因分型基于原位光刻技术，需要对DNA定量后进行PCR扩增；该方法采用30 mer的寡核苷酸探针退火，并同样利用特异荧光基团判定基因型。SNP基因分型芯片价格适中、准确度高、重复性好、试验流程标准化程度高，是目前最常用的全基因组SNP分型方法[2-3]。
近年来，研究人员利用SNP芯片在人类遗传、动植物遗传育种等领域开展了大量研究。人类单倍型图谱计划HapMap提供了大量的人类基因组变异信息[4]，根据其设计出的芯片为人类疾病的研究提供了极大的便利[5]。猪、鸡、牛、羊等重要农业动物都已经有商业化的高密度SNP芯片(表 1)，并已广泛应用于重要功能基因挖掘、人工选择信号鉴定、基因组选择育种等领域。Qiao等[6]利用山羊SNP芯片进行全基因组关联分析(genome-wide association study，GWAS)，定位到了多个与羊绒细度性状相关的候选基因；Gu等[7]利用60K鸡SNP芯片进行全基因组关联分析，定位了多个影响体重的主效QTLs区间；Zhu等[8]通过绵羊50K芯片对3个绵羊品种进行选择信号分析，检测到X染色体上多个近期的选择信号；Ostersen等[9]通过猪60K芯片对1 911头丹麦杜洛克猪进行基因分型，通过基因组最佳线性无偏估计(genomic best linear unbiased prediction, GBLUP)等方法证明，对于纯种猪的基因组选择来说，去除父母均值的估计育种值适合作为估计中的因变量。
|
表 1 主要畜禽商业化SNP芯片[11] Table 1 Commercial SNP genotyping chips for major livestock and poultry species |
随着研究的扩展，商业化SNP芯片的不足也逐渐显现出来。首先，芯片所含的标记数量有限，难以满足所有类型的研究需求。例如，在QTL的精细定位研究中，研究者们通常需要较多的标记来捕获特定区段所有的重组事件或单倍型多态，而一般标记数为几万到几十万的芯片很难达到QTL精细定位所需的标记密度。其次，SNP芯片只能检测已知突变，无法检测新生突变。最后，商业化SNP芯片一般根据某些知名品种设计，与大多数研究的地方品种或闭锁群体的遗传距离较大，造成芯片的部分标记位点在特定群体中失效[10]。为了解决这个问题，需要对特定品种特定群体进行芯片定制，这在一定程度上提高了芯片设计和制作的成本。这些不足以及测序方法的进步使得SNP芯片的应用受到限制，并逐步被测序技术所取代。
1.2 基于测序的标记开发策略近年来，高通量测序技术快速发展，测序通量显著增加，成本快速下降，使得研究人员能够通过测序的手段更加直接和全面地获得个体或群体的基因组序列信息。基因组测序理论上可以获得其基因组上几乎所有的变异，为遗传研究、动植物育种、疾病研究等提供最全面的材料。然而，在群体水平的研究中，对成百上千的样本进行全基因组深度测序需要极高的成本。于是，研究人员开发了简化基因组测序(reduced-representation genome sequencing，RRGS)以及全基因组低深度重测序的方法，在获得足够标记数目的情况下，大大降低测序成本。与SNP芯片相比，简化基因组测序和低深度重测序能够获得更高的标记密度，并克服样本特异性的缺陷，目前，已广泛应用于动植物群体遗传学研究。
1.2.1 简化基因组测序简化基因组测序首先通过酶切的方式打断基因组，之后富集特定酶切位点附近一定长度的基因组片段进行测序。酶切位点可以出现在基因组的任意位置，由于同一物种的不同个体或者近缘物种间的酶切位点位置相对保守，共享相同的酶切位点为酶切简化测序提供了先决条件[11]。RRGS能在一次试验中获得比中通量SNP芯片更高密度的、覆盖全基因组的遗传标记。由于被测序的是基因组的一部分，因此，在相同的测序通量和测序深度的情况下，能够测得更多的个体，这也为群体遗传学中大量样本的基因分型提供了可能。利用限制性内切酶来鉴定基因组变异渊源已久[12-13]，但早期主要关注的是酶切位点本身携带的突变。简化基因组测序则依托二代测序技术，获取酶切位点附近序列的标记信息[14-15]，目前，常用的RRGS方法流程如图 1所示。
|
深蓝色片段代表基因组插入片段，其余颜色代表测序接头的不同组成部分 The dark blue fragments represent the genome insertion and the remaining colors represent the different components of the sequencing adapter 图 1 目前常用的RRGS方法的建库流程示意图 Fig. 1 Experimental flow chart of several RRGS methods |
简化基因组测序的流程灵活性使其在短短几年间出现了很多不同的试验方案，从试验设计到应用，方法上的每一步调整都会直接影响试验的步骤和数据的产出。由于试验方案的调整，产生多个RRGS的命名方法。例如RAD-seq (restriction site-associated DNA sequencing)以方法导向命名[16]，而GBS (genotyping by sequencing)，其名称表达了该方法的目的性[17]，这实际上造成了RRGS方法命名一定程度的混乱。迄今为止，RRGS已经包含10余种不同命名方法，但目前广泛应用的RRGS技术主要是RAD-seq和GBS方法及其改良版本。各种简化基因组测序的通用之处是通过一个或者一组限制性内切酶进行酶切基因组，然后，与特定的测序接头或者双链寡聚核苷酸序列连接，最后，进行高通量测序。简化基因组测序不同方法的区别，主要体现在内切酶选择、酶切方式、接头连接方法、标签序列设计、片段大小筛选流程、测序数据类型、数据分析方法等环节。
RRGS能在一次试验中以更低的成本获得比中通量SNP芯片更高密度的、覆盖全基因组的遗传标记。在奶牛、猪、羊、鸡、鸭、鱼等重要畜禽、水产物种中，RRGS方法已有大量的应用[18]。美国农业部下属的玉米、小麦遗传改良实验室都将GBS技术作为重点研究和推广的技术，用于玉米、小麦的多样化研究以及全基因组关联分析研究和基因组选择[19]。对于特定物种来说，研发出一套成熟的RRGS流程，就意味着低成本、高通量、高准确性基因分型体系的建立，这是商业化SNP芯片所无法比拟的[20]。在我国，成熟的猪GBS流程已经应用在广东温氏公司大规模的杜洛克猪基因组选择育种研究中，并取得了显著的遗传进展[21]。
由于技术上的原因，简化基因组测序特有的误差和偏好也值得研究者注意。首先，当同一限制性酶切位点在不同的个体间出现多态时，就有可能造成部分个体在此位点无法正常酶切，进而导致某一种等位基因的丢失，最终表现为某一杂合子位点被误判为纯合子。如果某一点的酶切出现问题，其临近位点的酶切片段长度也会同时受到干扰，如果其长度超过片段筛选步骤中的预期长度，这一临近位点将可能被错误地丢弃[22]。其次，对于很多二代测序的基因组文库，基因组随机打断后不同DNA片段的起止位置几乎不可能完全相同，这样就可以通过信息学的手段过滤PCR所产生的重复DNA片段。然而，这种方式不适合大多数简化基因组测序文库。因为简化基因组方法产生的片段均位于酶切位点处，其首末位置理论上是完全相同的，使得其表现类似于PCR重复而无法区分。目前，已开发出一些方法来降低PCR的影响，比如利用简并碱基降低PCR干扰[23-24]，甚至在流程中去除PCR扩增步骤(例如ezRAD技术)[25]。
1.2.2 全基因组低深度重测序全基因组低深度重测序是继RRGS之后新一代低测序成本的标记开发方法。该策略首先对群体中所有个体进行全基因组低深度测序和变异查找，之后根据SNP位点间的连锁不平衡(linkage disequilibrium，LD)对缺失基因型进行推断，最终获取大规模样本全基因组水平的高密度遗传标记[26]，其大体流程如图 2所示。
|
蓝色方框和白色方框分别代表已知基因型和缺失基因型 The blue and white blocks represent the known and missing genotypes, respectively 图 2 基于低深度重测序数据的基因型填充 Fig. 2 Imputation based on low-coverage resequencing data |
对于单个个体低深度测序而言，有可能会出现将杂合子错判为纯合子的结果[27]。所以，低深度测序最为核心的步骤就是群体水平的基因型填充(imputation)。通过已知的单倍型参考数据或者大量样本中样本间的共享单倍型信息，使用统计学方法对样本中缺失基因型进行推断。目前，常见的基因型填充策略大多是基于参考单倍型，大体分为两个步骤：第一步对目标样本进行预定相(pre-phasing)，从而得到区域内的单倍型信息；第二步根据参考数据集的单倍型对缺失基因型进行填充[28]。目前，对人类基因组的研究通常使用SHAPEIT2获取等位基因的相位信息[29]，之后通过BEAGLE[30]、IMPUTE2[31]等工具结合参考数据集进行填充。不依赖于参考单倍型的填充而借助大量个体的测序数据中片段的共享对群体祖先单倍型进行推断，之后再根据推断出的单倍型对每个个体的缺失数据进行填充。Davies等[2]开发的STITCH (sequencing to imputation through constructing haplotypes)软件通过使用隐马尔可夫模型以及EM算法来对研究群体的祖先单倍型进行估计，之后根据祖先单倍型对低深度重测序数据中的缺失基因型进行填充，在一定程度上克服了非人类基因组研究中缺少高质量参考单倍型数据集的难题。2010年，Huang等[32]通过对517个水稻地方品种平均约1×/个体的低深度测序，利用数据填充的方法获得了较高准确性的SNP数据集，并以此构建了高密度的水稻基因组单体型(haplotype)图谱，利用填充后的高密度变异图谱，随后对籼稻亚种的14种农艺性状进行了全基因组关联分析，鉴定到的基因座平均解释约36%的表型变异。2018年，Liu等[33]完成了一项大型全基因组超低深度数据的关联分析，数据来自14万余中国人的无创产前基因检测(non-invasive prenatal testing，NIPT)，每一例样本只有(0.06~0.1)×的测序深度。试验采用自主开发的BaseVar软件替代GATK软件进行基因型鉴定，并使用STITCH软件进行基因型填充，运算速度和准确性均有非常出色的表现。在模式动物小鼠中，Nicod等[34-35]使用0.15×的超低深度测序对1 887只远缘杂交系小鼠进行了基因型分型和全基因组关联分析。针对92个性状定位到了156个独立的关联标记。在猪、鸡等重要畜禽中，基于低深度数据分析的各类研究结果也初见报道[36]，这些研究也为人们在畜禽物种中进行大样本低深度的基因分型提供了新的思路。
2 基因组预测方法随着SNP分型技术的快速发展，分型成本逐渐降低，使得利用全基因组SNP分型数据进行基因组预测成为可能。全基因组分型数据提供了更全面的遗传信息，理论上能解释更多的遗传方差。然而，全基因组数据标记数目远远大于分型的样本数，标记之间由于连锁不平衡或上位效应往往具有共线性，这导致了简单的回归模型不能很好拟合。因此，大量SNP分型数据的获得也对基因组预测方法的开发与完善提出了很大的挑战。目前，基于SNP分型数据的基因组预测的方法主要有正则化线性回归、基因组选择、机器学习等方法。
2.1 正则化线性回归线性回归模型为y=Xβ+ε，其中，y为n个个体的表型值向量，β是f维标记效应向量，X是所有标记所组成的已知n×f矩阵，ε是n维随机误差向量。正则化线性回归的模型求解基于最小二乘法，即最小化每个样本的估计值与真实值之间误差的平方，最小二乘法为每个标记分配一个较大的回归系数，使得对标记的取值非常敏感。当数据的样本量远远小于标记数或数据中存在多重共线性时，最小二乘法往往找不到较好的模型对数据进行拟合。正则化线性回归在代价函数后面加上正则化项，用于防止过拟合，提高模型的泛化能力。两种经典的正则化回归模型分别为嵌套法(the least absolute shrinkage and selection operator, LASSO)和岭回归(ridge regression, RR)。
2.1.1 嵌套法1996年，Tibshirani[37]首次提出了LASSO。该方法采用L1正则化，即对参数的绝对值进行压缩，该方法对标记效应的估计如下：
| $ \begin{array}{l} \hat \beta = \\ \mathop {{\rm{argmin}}}\limits_\beta \left\{ {\sum\nolimits_{i = 1}^n {{{\left( {{y_i} - {\beta _0} - \sum\nolimits_{j = 1}^f {{\beta _j}{x_{ij}}} } \right)}^2}} + \lambda \sum\nolimits_j {\left| {{\beta _j}} \right|} } \right\} \end{array} $ | (1) |
其中，yi表示样本i真实表型值，xij为样本i的第j个特征取值，∑j=1fβjxij表示样本i的表型估计值，λ表示正则化参数，λ越大表示正则化力度越大，
LASSO方法使得所有标记的系数绝对值之和小于一个常数，因此，当标记的效应较小时，回归系数被降为0，从而实现降维进而增加了模型的可解释性。LASSO回归方法适用于大幅压缩变量的情况，目前该方法在众多领域都得到了应用[38-39]。最近，Arbet等[40]使用LASSO进行GWAS，并针对正则化系数推断问题提出了几种优化方案，其中，LASSO-PL和LASSO-AL在应对高维数据时依然表现优秀。
2.1.2 岭回归早在20世纪70年代，RR就被用于积分方程，使得RR广为人知，因此，RR又名Tikhonov regularization。岭回归采用L2正则化，即对参数的平方进行压缩，该方法对标记效应的估计如下：
| $ \begin{array}{l} \hat \beta = \\ \mathop {argmin}\limits_\beta \left\{ {\sum\nolimits_{i = 1}^n {{{\left( {{y_i} - {\beta _0} - \sum\nolimits_{j = 1}^f {{\beta _j}{x_{ij}}} } \right)}^2}} + \lambda \sum\nolimits_j {\beta _j^2} } \right\} \end{array} $ | (2) |
其中，yi表示样本i真实表型值，xij为样本i的第j个特征取值，∑j=1fβjxij表示样本i的表型估计值，λ表示正则化参数，λ越大表示正则化力度越大，
RR方法中，参数的压缩程度随着标记效应进行调整。当标记的效应β非常小时，压缩的力度也会相应变小，该方法不会将原本效应非0的标记压缩为0。De Vlaming和Groenen[41]指出，在面临较大规模的数据时，相对于重复对单点进行回归的方法，RR将显著提高预测的准确度。
2.2 基因组选择(genomic selection, GS)在动物育种领域中，需要一种准确的方法对优良的种牛、种猪等个体进行早期选择，进而实现高效的品种选育。传统回归方法没有将固定效应与随机效应进行剖分，无法直接得到个体可遗传部分的育种值。另外，传统的回归方法未对标记的效应分布进行假设，以致在进行效应求解时的计算量巨大。GS利用覆盖全基因组高密度SNPs信息，通过标记效应的求解和加和或基因型关系矩阵的求解，得到个体基因组的估计育种值(genomic estimated breeding value, GEBV)，从而达到对样本个体进行打分排序和早期选择的目的。
2.2.1 基因组最佳线性无偏估计2008年，Vanraden[42]提出基因组最佳线性无偏估计，GBLUP使用SNP构建基因组关系矩阵(genomic relationship matrix, GRM，G)，然后，采用混合线性模型估计GEBV。该方法的模型：
| $ y = X\beta + Z\alpha + \varepsilon $ | (3) |
其中，y为n个个体的表型值向量；β是f维“固定”效应向量；X是所有固定效应所组成的已知n×f矩阵；α表示q维随机效应向量，α~N(0, Iσα2)，σα2服从逆卡方分布，即：σα2~χ-2(vα, Sδ2)；Z为随机效应的n×q矩阵；ε是n维随机误差向量。
Vanraden[42]提出了G矩阵的构建方法。首先，对SNP的基因型进行编码，编码的方式可以是-1、0、1，即纯合基因型编码为-1和1，杂合基因型编码为0；或0、1、2，即低频基因型编码为2，杂合基因型编码为1，高频基因型编码为0。为了进行中心化，将每个基因型编码后的值减去相应的处于哈代-温伯格平衡(hardy-weinberg equilibrium, HWE)时的均值2×(pj-0.5), pj为j位点的最小等位基因频率，得到Z矩阵，并计算得到G矩阵：
| $ G = ZZ'/2\Sigma {p_i}\left( {1 - {p_i}} \right) $ | (4) |
随后，根据混合模型方程组(mixed model equation，MME)对方程进行求解：
| $ \left[ {\begin{array}{*{20}{c}} {\hat \beta }\\ {\hat \alpha } \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {X{X^\prime }}&{{X^\prime }Z}\\ {{Z^\prime }X}&{{Z^\prime }Z + \lambda {G^{ - 1}}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{X^\prime }y}\\ {{Z^\prime }y} \end{array}} \right] $ | (5) |
其中，G-1是基因组关系矩阵G的逆矩阵，λ=σe2/σα2。
迄今为止，GBLUP方法已经在荷斯坦牛育种中取得了较大的成效。自2018年GS技术引入美国以来，公牛和母牛的世代间隔显著降低，低遗传力性状的选择强度显著提高，生育能力、寿命、乳房健康等性状也得到显著改良[43]。GBLUP方法使用混合线性模型，更准确的得到可遗传方差组分。2019年，Teng等[44]将GBLUP用于杂交鸡的20个生长性状的预测，结果表明，对于大多数性状GBLUP的预测准确性高于传统的BLUP方法。相较于传统的回归方法，GBLUP方法使用混合线性模型，更准确的得到可遗传方差组分。然而，大多数情况下，基因组中大多数标记具有小效应，少部分标记具有大效应，这与GS的假设不一致[45]。因此，GBLUP方法的预测准确度还有待进一步提升。
2.2.2 一步法基因组最佳线性无偏估计(single-step GBLUP, ssGBLUP)传统的BLUP方法使用系谱信息与表型信息进行预测。为了将系谱信息有效的利用起来，以达到更好的估计效果，2010年，Aguilar等[46]在GBLUP的基础上提出了ssGBLUP。ssGBLUP使用一个联合系谱信息和基因型信息的矩阵H代替传统BLUP方法中的亲缘关系矩阵A。
| $ H = \omega A + (1 - \omega )G $ | (6) |
其中，ω表示系谱矩阵的尺度，A表示系谱矩阵，G表示基因型关系矩阵。
| $ \begin{array}{l} A = \left[ {\begin{array}{*{20}{l}} {{A_{11}}}&{{A_{12}}}\\ {{A_{21}}}&{{A_{22}}} \end{array}} \right]\\ H = \left[ {\begin{array}{*{20}{c}} {{A_{11}}}&{{A_{12}}}\\ {{A_{21}}}&G \end{array}} \right] = A + \left[ {\begin{array}{*{20}{c}} 0&0\\ 0&{G - {A_{22}}} \end{array}} \right]\\ {H^{ - 1}} = {A^{ - 1}} + \left[ {\begin{array}{*{20}{c}} 0&0\\ 0&{{G^{ - 1}} - A_{22}^{ - 1}} \end{array}} \right] \end{array} $ | (7) |
其中，H-1是H的逆矩阵。A的下标1和2分别代表没有进行基因分型的群体和进行基因分型的群体，A-1是A的逆矩阵。
Fragomeni等[47]使用加权和不加权的SNP进行GBLUP和ssGBLUP预测。结果显示，加权的SNP对GBLUP有较大的影响，而ssGBLUP方法不受到SNP是否加权的影响，相较于GBLUP可以解释更多的遗传方差。2018年，Teissier等[48]使用ssGBLUP对法国山羊产奶蛋白含量性状进行基因组育种值估计，结果表明，基于ssGBLUP得到的GEBV的准确性比基于BLUP得到的准确性高5%~7%。
2.2.3 贝叶斯(Bayes)估计GBLUP及其衍生方法假设所有的标记解释相等的方差，这种假设很多情况下与事实不符, 自Meuwissen等[49]第一次提出BayesA、BayesB以来，目前，应用于育种值估计的贝叶斯方法按照先验假设的不同共分为BayesA、BayesB、BayesCπ、BayesDπ等多种方法。该方法的通用模型：
| $ y = X\beta + \sum\nolimits_{j = 1}^f {{z_j}{g_j} + \varepsilon } $ | (8) |
其中，y为n个个体的表型值向量，zj是n个个体在第j个SNP的基因型向量，对应于SNP的k种基因型，Zki取值为-1、0、1。
BayesA的先验分布：SNP效应gk服从均值为0，方差为σgk2的正态分布; SNP效应方差服从自由度为ν，尺度参数为S的逆卡方分布；残差效应方差σe2服从自由度为-2，尺度参数为0的逆卡方分布。使用MCMC算法，对各个变量进行Gibbs抽样，获得变量的估计值[49]。BayesB的先验假设：概率π，σgk2=0；概率1-π，p(σgk2)=x-2(ν, S), π人为给定[49]。其他先验假设与BayesA相同。在对变量进行估计时，除σgk2使用Metropolis-Hastings算法，其他与BayesA相同。BayesCπ和BayesDπ假设π是未知参数，服从Uniform(0, 1)；BayesCπ将概率1-π的非零效应的SNP的方差设为一个共同方差，BayesDπ假设尺度参数S为未知参数，服从Gamma(1, 1)[50]。
总体来说，Bayes方法估计GEBVs的准确性优于BLUP方法，在具有较大QTL的性状上优势更为明显[45]。然而由于其计算效率相对较低，过程较为复杂，目前应用相对较少。
2.3 机器学习(machine learning)无论是传统的线性回归方法还是基因组选择方法，更多关注于单个标记与表型之间的关联，而对于复杂性状来说，标记之间的上位效应也是不可忽视的一部分。另外，随着测序技术的发展，标记数目大量增加，给传统方法在计算上也带来挑战。机器学习通过采用已有的“学习算法” (learning algorithm)和经验数据，基于已有数据训练出模型，进而对新的数据进行预测。基因组预测致力于构建基因型和数量性状表型值之间的关联，属于监督学习中的回归问题。目前，几个主流方法有支持向量回归(support vector regression, SVR)、随机森林回归(random forest regression, RFR)、人工神经网络(artificial neural network, ANN)等。
2.3.1 支持向量回归1996年，支持向量机(support vector machine, SVM)的回归版本SVR被提出。SVR改变了SVM结构经验最小化的原则，以置信范围最小化作为优化问题的目标。SVR忽略小于常数ε的损失，得到较为稀疏的支持向量，构建出这部分样本的损失函数并进行最小化求解。SVR的模型：
| $ f\left( x \right) = {w^T}x + b $ | (9) |
其中，训练样本为D={(x1, y1), (x2, y2), …(xn, yn)}，x是q维输入向量，w、b是未知参数。
引入ε函数作为损失函数，松弛因子ξi，
| $ \begin{array}{l} \min \frac{1}{2}{\left\| \omega \right\|^2} + C\sum\limits_{i = 1}^m {\left( {{\xi _i} + {{\hat \xi }_i}} \right)} \\ {\rm{s}}.\;{\rm{t}}.\;f\left( {{x_i}} \right) - {y_i} \le \varepsilon + {\xi _i},\\ {y_i} - f\left( {{x_i}} \right) \le \varepsilon + {{\hat \xi }_i}\\ {\xi _i} \ge 0,{{\hat \xi }_i} \ge 0,i = 1,2 \cdots ,{\rm{m}}{\rm{.}} \end{array} $ | (10) |
其中，C是正则化参数，ε表示被容忍的f(x)与y之间的偏差；ξi，
| $ f\left( x \right) = \sum\nolimits_{i = 1}^m {\left( {{\alpha _i} - {{\hat \alpha }_i}} \right)\kappa \left( {x,{x_i}} \right) + b} $ | (11) |
其中，αi，
SVR非常适合高维小样本数据的分析，能很好的解决小样本、非线性、维度灾难的问题。2019年，Chen等[51]提出PredPSI-SVR，该方法使用SVR稳健的对外显子剪接进行预测。
2.3.2 随机森林回归1995年，Ho[52]在决策树的基础上提出了随机森林。RFR是一种集成学习(ensemble learning)方法，通过组合多个CART回归决策树，使得结果有较高的精确度和泛化能力。训练一个RFR模型首先对原始数据集(样本含量为n)有放回的进行采样，每次采集b个样本，得到k个数据集。没有采集到的数据称为袋外数据(out of bag, OOB)，用来检测模型的泛化能力。之后，使用CART回归树进行训练，每棵树选取p个特征进行训练，p < < f，f为总的特征数目。CART回归决策树采用Gain_σ作为分裂的评价指标，Gain_σ越小，说明划分后的左右两个子集之间的差异越小。对于样本集S，选择所有属性中划分后Gain_σ的最小值，作为S的最优划分方案。最后使用OOB对回归效果进行评价，预测值
相比SVR方法，随机森林具有较好的结果可视性，输出的随机森林模型可以给出每个标记的重要程度。2018年，Brennan等[53]分别使用RFR和全基因组关联分析方法对基因型和表型进行关联，作者认为RFR和GWAS方法可以识别到不同类型的遗传变异，如GWAS倾向于识别具有大效应的单个位点，而RFR方法可以识别到具有复杂遗传互作特征的位点，如上位效应。因此，二者共同识别到的遗传变异与表型有强关联效应。
2.3.3 人工神经网络早在1943年，Mcculloch和Pitts[54]就提出了神经网络的计算模型。ANN构建的想法来源于人类大脑分层次学习的架构。相对于SVR等“浅层学习”方法，ANN可对输入层的数据逐层处理，从而将特征逐层映射到新的特征空间。
目前，已经出现的ANN有数百种之多，其中，3种常见的神经网络架构为前馈神经网络(feedforward neural network, FNN)、卷积神经网络(convolutional neural network, CNN)和循环神经网络(recurrent neural network, RNN)[55]。
2018年，Waldmann[1]提出了用于基因组预测的前馈神经网络，该方法的模型：
| $ \hat y = \sigma \left( {X{W_1} + b} \right){W_2} $ | (12) |
其中，
模型的优化目标为：
| $ \begin{array}{*{20}{c}} {{\rm{minimize}}\left\{ {J\left( \theta \right) = } \right.}\\ {\frac{1}{{2n}}\sum\nolimits_{i = 1}^n {\left\| {{y_i} - {{\hat y}_i}} \right\|_2^2 + {\lambda _1}\left\| {{W_1}} \right\|_2^2 + } }\\ {\left. {{\lambda _2}\left\| {{W_2}} \right\|_2^2 + {\lambda _3}\left\| b \right\|_2^2} \right\}} \end{array} $ | (13) |
其中，θ表示参数的集合，λ1、λ2、λ3分别为对应于W1、W2、b的正则化参数。
在进行模型的训练时，首先，对W1、W2、b进行初始化，输入数据与现有的参数矩阵沿着神经网络逐层向前传播，到达输出层得到代价函数。代价函数沿着该网络进行反向传播，经过数次迭代获得最优参数。使用单隐藏层近贝叶斯神经网络对猪的5个匿名性状进行基因组预测，结果显示，该方法的准确性优于GBLUP和BayesLasso[56]。2018年，Ma等[56]提出使用卷积神经网络进行深度基因组选择(DeepGS)，作者将该方法与RRBLUP、GBLUP进行对比，结果显示，两种方法在高质量个体的捕获上具有互补性。
浅层的学习方法需要对问题本身具有一定的了解，即需要人为的先验知识，而深度学习方法通过一层层的累加自动实现对问题的认知。在面对极大数据的处理时，深度学习方法具有较大的优势。在进行深度学习的过程中，为了避免过拟合的情况，要适当把握网络的复杂度，适当调节正则化的方法[57]。
2.4 基因组预测的评估方法进行基因组预测后，需进行预测效果的评估，判断方法的有效性。常用的评估方法：1)相关性，两种相关性的评价分别为真实值y与估计值
近些年全基因组遗传标记分型方法的进步与发展极大地推进了群体遗传学的研究。总体而言, 本文介绍的3种全基因组基因分型的技术各有优缺点。首先，3种方法获取的标记密度存在数量级上的差别。SNP分型芯片和简化基因组测序能够得到全基因组数万到数百万的标记信息，而低深度重测序的方法理论上能得到全基因组的所有变异。其次，就每个个体的分型成本而言，SNP芯片最高，但流程标准化程度高；全基因组低深度重测序的费用最低，但需要以大样本量和高计算量为基础。从目前发展趋势和方法的灵活性来讲, 基于测序相关技术的前景要优于芯片技术，尤其是低深度重测序策略已被人们广泛关注和应用。
但是，在选择不同方法时需要考虑在标记密度和总的分型成本上达到平衡。例如，遗传图谱的构建可以采用SNP芯片或者简化基因组测序方法；而在功能基因精细定位和致因突变鉴定方面则低深度重测序的方法更佳。此外，样本量越大、群体杂合度越低、样本间的亲缘关系越近则越适合使用简化基因组测序和低深度重测序的方式，这种模式也非常适合大型育种群体基因组选择的相关研究。
本文总结的3种基因组预测方法也各有特点。Lasso将不重要的标记的效应降为0，从而实现降维的效果，而RR对于小效应的标记压缩的程度也相应降低，不能实现特征的选择。不同的GS方法对标记进行了不同的先验假设，其中BayesB、BayesCπ、BayesDπ均假设概率π的标记效应为0，从而实现降维。ML方法中，RFR使用所有标记进行训练，未进行降维操作。
正则化线性回归方法对标记的效应进行正则化，从而避免模型过拟合。GS方法将固定效应与随机效应分开，直接得到个体的育种值，GBLUP方法由于假设所有的标记具有相同的方差，大大加快了运算速度，由于能很好地利用全基因组标记信息获得相对准确的估计值，成为目前育种中的标杆方法。ML方法中，RFR与特定的ANN方法考虑到标记间的上位效应，ANN在所有的方法中尤其适合大样本数据的处理，未来在基因组预测领域中的应用有待进一步发掘。不同的性状具有不同的遗传结构，在基因组预测领域，目前没有任何一种方法是适合于所有的物种或性状，因此，建议在进行的研究的过程中，将多种方法进行比较，找到最优的解决方案。
| [1] | WALDMANN P. Approximate Bayesian neural networks in genomic prediction[J]. Genet Sel Evol, 2018, 50: 70. DOI: 10.1186/s12711-018-0439-1 |
| [2] | DAVIES R W, FLINT J, MYERS S, et al. Rapid genotype imputation from sequence without reference panels[J]. Nat Genet, 2016, 48(8): 965–969. DOI: 10.1038/ng.3594 |
| [3] | PASANIUC B, ROHLAND N, MCLAREN P J, et al. Extremely low-coverage sequencing and imputation increases power for genome-wide association studies[J]. Nat Genet, 2012, 44(6): 631–635. DOI: 10.1038/ng.2283 |
| [4] | FRAZER K A, BALLINGER D G, COX D R, et al. A second generation human haplotype map of over 3.1 million SNPs[J]. Nature, 2007, 449(7164): 851–861. DOI: 10.1038/nature06258 |
| [5] | STEEMERS F J, GUNDERSON K L. Whole genome genotyping technologies on the BeadArrayTM platform[J]. Biotechnol J, 2007, 2(1): 41–49. DOI: 10.1002/biot.200600213 |
| [6] | QIAO X, SU R, WANG Y, et al. Genome-wide target enrichment-aided chip design:a 66 K SNP chip for cashmere goat[J]. Sci Rep, 2017, 7(1): 8621. DOI: 10.1038/s41598-017-09285-z |
| [7] | GU X R, FENG C G, MA L, et al. Genome-wide association study of body weight in chicken F2 resource population[J]. PLoS One, 2011, 6(7): e21872. DOI: 10.1371/journal.pone.0021872 |
| [8] | ZHU C Y, FAN H Y, YUAN Z H, et al. Detection of selection signatures on the X chromosome in three sheep breeds[J]. Int J Mol Sci, 2015, 16(9): 20360–20374. DOI: 10.3390/ijms160920360 |
| [9] | OSTERSEN T, CHRISTENSEN O F, HENRYON M, et al. Deregressed EBV as the response variable yield more reliable genomic predictions than traditional EBV in pure-bred pigs[J]. Genet Sel Evol, 2011, 43: 38. DOI: 10.1186/1297-9686-43-38 |
| [10] | DIDION J P, YANG H, SHEPPARD K, et al. Discovery of novel variants in genotyping arrays improves genotype retention and reduces ascertainment bias[J]. BMC Genomics, 2012, 13: 34. DOI: 10.1186/1471-2164-13-34 |
| [11] |
王宇哲.利用简化基因组测序精细定位鸡生长性状功能基因[D].北京: 中国农业大学, 2017.
WNG Y Z.Fine mapping the QTL for growth traits in chicken by improves genotyping by sequencing[D].Beijing: China Agricultural University, 2017.(in Chinese) |
| [12] | AVISE J C, LANSMAN R A, SHADE R O. The use of restriction endonucleases to measure mitochondrial dna sequence relatedness in natural populations.I.population structure and evolution in the genus peromyscus[J]. Genomics, 1979, 92(1): 279–295. |
| [13] | BROWN W M. Polymorphism in mitochondrial DNA of humans as revealed by restriction endonuclease analysis[J]. Proc Natl Acad Sci U S A, 1980, 77(6): 3605–3609. DOI: 10.1073/pnas.77.6.3605 |
| [14] | VAN TASSELL C P, SMITH T P L, MATUKUMALLI L K, et al. SNP discovery and allele frequency estimation by deep sequencing of reduced representation libraries[J]. Nat Methods, 2008, 5(3): 247–252. DOI: 10.1038/nmeth.1185 |
| [15] | WIEDMANN R T, SMITH T P, NONNEMAN D J. SNP discovery in swine by reduced representation and high throughput pyrosequencing[J]. BMC Genet, 2008, 9: 81. |
| [16] | MILLER M R, DUNHAM J P, AMORES A, et al. Rapid and cost-effective polymorphism identification and genotyping using restriction site associated DNA (RAD) markers[J]. Genome Res, 2007, 17(2): 240–248. DOI: 10.1101/gr.5681207 |
| [17] | ELSHIRE R J, GLAUBITZ J C, SUN Q, et al. A robust, simple genotyping-by-sequencing (GBS) approach for high diversity species[J]. PLoS One, 2011, 6(5): e19379. DOI: 10.1371/journal.pone.0019379 |
| [18] | ANDREWS K R, GOOD J M, MILLER M R, et al. Harnessing the power of RADseq for ecological and evolutionary genomics[J]. Nat Rev Genet, 2016, 17(2): 81–92. DOI: 10.1038/nrg.2015.28 |
| [19] | HE J F, ZHAO X Q, LAROCHE A, et al. Genotyping-by-sequencing (GBS), an ultimate marker-assisted selection (MAS) tool to accelerate plant breeding[J]. Front Plant Sci, 2014, 5: 484. |
| [20] | WANG Y Z, CAO X M, ZHAO Y Q, et al. Optimized double-digest genotyping by sequencing (ddGBS) method with high-density SNP markers and high genotyping accuracy for chickens[J]. PLoS One, 2017, 12(6): e0179073. DOI: 10.1371/journal.pone.0179073 |
| [21] | TAN C, WU Z F, REN J L, et al. Genome-wide association study and accuracy of genomic prediction for teat number in Duroc pigs using genotyping-by-sequencing[J]. Genet Sel Evol, 2017, 49: 35. DOI: 10.1186/s12711-017-0311-8 |
| [22] | GAUTIER M, GHARBI K, CEZARD T, et al. The effect of RAD allele dropout on the estimation of genetic variation within and between populations[J]. Mol Ecol, 2013, 22(11): 3165–3178. DOI: 10.1111/mec.12089 |
| [23] | SCHWEYEN H, ROZENBERG A, LEESE F. Detection and removal of PCR duplicates in population genomic ddRAD studies by addition of a degenerate base region (DBR) in sequencing adapters[J]. Boil Bull, 2014, 227(2): 146–160. DOI: 10.1086/BBLv227n2p146 |
| [24] | TIN M M Y, RHEINDT F E, CROS E, et al. Degenerate adaptor sequences for detecting PCR duplicates in reduced representation sequencing data improve genotype calling accuracy[J]. Mol Ecol Resour, 2015, 15(2): 329–336. DOI: 10.1111/1755-0998.12314 |
| [25] | TOONEN R J, PURITZ J B, FORSMAN Z H, et al. ezRAD:a simplified method for genomic genotyping in non-model organisms[J]. PeerJ, 2013, 1: e203. DOI: 10.7717/peerj.203 |
| [26] | LITI G, CARTER D M, MOSES A M, et al. Population genomics of domestic and wild yeasts[J]. Nature, 2009, 458(7236): 337–341. DOI: 10.1038/nature07743 |
| [27] | NIELSEN R, PAUL J S, ALBRECHTSEN A, et al. Genotype and SNP calling from next-generation sequencing data[J]. Nat Rev Genet, 2011, 12(6): 443–451. DOI: 10.1038/nrg2986 |
| [28] | HOWIE B, FUCHSBERGER C, STEPHENS M, et al. Fast and accurate genotype imputation in genome-wide association studies through pre-phasing[J]. Nat Genet, 2012, 44(8): 955–959. DOI: 10.1038/ng.2354 |
| [29] | DELANEAU O, ZAGURY J F, MARCHINI J. Improved whole-chromosome phasing for disease and population genetic studies[J]. Nat Methods, 2013, 10(1): 5–6. DOI: 10.1038/nmeth.2307 |
| [30] | BROWNING B L, BROWNING S R. Genotype imputation with millions of reference samples[J]. Am J Hum Genet, 2016, 98(1): 116–126. DOI: 10.1016/j.ajhg.2015.11.020 |
| [31] | HOWIE B N, DONNELLY P, MARCHINI J. A flexible and accurate genotype imputation method for the next generation of genome-wide association studies[J]. PLoS Genet, 2009, 5(6): e1000529. DOI: 10.1371/journal.pgen.1000529 |
| [32] | HUANG X H, WEI X H, SANG T, et al. Genome-wide association studies of 14 agronomic traits in rice landraces[J]. Nat Genet, 2010, 42(11): 961–967. DOI: 10.1038/ng.695 |
| [33] | LIU S Y, HUANG S J, CHEN F, et al. Genomic analyses from non-invasive prenatal testing reveal genetic associations, patterns of viral infections, and chinese population history[J]. Cell, 2018, 175(2): 347–359. DOI: 10.1016/j.cell.2018.08.016 |
| [34] | NICOD J, DAVIES R W, CAI N, et al. Genome-wide association of multiple complex traits in outbred mice by ultra-low-coverage sequencing[J]. Nat Genet, 2016, 48(8): 912–918. DOI: 10.1038/ng.3595 |
| [35] | YANG R F, GUO X L, ZHU D, et al. Genome-wide association analyses of multiple traits in Duroc pigs using low-coverage whole-genome sequencing strategy[J]. Preprint Server Biol, 2019. DOI: 10.1101/754671 |
| [36] | ROS-FREIXEDES R, BATTAGIN M, JOHNSSON M, et al. Impact of index hopping and bias towards the reference allele on accuracy of genotype calls from low-coverage sequencing[J]. Genet Sel Evol, 2018, 50: 64. DOI: 10.1186/s12711-018-0436-4 |
| [37] | TIBSHIRANI R. Regression shrinkage and selection via the lasso[J]. J Roy Stat B, 1996, 58(1): 267–288. |
| [38] | VASQUEZ M M, HU C C, ROE D J, et al. Least absolute shrinkage and selection operator type methods for the identification of serum biomarkers of overweight and obesity:simulation and application[J]. BMC Med Res Methodol, 2016, 16(1): 154. DOI: 10.1186/s12874-016-0254-8 |
| [39] | VASQUEZ M M, HU C C, ROE D J, et al. Measurement error correction in the least absolute shrinkage and selection operator model when validation data are available[J]. Stat Methods Med Res, 2019, 28(3): 670–680. DOI: 10.1177/0962280217734241 |
| [40] | ARBET J, MCGUE M, CHATTERJEE S, et al. Resampling-based tests for Lasso in genome-wide association studies[J]. BMC Genet, 2017, 18: 70. DOI: 10.1186/s12863-017-0533-3 |
| [41] | DE VLAMING R, GROENEN P J F. The current and future use of ridge regression for prediction in quantitative genetics[J]. Biomed Res Int, 2015, 2015: 143712. |
| [42] | VANRADEN P M. Efficient methods to compute genomic predictions[J]. J Dairy Sci, 2008, 91(11): 4414–4423. DOI: 10.3168/jds.2007-0980 |
| [43] | GARCÍA-RUIZ A, COLE J B, VANRADEN P M, et al. Changes in genetic selection differentials and generation intervals in US Holstein dairy cattle as a result of genomic selection[J]. Proc Natl Acad Sci U S A, 2016, 113(28): E3995-–E4004. DOI: 10.1073/pnas.1519061113 |
| [44] | TENG J Y, GAO N, ZHANG H B, et al. Performance of whole genome prediction for growth traits in a crossbred chicken population[J]. Poult Sci, 2019, 98(5): 1968–1975. DOI: 10.3382/ps/pey604 |
| [45] | WANG X, XU Y, HU Z L, et al. Genomic selection methods for crop improvement:current status and prospects[J]. Crop J, 2018, 6(4): 330–340. DOI: 10.1016/j.cj.2018.03.001 |
| [46] | AGUILAR I, MISZTAL I, JOHNSON D L, et al. Hot topic:a unified approach to utilize phenotypic, full pedigree, and genomic information for genetic evaluation of Holstein final score[J]. J Dairy Sci, 2010, 93(2): 743–752. DOI: 10.3168/jds.2009-2730 |
| [47] | FRAGOMENI B O, LOURENCO D A L, LEGARRA A, et al. Alternative SNP weighting for single-step genomic best linear unbiased predictor evaluation of stature in US Holsteins in the presence of selected sequence variants[J]. J Dairy Sci, 2019, 102(11): 10012–10019. DOI: 10.3168/jds.2019-16262 |
| [48] | TEISSIER M, LARROQUE H, ROBERT-GRANIÉ C. Weighted single-step genomic BLUP improves accuracy of genomic breeding values for protein content in French dairy goats:a quantitative trait influenced by a major gene[J]. Genet Sel Evol, 2018, 50: 31. DOI: 10.1186/s12711-018-0400-3 |
| [49] | MEUWISSEN T H, HAYES B J, GODDARD M E. Prediction of total genetic value using genome-wide dense marker maps[J]. Genetics, 2001, 157(4): 1819–1829. |
| [50] | HABIER D, FERNANDO R L, KIZILKAYA K, et al. Extension of the bayesian alphabet for genomic selection[J]. BMC Bioinformatics, 2011, 12: 186. DOI: 10.1186/1471-2105-12-186 |
| [51] | CHEN K, LU Y T, ZHAO H Y, et al. Predicting the change of exon splicing caused by genetic variant using support vector regression[J]. Hum Mutat, 2019, 40(9): 1235–1242. DOI: 10.1002/humu.23785 |
| [52] | HO T K.Random decision forests[C]//Proceedings of 3rd International Conference on Document Analysis and Recognition.Montreal: IEEE, 1995. |
| [53] | BRENNAN R S, HEALY T M, BRYANT H J, et al. Integrative population and physiological genomics reveals mechanisms of adaptation in killifish[J]. Mol Biol Evol, 2018, 35(11): 2639–2653. |
| [54] | MCCULLOCH W S, PITTS W. A logical calculus of the ideas immanent in nervous activity.1943[J]. Bull Math Biol,, 1990, 52(1-2): 99–115. DOI: 10.1016/S0092-8240(05)80006-0 |
| [55] | WU J, ZHAO Y Q. Machine learning technology in the application of genome analysis:a systematic review[J]. Gene, 2019, 705: 149–156. DOI: 10.1016/j.gene.2019.04.062 |
| [56] | MA W L, QIU Z X, SONG J, et al. A deep convolutional neural network approach for predicting phenotypes from genotypes[J]. Plnata, 2018, 248(5): 1307–1318. |
| [57] | ZOU J, HUSS M, ABID A, et al. A primer on deep learning in genomics[J]. Nat Genet, 2019, 51(1): 12–18. DOI: 10.1038/s41588-018-0295-5 |