2023年1月30日 收稿; 2023年5月5日 收修改稿 摘要 : 在育种中，常常通过利用单核苷酸多态性(SNPs)来预测表型以辅助育种，提高育种效率。传统的统计分析方法受到数据缺失等诸多因素的限制，在一些情况下效果不佳。针对此问题，提出一种利用多尺度特征进行植物性状预测的卷积神经网络模型(MSF-CNN)，该模型通过卷积提取3个不同尺度的SNPs特征，对植物性状数值进行回归预测，并通过对模型中SNPs的权重分析SNP位点的显著性。测试结果表明，与目前已知的其他方法相比，MSF-CNN模型在有基因型数据缺失值的数据集上表型预测的准确性更高。此外，通过显著性图研究基因型对性状的贡献，发现数个较显著的SNP位点。说明该深度学习模型可以更准确地预测定量表型，并能够高效识别与全基因组关联研究相关的SNP位点。 Multi-scale featured convolution neural network-based soybean phenotypic prediction Abstract : In breeding, single nucleotide polymorphisms (SNPs) in the genome are often used to predict quantitative phenotypes to assist breeding, thereby improving breeding efficiency. The traditional statistical analysis method is limited by many factors including missing data, and its performance sometimes can not meet the requirements. In this paper, we proposed a multi-scale feature convolutional neural network model (MSF-CNN) to predict plant traits. The model extracted SNP features at three different scales through convolution and analyzed the significance of SNP sites through the weight of the SNPs input into the model. The test results showed that MSF-CNN model performed with higher accuracy than the known methods and other deep learning models in phenotype prediction on the datasets with missing genotypic data. This paper also studied the contribution of genotype to traits through saliency map, and discovered several significant SNP loci. These results showed that, compared with other known methods available at present, the deep learning model proposed in this paper can obtain more accurate prediction results of quantitative phenotypes, and can also effectively and efficiently identify SNPs associated with genome-wide association research.

在研究基因型，特别是单核苷酸多态性(single nucleotide polymorphism, SNP)与表型之间的关联以及指导育种时，全基因组关联分析(genome-wide association study, GWAS)对数量性状的表型预测有着明显的作用。GWAS目前已广泛应用于许多主要农作物的育种工作中 ^{[ 1 - 3 ]} ，如大豆( Glycine max )、水稻( Oryza sativa )和玉米( Zea mays )等。传统的统计方法，如最佳线性无偏预测、贝叶斯模型 ^{[ 4 ]} 等，目前被广泛用于预测基因型效应和预测表型。传统统计方法通常假设基因型随机效应遵循高斯分布等先验分布，同时将每个基因型对相关表型的贡献视为相互独立的因素。这些先验假设需要足够大的训练样本来修正和调整。然而，在实践中，单个基因的效应是未知的，可能不严格遵循先验分布。此外，SNPs还可与其他SNPs发生相互作用，通过上位性效应，导致复杂的疾病或性状 ^{[ 5 ]} 。因此，传统统计方法的准确性往往受到数据量和各种复杂遗传因素的限制。

此外，对于统计方法而言，基因型数据中的缺失值会带来巨大制约。通常，这些缺失值会在预处理过程中筛选出来，并通过插补数据进行填充 ^{[ 6 ]} 。插补是一种用于估计模板群体中基因型的缺失值的算法，目前已有几种具有或不具有参考基因组信息的基因组插补方法。插补精度高度依赖于观察到的非缺失基因型和整个群体的缺失率，这直接影响表型预测模型的性能 ^{[ 6 ]} 。通过统计方法进行表型预测的模型，需要将基因型矩阵一起估算，然后将其划分为用于模型训练和测试的训练和测试数据集。在某种程度上，测试集并非完全独立于训练集，因为在这种情况下，训练集可能包含从测试集数据中进行估计和插补产生的基因型。不准确的插补方法也可能带来误差和不确定性。因此，这些插补方法可能无法有效推断隐藏在整个基因组中的信息。

