Julia語言在動(dòng)物育種中的應(yīng)用

梅步俊　王志華

摘 要：Julia語言是高性能、動(dòng)態(tài)編譯的高級(jí)計(jì)算機(jī)語言，具有極強(qiáng)的靈活性，適合于解決數(shù)值和科學(xué)計(jì)算問題，擁有與傳統(tǒng)的靜態(tài)型語言相媲美的執(zhí)行速度，Julia語言的開發(fā)目的是創(chuàng)建一個(gè)功能強(qiáng)大、易用性好和高效的單一語言環(huán)境。在動(dòng)物育種中使用Julia語言可以編寫語法簡潔，且運(yùn)行速度接近于Fortran或C++編寫的程序。該文提供了7個(gè)動(dòng)物育種中常用程序的Julia代碼及示例，包括計(jì)算分子血緣相關(guān)矩陣（A）、分子血緣相關(guān)逆矩陣（A-）、近交系數(shù)、設(shè)計(jì)矩陣、分塊矩陣，混合模型方程組（MME）和基因組關(guān)系矩陣（G）。這些代碼可以為編寫動(dòng)物育種實(shí)用程序及相關(guān)教學(xué)活動(dòng)奠定基礎(chǔ)。

關(guān)鍵詞：Julia語言；動(dòng)物育種；計(jì)算生物學(xué)

中圖分類號(hào) Q819 文獻(xiàn)標(biāo)識(shí)碼 A 文章編號(hào) 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

隨著基因組測序技術(shù)的發(fā)展，基因組測序成本不斷降低，基因組測序數(shù)據(jù)逐漸在動(dòng)物育種中廣泛使用，這些進(jìn)展增加了我們對(duì)分子水平數(shù)量性狀的遺傳機(jī)理的認(rèn)識(shí)，為進(jìn)一步提高育種效率奠定了基礎(chǔ)，特別是對(duì)那些使用現(xiàn)行的育種方法效率不高或不能獲得理想改良效果的性狀。然而，新的理論和方法一般都會(huì)涉及大量復(fù)雜的運(yùn)算，這一方面有賴于高性能的計(jì)算機(jī)硬件設(shè)備的發(fā)展，另一方面也需要有適應(yīng)動(dòng)物育種特點(diǎn)的先進(jìn)的計(jì)算方法。同時(shí)，伴隨著遺傳育種理論和方法的不斷發(fā)展，新的計(jì)算方法也不斷出現(xiàn)。一方面，動(dòng)物育種理論和方法的發(fā)展產(chǎn)生了新的計(jì)算問題；另一方面，不斷涌現(xiàn)的新的計(jì)算方法又催生了動(dòng)物育種理論和方法的新發(fā)展。因此，計(jì)算技術(shù)、方法的研究一直是動(dòng)物育種理論研究和應(yīng)用研究中不可或缺的關(guān)鍵技術(shù)領(lǐng)域。不掌握這些技術(shù)方法，就不具備真正理解現(xiàn)代動(dòng)物育種理論和方法的基礎(chǔ)，也就難以開展較為深入的研究。

雖然有許多現(xiàn)成的軟件或程序可以解決動(dòng)物育種中的諸多問題，但是由于實(shí)踐中會(huì)出現(xiàn)林林總總的計(jì)算問題，編寫程序仍然是育種工作者或育種理論研究者的必備技能。同時(shí)，由于新的計(jì)算理論、技術(shù)和算法層出不窮，所以在很多情況下，沒有現(xiàn)成的軟件或程序可以實(shí)現(xiàn)動(dòng)物育種學(xué)研究者創(chuàng)造新的方法或改善現(xiàn)有計(jì)算效果的意圖。因此，掌握若干計(jì)算機(jī)語言，并能用其解決育種學(xué)問題，往往是研究動(dòng)物育種學(xué)前沿問題的基礎(chǔ)。據(jù)統(tǒng)計(jì)，目前廣泛使用的計(jì)算機(jī)語言約有91種，依據(jù)這些語言的不同特點(diǎn)及不同研究領(lǐng)域的傳統(tǒng)，特定領(lǐng)域會(huì)使用不同的計(jì)算機(jī)語言。美國農(nóng)業(yè)部資助的“Animal Genome”數(shù)據(jù)庫項(xiàng)目（http：//www.animalgenome.org/）收集了329種遺傳學(xué)分析軟件。在動(dòng)物育種中，廣泛使用的計(jì)算機(jī)語言有C++（包括C）、Fortran、Java、MATLAB、AWK、Python、Visual Basic、R、Perl等，這些語言可初略的分為編譯型語言和解釋性語言。前者程序執(zhí)行速度快，但對(duì)于一般的動(dòng)物育種學(xué)研究者而言學(xué)習(xí)及編寫程序的難度較大，開發(fā)周期也相對(duì)較長。后者對(duì)不同系統(tǒng)平臺(tái)的兼容性較好，借助特定的函數(shù)庫，開發(fā)特定程序的周期較短，但此類語言在運(yùn)行程序時(shí)需要專門有一個(gè)解釋器，每個(gè)語句都是在執(zhí)行的時(shí)候才翻譯，執(zhí)行一次就要翻譯一次，因此效率低。但這些區(qū)別也不能一概而論，部分解釋型語言的解釋器，通過在運(yùn)行時(shí)動(dòng)態(tài)優(yōu)化代碼，甚至能夠獲的超過編譯型語言的性能。

