Julia语言在动物育种中的应用

梅步俊　王志华

摘 要：Julia语言是高性能、动态编译的高级计算机语言，具有极强的灵活性，适合于解决数值和科学计算问题，拥有与传统的静态型语言相媲美的执行速度，Julia语言的开发目的是创建一个功能强大、易用性好和高效的单一语言环境。在动物育种中使用Julia语言可以编写语法简洁，且运行速度接近于Fortran或C++编写的程序。该文提供了7个动物育种中常用程序的Julia代码及示例，包括计算分子血缘相关矩阵（A）、分子血缘相关逆矩阵（A-）、近交系数、设计矩阵、分块矩阵，混合模型方程组（MME）和基因组关系矩阵（G）。这些代码可以为编写动物育种实用程序及相关教学活动奠定基础。

关键词：Julia语言；动物育种；计算生物学

中图分类号 Q819 文献标识码 A 文章编号 1007-7731（2015）18-12-04

Application of Julia Language in Animal Breeding

Mei Bujun1，3 et al.

（1Agriculture department，Hetao College，Bayannur 015000，China；3Department of Animal Science，Iowa State University，Iowa 50010，USA）

Abstract：Julia is a high-performance，high-level dynamic programming language for technical computing.It is a flexible dynamic language，appropriating for numerical and scientific computing，with execution comparable to traditional statically-typed languages.Julia aims to create an quite unusual combination of power，ease-of-use，and efficiency in a single language.The reasons why I have a preference for Julia in animal breeding are：speed and nice syntax.The paper offers 7 source codes in research of animal breeding，including calculating additive relationship matrix or numerator relationship matrix（A），inverse of additive relationship matrix（A-），inbreeding coefficient，design matrix，block or direct matrix，mixed model equation（MME），and genomic relationship matrix coefficients（G）.These programs can be used as a basis for practices of animal breeding and also used for education in animal science.

Key words：Julia language；Animal breeding；Computational biology

随着基因组测序技术的发展，基因组测序成本不断降低，基因组测序数据逐渐在动物育种中广泛使用，这些进展增加了我们对分子水平数量性状的遗传机理的认识，为进一步提高育种效率奠定了基础，特别是对那些使用现行的育种方法效率不高或不能获得理想改良效果的性状。然而，新的理论和方法一般都会涉及大量复杂的运算，这一方面有赖于高性能的计算机硬件设备的发展，另一方面也需要有适应动物育种特点的先进的计算方法。同时，伴随着遗传育种理论和方法的不断发展，新的计算方法也不断出现。一方面，动物育种理论和方法的发展产生了新的计算问题；另一方面，不断涌现的新的计算方法又催生了动物育种理论和方法的新发展。因此，计算技术、方法的研究一直是动物育种理论研究和应用研究中不可或缺的关键技术领域。不掌握这些技术方法，就不具备真正理解现代动物育种理论和方法的基础，也就难以开展较为深入的研究。

虽然有许多现成的软件或程序可以解决动物育种中的诸多问题，但是由于实践中会出现林林总总的计算问题，编写程序仍然是育种工作者或育种理论研究者的必备技能。同时，由于新的计算理论、技术和算法层出不穷，所以在很多情况下，没有现成的软件或程序可以实现动物育种学研究者创造新的方法或改善现有计算效果的意图。因此，掌握若干计算机语言，并能用其解决育种学问题，往往是研究动物育种学前沿问题的基础。据统计，目前广泛使用的计算机语言约有91种，依据这些语言的不同特点及不同研究领域的传统，特定领域会使用不同的计算机语言。美国农业部资助的“Animal Genome”数据库项目（http：//www.animalgenome.org/）收集了329种遗传学分析软件。在动物育种中，广泛使用的计算机语言有C++（包括C）、Fortran、Java、MATLAB、AWK、Python、Visual Basic、R、Perl等，这些语言可初略的分为编译型语言和解释性语言。前者程序执行速度快，但对于一般的动物育种学研究者而言学习及编写程序的难度较大，开发周期也相对较长。后者对不同系统平台的兼容性较好，借助特定的函数库，开发特定程序的周期较短，但此类语言在运行程序时需要专门有一个解释器，每个语句都是在执行的时候才翻译，执行一次就要翻译一次，因此效率低。但这些区别也不能一概而论，部分解释型语言的解释器，通过在运行时动态优化代码，甚至能够获的超过编译型语言的性能。

1 Julia语言的特点

Julia语言受NumFOCUS资助，其创始人为若干精通Matlab科学计算的编程人员，创立此项目的初衷据称是由于不满意现有的编程工具。该项目大约于2009年开始，目前的版本为v0.3.11，其源代码，及各种平台的可执行文件及专业编译器Juno可在http：//julialang.org网站下载[1]。Julia语言可以通过基于网页的Jupyter（IJulia）交互环境执行，方便在教学等情景下展示执行结果。Julia是新的高性能、编译型、动态交互式的高级编程语言，集中了许多计算机语言的优点，“它拥有类似于C语言一样的执行速度，拥有如同Ruby语言的动态性，又有Matlab般熟悉的数学记号和线性代数运算能力，兼具像Python般的通用性，又像R语言一样擅长于统计分析，并有Perl般处理字符串的能力和shell等胶水语言的特点，并易于学习”。目前已有多所国际知名高校的数值计算或统计学课程结合Julia语言进行讲解，如斯坦福大学的“应用矩阵方法（Introduction to Matrix Methods；课程代码：EE103）”和麻省理工学院的“线性代数（Linear Algebra；课程代码：18.06）”。爱荷华州立大学动物科学系2015年5月在其开设的“家畜基因组预测（Genomic Prediction in Livestock）”短期课程中结合Julia语言进行了讲解。使用7种标准检查程序，Julia语言的运行速度接近于C及Fortran语言（见图1），但其编写数值计算程序的速度却快得多。一般情况下，Julia语言运行数值计算程序时的速度也接近于C++，是R语言速度的100倍，MATLAB语言的1 000倍。