近年来，深度学习在生物学领域逐渐得到广泛的应用，在计算生物学中的一些重大挑战问题上取得了前所未有的进展。目前，深度学习已经在蛋白质结构预测、蛋白质功能预测、基因组工程、系统生物学和数据集成以及系统发育预测等5个领域得到广泛应用 ^{[ 7 ]} 。在部分领域，如在蛋白质功能预测方面，深度学习的性能明显超过其他机器学习模型和经典方法，并产生了广泛影响。在系统生物学和基因型关联研究中，与现有方法相比，深度学习可以整合和集成多维度的数据信息，发现具有不同数据模式特征的有意义的亚组 ^{[ 8 ]} ，识别SNP相互作用 ^{[ 9 ]} 以及对基因组变异进行分类 ^{[ 10 ]} 。深度学习可以利用各种不同来源数据，对其进行联合建模。对于表型预测问题，神经网络可以从原始测序数据或基因型数据中直接获取信息。目前，已有多篇文献在表型预测领域中应用了深度学习，如通过癌细胞的多组学信息预测乳腺癌患者的5年生存期 ^{[ 8 ]} 、通过生物标志物信息预测患者的阿尔兹海默症进展 ^{[ 11 ]} 等，但多见于医学领域。在植物遗传学和育种方面尚未得到有效的应用，仅有少量使用神经网络预测进行辅助育种的例子 ^{[ 12 ]} ，但很少有使用大量SNP数据的研究，且并未采用较为先进的深度学习网络。然而，上述深度学习方法同样有局限性，它们并没有有效地解决数据缺失和上位效应等问题。具体而言，深度学习可能会过度拟合数据，并且容易在训练集和测试集之间产生“污染”，这些都可能导致对性能的高估 ^{[ 13 ]} 。

针对传统统计分析方法和已知深度学习方法存在的问题，本文提出一个多尺度特征卷积神经网络(muti-scale feature convolutional neural network, MSF-CNN)模型，通过提取SNPs位点分布的多尺度特征预测对应植物的表型，并根据训练后模型的权重评估SNPs位点的显著性。本文利用大豆数据集分别训练株高、含水量、含油量、蛋白质含量、产量等5个表型所对应的模型，并进行模型评估。同时，分析对含油量和蛋白质含量等表型有贡献的SNPs的权重，并与Wald检验 ^{[ 14 ]} 的结果进行对照。结果表明，该模型与传统统计方法以及现有深度学习相比，在未插补数据集上的性状预测性能更优，在SNPs显著性分析方面也有较好的效果。

1 实验方法和材料 1.1 数据集

本文以大豆作为研究对象，使用大豆的实验数据集作为基准来评估所提出模型的性能。大豆数据集来自SoyBase数据库中的SoyNAM项目，该项目是目前最完善的大豆SNP数据集 ^{[ 14 - 15 ]} 。通过对其中数据的质量控制和筛选，最终得到一个在包含5 000多个重组自交系中发现的4 236个常见SNP，以及相应的大豆株系的株高、含水量、含油量、蛋白质含量和产量等性状表型的数据集。数据集中缺失的基因型数据通过MaCH软件进行估计和补齐，与随机森林回归法和期望最大化插补法等其他插补方法相比，该插补方法误差更小、预测的准确性更高 ^{[ 16 ]} 。基因组选择的效果依赖于基因标记的覆盖率，因此为了以低成本获得大量基因标记，通常会使用较低的测序深度，难免产生大量位点缺失的数据，而大多数分析需要完整的数据。因此，在将测序基因分型的数据用于分析之前，插补往往是必要的步骤。研究表明，对缺失数据的插补可以提高遗传关联研究的能力 ^{[ 6 ]} 。而在深度学习网络的训练中，使用一个独立的编码代表数据缺失，因此可以使用未经插补的原始数据直接进行训练，而不需经过插补。我们在训练模型时，使用经插补的数据和未经插补的数据分别进行训练，并比较了两者的训练效果。

1.1.1 SNP表示

本文对性状数据进行了规整化以便训练，将其均转换为(-1，1)区间内的数值。模型中，输入的性状数据为一个列表，其内容包含植株个体的编号、每一个个体在4 236个常见SNP位点的基因型，及其性状数据。SNP按其所在染色体及染色体上的位置按顺序进行排序，依次由1号染色体列至20号染色体。每一个SNP可被视为输入矩阵中独立的一行，因此SNP的顺序不影响训练的结果。对于每一个性状，输入的SNP数据均相同。使用多个二进制编码构成的向量代表一个SNP位点的3种基因型(0、1、2)以及数据缺失(-1)，并以此作为输入向量的一部分。每一种基因型由一个对应位置为1的4维向量表示，其余的设置为0。例如，3种基因型[AA，Aa，aa]分别表示为[0, 1, 0, 0]，[0, 0, 0, 1]和[0, 0, 1, 0]。缺失的基因型表示为[1, 0, 0, 0]。由于SNP数据已通过对基因组序列的预处理进行了提炼，SNP数据本身并不包含基因序列信息，仅为区分3种不同基因型[AA，Aa，aa]的标记，因此不会受基因序列的同源性影响，避免了训练集、验证集、测试集之间的样本交叉。