1 Julia語言的特點(diǎn)

Julia語言受NumFOCUS資助，其創(chuàng)始人為若干精通Matlab科學(xué)計(jì)算的編程人員，創(chuàng)立此項(xiàng)目的初衷據(jù)稱是由于不滿意現(xiàn)有的編程工具。該項(xiàng)目大約于2009年開始，目前的版本為v0.3.11，其源代碼，及各種平臺(tái)的可執(zhí)行文件及專業(yè)編譯器Juno可在http：//julialang.org網(wǎng)站下載[1]。Julia語言可以通過基于網(wǎng)頁的Jupyter（IJulia）交互環(huán)境執(zhí)行，方便在教學(xué)等情景下展示執(zhí)行結(jié)果。Julia是新的高性能、編譯型、動(dòng)態(tài)交互式的高級(jí)編程語言，集中了許多計(jì)算機(jī)語言的優(yōu)點(diǎn)，“它擁有類似于C語言一樣的執(zhí)行速度，擁有如同Ruby語言的動(dòng)態(tài)性，又有Matlab般熟悉的數(shù)學(xué)記號(hào)和線性代數(shù)運(yùn)算能力，兼具像Python般的通用性，又像R語言一樣擅長于統(tǒng)計(jì)分析，并有Perl般處理字符串的能力和shell等膠水語言的特點(diǎn)，并易于學(xué)習(xí)”。目前已有多所國際知名高校的數(shù)值計(jì)算或統(tǒng)計(jì)學(xué)課程結(jié)合Julia語言進(jìn)行講解，如斯坦福大學(xué)的“應(yīng)用矩陣方法（Introduction to Matrix Methods；課程代碼：EE103）”和麻省理工學(xué)院的“線性代數(shù)（Linear Algebra；課程代碼：18.06）”。愛荷華州立大學(xué)動(dòng)物科學(xué)系2015年5月在其開設(shè)的“家畜基因組預(yù)測（Genomic Prediction in Livestock）”短期課程中結(jié)合Julia語言進(jìn)行了講解。使用7種標(biāo)準(zhǔn)檢查程序，Julia語言的運(yùn)行速度接近于C及Fortran語言（見圖1），但其編寫數(shù)值計(jì)算程序的速度卻快得多。一般情況下，Julia語言運(yùn)行數(shù)值計(jì)算程序時(shí)的速度也接近于C++，是R語言速度的100倍，MATLAB語言的1 000倍。

注：設(shè)C語言的運(yùn)行時(shí)間為1.0；C和Fortran語言使用gcc 4.8.2進(jìn)行編譯；C、Fortran和Julia使用OpenBLAS v0.2.13；Python運(yùn)行rand_mat_stat和rand_mat_mul使用NumPy v1.9.2庫函數(shù)。

2 Julia語言在動(dòng)物育種中的應(yīng)用

2.1 分子血緣相關(guān)矩陣（A）及其逆矩陣計(jì)算 分子血緣相關(guān)矩陣或加性遺傳相關(guān)矩陣及其逆矩陣計(jì)算在動(dòng)物育種學(xué)教學(xué)及種畜的遺傳評(píng)估中有重要意義。使用Julia語言計(jì)算分子血緣相關(guān)矩陣需構(gòu)建“Pedigree”和“PedNode”類型，使用“mkPedigree”函數(shù)對(duì)系譜中的個(gè)體進(jìn)行排序，使子代位于親代之后。使用“Amatrix”函數(shù)計(jì)算分子血緣相關(guān)矩陣[2]，其輸出結(jié)果使用稀疏矩陣存儲(chǔ)，可使用“round（full（A），2）”轉(zhuǎn)化為保留2位小數(shù)的滿矩陣。以圖2所示系譜為例，其重新編碼的系譜如下：