注：设C语言的运行时间为1.0；C和Fortran语言使用gcc 4.8.2进行编译；C、Fortran和Julia使用OpenBLAS v0.2.13；Python运行rand_mat_stat和rand_mat_mul使用NumPy v1.9.2库函数。

2 Julia语言在动物育种中的应用

2.1 分子血缘相关矩阵（A）及其逆矩阵计算 分子血缘相关矩阵或加性遗传相关矩阵及其逆矩阵计算在动物育种学教学及种畜的遗传评估中有重要意义。使用Julia语言计算分子血缘相关矩阵需构建“Pedigree”和“PedNode”类型，使用“mkPedigree”函数对系谱中的个体进行排序，使子代位于亲代之后。使用“Amatrix”函数计算分子血缘相关矩阵[2]，其输出结果使用稀疏矩阵存储，可使用“round（full（A），2）”转化为保留2位小数的满矩阵。以图2所示系谱为例，其重新编码的系谱如下：

2.2 近交系数计算 近交系数为形成合子的2个配子源于同一共同祖先的概率。由通径系数原理可知个体X的近交系数即为形成X个体的两个配子间的相关系数[4]，用一般用Fx表示。使用“Inbreeding”函数计算10个个体的近交系数为：

2.3 设计矩阵 统计学中，设计矩阵（Design matrix）的元素为解释变量（Explanatory variable），如在方差分析中用指示变量（1和0）表示连续变量在模型中的位置[5]。用“design”函数可以构建设计矩阵，如5头奶牛分别养殖在3个牧场（1，1，2，3，2），则由“design”函数构建的设计矩阵以稀疏矩阵形式存储，转换满矩阵为：

2.5 混合模型方程组（MME）建立 1953年，C.R.Henderson以混合模型为基础建立线性方程组。由此方程组求解，可得到固定效应的最佳线性无偏估计和随机效应的最佳线性无偏预测[7]。混合模型方程组是解决许多动物育种学问题的基础。由“MME”函数可以建立混合模型方程组的“左手项”和“右手项”，并得出相应参数的估计量。

[X'R-XX'R-1ZZ'R-1XZ'R-1Z+G-1bu=X'R-1yZ'R-1y]

如：X=[1.0，1.0，1.0，1.0，1.0，1.0]，；Z=[1，2，3，1，2，1]，可由上述“design”函数转化为设计矩阵；GI为3×3单位矩阵；RI为6×6单位矩阵；y=[2.0，1.5，2.0，1.2，0.89，1.2]，由“MME”函数计算得到，剩余方差（SSR）为13.121 8；参数[b]（[u]），“右手项（RHS）”，“左手项（LHS）”分别为：

2.6 基因组关系矩阵（G）的建立 “computeG”函数可以计算VanRaden（2008）或Yang等（2010）定义的G矩阵，并可以使用等位基因频率平均值或2倍基因频率（需满足Hardy-Weinberg平衡）进行中心化[8]。

3 结论

实践中，由于R语言开发育种程序的速度较快，但程序运行速度较慢，其编写的程序很难用于实际育种工作。笔者只是用R语言测试统计方法的可行性，在此基础上再将R程序转变为Fortran或C++等程序，整个开发过程实际经历了2个内容基本相同的阶段。Julia语言兼具有R语言和Fortran（或C++）的优点，有效地缩短了程序开发时间。由Julia语言编写的7个动物育种中常用的程序（A矩阵计算、A逆矩阵计算、近交系数计算、设计矩阵、块矩阵和、MME矩阵、G矩阵计算），可用于动物育种学应用程序的编写和教学，源代码可以通过邮件（meibujun@163.com或meibujun@iastate.edu）向作者索取。经过测试，7个程序的计算结果可靠。这些程序可以作为广大动物育种工作者掌握Julia语言应用于动物育种学的基础。

致谢：本研究部分灵感及部分计算设备由中国农业大学动物科技学院张勤教授课题组提供；感谢美国爱荷华州立大学动物科学系Rohan L.Fernando教授、Hao Chen博士和Jian Zeng博士在编写程序过程中的帮助。

参考文献

[1]Balbaert，I.，Getting Started with Julia Programming[J].2015：PacktLib.214.

[2]Ter Heijden，E.，J.P.Chesnais，C.G.Hickman，An efficient method of computing the numerator relationship matrix and its inverse matrix with inbreeding for large sets of animals[J].Theor Appl Genet.，1977，49（5）：237-241.

[3]Oikawa，T.and K.Yasuda，Inclusion of genetically identical animals to a numerator relationship matrix and modification of its inverse[J].Genet Sel Evol.，2009，41：25.

[4]Boucher，W.，Calculation of the inbreeding coefficient[J].J Math Biol.，1988，26（1）：p.57-64.

[5]Matsumoto，Y.and Y.Kuroyanagi，Design of a matrix for cultured dermal substitute suitable for simultaneous transplantation with auto-skin graft：evaluation in animal test[J].J Biomater Sci Polym Ed.，2010，21（1）：83-94.

[6]Drake，M.D.，PLZT Matrix-Type Block Data Composerst[J].Appl Opt.1974，13（2）：347-352.

[7]Cheung，M.W.，A model for integrating fixed-，random-，and mixed-effects meta-analyses into structural equation modeling[J].Psychol Methods，2008，13（3）：182-202.

[8]Meyer，K.，B.Tier，and H.U.Graser，Technical note：updating the inverse of the genomic relationship matrix[J].J Anim Sci.，2013，91（6）：2583-2586. （责编：张宏民）