1.1.2 性状数据标准化

对表型数据进行标准化，将其转化为 z 分数。其公式如下

1.2 MSF-CNN 1.2.1 模型结构

密集卷积网络(dense convolutional network, DenseNet)是一种较新的网络结构 ^{[ 17 ]} 。其特征为其每一层都与其他层相连。每一层使用所有先前层的特征图作为输入，其自身的特征图作为所有后续层的输入。因为卷积网络在靠近输入的层和靠近输出的层之间包含较短的连接，因此卷积网络可以更深入、更准确、更有效地训练。

本文借鉴DenseNet结构的特点，构建MSF-CNN，其结构如 图 1 所示。模型包括1个输入层、4个卷积层、2个拼接层、1个平坦层、1个全连接层和1个输出。输入数据向量包含4 236个SNP；每个SNP位点对应1个如2.1.1节所述的4维向量，因此输入层所输入的数据为1个4 236×4的矩阵。在经过第1个卷积层的处理后，转变为4 236×12的矩阵。随后，在卷积层后又加入池化层，得到大小为2 118×12的矩阵。以上2个矩阵通过接合层进行接合，最终在第1个拼接层产生1个6 354×12的矩阵。同理，当这一矩阵再次进行后续的卷积、池化、拼接过程后，在第2个拼接层生成1个9 531×12的矩阵。该矩阵由多层次之间的数据接合形成，因此具备多尺度下的数据特征，可以比过往模型更有效地挖掘和预测数据之间存在的多尺度关联与相互作用。在2次拼接之后，数据再次通过一次卷积进行后处理，并输出到平坦层。平坦层将特征矩阵展开为1维向量并传递到最后一个全连接层，输出最终预测的表型。

\begin{equation*} \text { MSE }=\frac{1}{n} \sum _{i=1}^{n}\left(\widetilde{Y_{i}}-Y_{i}\right)^{2} . \end{equation*} \begin{equation*} r=\frac{\sum\limits_{i=1}^{n}\left(X_{i}-\bar{X}\right)\left(Y_{i}-\bar{Y}\right)}{\sqrt{\sum\limits_{i=1}^{n}\left(X_{i}-\bar{X}\right)^{2}} \sqrt{\sum\limits_{i=1}^{n}\left(Y_{i}-\bar{Y}\right)^{2}}} . \end{equation*} $\begin{equation*} S_{c}(I)=\boldsymbol{\omega}_{c}^{\mathrm{T}} \boldsymbol{I}+b_{c} . \end{equation*}

在以已识别SNP为中心的20 kbp区域内找到了最近的基因，并标注了其名称和功能。根据基因模型“Glyma.Wm82.a1.v.1” ^{[ 22 ]} ，使用Soybase数据库 ^{[ 23 ]} 和SoyKB ^{[ 24 ]} 的蛋白质家族(PFAM)数据库 ^{[ 25 ]} 进行注释。基因注释表明，图中标注的SNP位点及其附近区域与它们的性状高度相关。另外，也检测到几个新的标记和区域。

通过对大豆产量预测模型的显著性分析，发现Gm02_6396340、Gm03_5730188、Gm03_5816390、Gm04_2571383等几个SNP在模型中对性状的影响最为显著。考察与这些SNP位点距离最接近的基因，发现Gm02_6396340的位点与大豆基因 Glyma.02g073100 重叠，该基因编码一个包含成簇蛋白结合域的高GC区DNA结合蛋白(PF07842)；Gm03_5730188的位点接近基因 Glyma.03g045300 ，其编码一个包含NB-ARC区域的蛋白(PF00931)；Gm03_5816390的位点与基因 Glyma.03g045700 重叠，其编码一个包含NB-ARC区域以及LRR的蛋白(PF00560)；Gm04_2571383的位点接近基因 Glyma.04g032000 ，其编码一个亮氨酸羧基甲基转移酶家族的蛋白(PF04072)。

通过对大豆株高预测模型的显著性分析，发现Gm10_44287415、Gm10_44500915、Gm10_ 44669893、Gm11_1266080等几个SNP在模型中对性状的影响最为显著。位于大豆10号染色体上的一个区域含有多个显著位点。考察与该SNP位点距离最接近的基因，发现Gm10_44287415的位点与基因 Glyma.10g210600 重叠，其编码一个生长素响应因子(PF06507)；Gm10_44500915的位点与基因 Glyma.10g212500 重叠，其编码一个丝氨酸羧肽酶家族蛋白(PF00450)；Gm10_ 44669893接近基因 S58482.1 ，其表达被生长素下调；Gm11_1266080的位点接近基因 AP2-9 ，其编码一个AP2-EREBP转录因子，AP2家族基因据报道可调节大豆的生长素合成，当其过表达时会提高大豆株高 ^{[ 26 ]} 。