2.2 近交系數(shù)計(jì)算 近交系數(shù)為形成合子的2個(gè)配子源于同一共同祖先的概率。由通徑系數(shù)原理可知個(gè)體X的近交系數(shù)即為形成X個(gè)體的兩個(gè)配子間的相關(guān)系數(shù)[4]，用一般用Fx表示。使用“Inbreeding”函數(shù)計(jì)算10個(gè)個(gè)體的近交系數(shù)為：

2.3 設(shè)計(jì)矩陣 統(tǒng)計(jì)學(xué)中，設(shè)計(jì)矩陣（Design matrix）的元素為解釋變量（Explanatory variable），如在方差分析中用指示變量（1和0）表示連續(xù)變量在模型中的位置[5]。用“design”函數(shù)可以構(gòu)建設(shè)計(jì)矩陣，如5頭奶牛分別養(yǎng)殖在3個(gè)牧場（1，1，2，3，2），則由“design”函數(shù)構(gòu)建的設(shè)計(jì)矩陣以稀疏矩陣形式存儲(chǔ)，轉(zhuǎn)換滿矩陣為：

2.5 混合模型方程組（MME）建立 1953年，C.R.Henderson以混合模型為基礎(chǔ)建立線性方程組。由此方程組求解，可得到固定效應(yīng)的最佳線性無偏估計(jì)和隨機(jī)效應(yīng)的最佳線性無偏預(yù)測[7]?；旌夏Ｐ头匠探M是解決許多動(dòng)物育種學(xué)問題的基礎(chǔ)。由“MME”函數(shù)可以建立混合模型方程組的“左手項(xiàng)”和“右手項(xiàng)”，并得出相應(yīng)參數(shù)的估計(jì)量。

[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”函數(shù)轉(zhuǎn)化為設(shè)計(jì)矩陣；GI為3×3單位矩陣；RI為6×6單位矩陣；y=[2.0，1.5，2.0，1.2，0.89，1.2]，由“MME”函數(shù)計(jì)算得到，剩余方差（SSR）為13.121 8；參數(shù)[b]（[u]），“右手項(xiàng)（RHS）”，“左手項(xiàng)（LHS）”分別為：

2.6 基因組關(guān)系矩陣（G）的建立 “computeG”函數(shù)可以計(jì)算VanRaden（2008）或Yang等（2010）定義的G矩陣，并可以使用等位基因頻率平均值或2倍基因頻率（需滿足Hardy-Weinberg平衡）進(jìn)行中心化[8]。

3 結(jié)論

實(shí)踐中，由于R語言開發(fā)育種程序的速度較快，但程序運(yùn)行速度較慢，其編寫的程序很難用于實(shí)際育種工作。筆者只是用R語言測試統(tǒng)計(jì)方法的可行性，在此基礎(chǔ)上再將R程序轉(zhuǎn)變?yōu)镕ortran或C++等程序，整個(gè)開發(fā)過程實(shí)際經(jīng)歷了2個(gè)內(nèi)容基本相同的階段。Julia語言兼具有R語言和Fortran（或C++）的優(yōu)點(diǎn)，有效地縮短了程序開發(fā)時(shí)間。由Julia語言編寫的7個(gè)動(dòng)物育種中常用的程序（A矩陣計(jì)算、A逆矩陣計(jì)算、近交系數(shù)計(jì)算、設(shè)計(jì)矩陣、塊矩陣和、MME矩陣、G矩陣計(jì)算），可用于動(dòng)物育種學(xué)應(yīng)用程序的編寫和教學(xué)，源代碼可以通過郵件（meibujun@163.com或meibujun@iastate.edu）向作者索取。經(jīng)過測試，7個(gè)程序的計(jì)算結(jié)果可靠。這些程序可以作為廣大動(dòng)物育種工作者掌握J(rèn)ulia語言應(yīng)用于動(dòng)物育種學(xué)的基礎(chǔ)。

致謝：本研究部分靈感及部分計(jì)算設(shè)備由中國農(nóng)業(yè)大學(xué)動(dòng)物科技學(xué)院張勤教授課題組提供；感謝美國愛荷華州立大學(xué)動(dòng)物科學(xué)系Rohan L.Fernando教授、Hao Chen博士和Jian Zeng博士在編寫程序過程中的幫助。

參考文獻(xiàn)

[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. （責(zé)編：張宏民）