通过对大豆含水量预测模型的显著性分析，发现Gm10_44124696、Gm10_44287415、Gm15_9530676、Gm17_6781998等几个SNP在模型中对性状的影响最为显著。考察与该SNP位点距离最接近的基因，其中Gm10_44124696接近基因 Glyma.10g209200 ，其编码一个ATP结合盒转运蛋白(PF00005)；Gm10_44287415已在株高一节中提及；Gm15_9530676接近基因 Glyma.15g120300 ，其编码一个姐妹染色单体凝聚蛋白；最接近Gm17_6781998位点的基因 Glyma17g09165 属于蛋白激酶结构域(PF00069)，并参与对寒冷、创伤、盐胁迫和甘露醇刺激等反应的生物学过程。

通过对大豆蛋白质含量预测模型的显著性分析，发现Gm07_7832406、Gm18_2444660、Gm20_29976653等几个SNP在模型中对性状的影响最为显著。同时，它们在Wald检验中同样显示有较强的显著性。考察与上述SNP位点距离最接近的基因，其中Gm07_7832406位点与基因 Glyma.07g085000 重叠，其编码一个植物特异性的Rop蛋白(PF03759)；Gm18_2444660位点接近基因 Glyma.18g031700 ；Gm20_29976653位点接近基因 Glyma.20g080300 ，其编码一个磷酸酶家族蛋白(PF03372)。

通过对大豆含油量预测模型的显著性分析，发现Gm02_7439248、Gm04_8184443、Gm18_1685024等几个SNP在模型中对性状的影响最为显著。同时，它们在Wald检验中也显示有较强的显著性。Gm02_7439248位点与基因 Glyma.02g085400 重叠，其编码一个起RNA代谢功能的金属β-内酰胺酶(PF07521)；Gm04_8184443与基因 Glyma04g09900 接近，该基因属于蛋白酪氨酸激酶家族(PF07714)，涉及蛋白磷酸化过程和寡肽转运过程；Gm18_1685024位点接近基因 Glyma.18g023100 ，其编码一个预苯酸脱氢酶(PF02153)。

3 结论

本文提出一种基于密集卷积网络的深度学习模型MSF-CNN，可以通过SNP标记准确预测表型，而不需要对原始数据中缺失的基因型进行插补。探索几种不同的深度学习架构，并对原本的密集卷积网络模型进行简化和优化，最终得到优化的模型结构。与目前常用的神经网络模型相比，该模型在真实实验数据集上具有最佳的预测性能。通过显著性图研究不同位点和基因对性状的贡献，发现了数个较显著的SNP位点。该模型还需要继续改进，使其具备明确所研究基因型间的相互作用对表型的影响、解释预测结果的潜在生物学意义的能力，为表型预测和SNP选择发挥作用。

Zhao Y S, Gowda M, Liu W X, et al. Accuracy of genomic selection in European maize elite breeding populations[J]. Theoretical and Applied Genetics, 2012, 124(4): 769-776. Doi:10.1007/s00122-011-1745-y Spindel J, Begum H, Akdemir D, et al. Genomic selection and association mapping in rice ( Oryza sativa ): effect of trait genetic architecture, training population composition, marker number and statistical model on accuracy of rice genomic selection in elite, tropical rice breeding lines[J]. PLoS Genetics, 2015, 11(2): e1004982. Doi:10.1371/journal.pgen.1004982 Xavier A, Jarquin D, Howard R, et al. Genome-wide analysis of grain yield stability and environmental interactions in a multiparental soybean population[J]. G3-Genes Genomes Genetics, 2018, 8(2): 519-529. Doi:10.1534/g3.117.300300 Endelman J B. Ridge regression and other kernels for genomic selection with R package rrBLUP[J]. The Plant Genome, 2011, 4(3): 250-255. Doi:10.3835/plantgenome2011.08.0024 Wang J X, Joshi T, Valliyodan B, et al. A Bayesian model for detection of high-order interactions among genetic variants in genome-wide association studies[J]. BMC Genomics, 2015, 16: 1011. Doi:10.1186/s12864-015-2217-6 Rutkoski J E, Poland J, Jannink J L, et al. Imputation of unordered markers and the impact on genomic selection accuracy[J]. G3 Genes|Genomes|Genetics, 2013, 3(3): 427-439. Doi:10.1534/g3.112.005363 Sapoval N, Aghazadeh A, Nute M G, et al. Current progress and open challenges for applying deep learning across the biosciences[J]. Nature Communications, 2022, 13(1): 1-12. Doi:10.1038/s41467-022-29268-7 Tong L, Mitchel J, Chatlin K, et al. Deep learning based feature-level integration of multi-omics data for breast cancer patients survival analysis[J]. BMC Medical Informatics and Decision Making, 2020, 20(1): 225. Doi:10.1186/s12911-020-01225-8 Uppu S, Krishna A, Gopalan R P. A deep learning approach to detect SNP interactions[J]. Journal of Software, 2016, 11(10): 965-975. Doi:10.17706/jsw.11.10.965-975 Liang Z H, Huang J X, Zeng X, et al. DL-ADR: a novel deep learning model for classifying genomic variants into adverse drug reactions[J]. BMC Medical Genomics, 2016, 9(S2): 48. Doi:10.1186/s12920-016-0207-4 Lee G, Nho K, Kang B, et al. Predicting Alzheimer's disease progression using multi-modal deep learning approach[J]. Scientific Reports, 2019, 9(1): 1-12. Doi:10.1038/s41598-018-37769-z Zingaretti L M, Gezan S A, Ferrão L F V, et al. Exploring deep learning for complex trait genomic prediction in polyploid outcrossing species[J]. Frontiers in Plant Science, 2020, 11: 25. Doi:10.3389/fpls.2020.00025 Whalen S, Schreiber J, Noble W S, et al. Navigating the pitfalls of applying machine learning in genomics[J]. Nature Reviews Genetics, 2022, 23(3): 169-181. Doi:10.1038/s41576-021-00434-9 Xavier A, Beavis W D, Specht J E, et al. SoyNAM: soybean nested association mapping dataset[DB]. R package version, 2015, 1. Song Q J, Yan L, Quigley C, et al. Genetic characterization of the soybean nested association mapping population[J]. The Plant Genome, 2017, 10(2): 10. Doi:10.3835/plantgenome2016.10.0109 Li Y, Willer C J, Ding J, et al. MaCH: using sequence and genotype data to estimate haplotypes and unobserved genotypes[J]. Genetic Epidemiology, 2010, 34(8): 816-834. Doi:10.1002/gepi.20533 Huang G, Liu Z, Van Der Maaten L, et al. Densely connected convolutional networks[C]. 2017 IEEE Conference on Computer Vision and Pattern Recognition (CVPR). July 21-26, 2017, Honolulu, HI, USA. IEEE, 2017: 2261-2269. DOI: 10.1109/CVPR.2017.243 . Zhang X Z, Wang Y F, Shi W S. pCAMP: performance comparison of machine learning packages on the edges[EB/OL]. (2019-06-05)[2023-03-24]. https://arxiv.org/abs/1906.01878 . Simonyan K, Zisserman A. Very deep convolutional networks for large-scale image recognition[EB/OL]. (2014-09-04)[2023-03-24]. https://arxiv.org/abs/1409.1556 . Srivastava N, Hinton G, Krizhevsky A, et al. Dropout: a simple way to prevent neural networks from overfitting[J]. The Journal of Machine Learning Research, 2014, 15(1): 1929-1958. Liu Y, Wang D L, He F, et al. Phenotype prediction and genome-wide association study using deep convolutional neural network of soybean[J]. Frontiers in Genetics, 2019, 10: 1091. Doi:10.3389/fgene.2019.01091 Schmutz J, Cannon S B, Schlueter J, et al. Erratum: genome sequence of the palaeopolyploid soybean[J]. Nature, 2010, 465(7294): 120. Doi:10.1038/nature08957 Grant D, Nelson R T, Cannon S B, et al. SoyBase, the USDA-ARS soybean genetics and genomics database[J]. Nucleic Acids Research, 2010, 38(Suppl 1): D843-D846. Doi:10.1093/nar/gkp798 Joshi T, Patil K, Fitzpatrick M R, et al. Soybean knowledge base (SoyKB): a web resource for soybean translational genomics[J]. BMC Genomics, 2012, 13(Suppl 1): S15. Doi:10.1186/1471-2164-13-S1-S15 Bateman A, Coin L, Durbin R, et al. The Pfam protein families database[J]. Nucleic Acids Research, 2004, 32(Suppl 1): D138-D141. Doi:10.1093/nar/gkh121 Xu Z Y, Wang R K, Kong K K, et al. An APETALA2/ethylene responsive factor transcription factor GmCRF4a regulates plant height and auxin biosynthesis in soybean[J]. Frontiers in Plant Science, 2022, 13: 983650. Doi:10.3389/fpls.2022.983650