水產(chǎn)動(dòng)物基因組近交分析軟件的開(kāi)發(fā)和應(yīng)用

欒培賢，張曉峰，戶國(guó)

（中國(guó)水產(chǎn)科學(xué)研究院黑龍江水產(chǎn)研究所，黑龍江 哈爾濱 150070）

水產(chǎn)動(dòng)物遺傳資源是漁業(yè)生產(chǎn)和可持續(xù)發(fā)展的基礎(chǔ)，也是滿足未來(lái)不可預(yù)見(jiàn)需求的重要基因庫(kù)。因此，世界主要國(guó)家均極為重視水產(chǎn)動(dòng)物遺傳資源的保護(hù)。在制定保護(hù)計(jì)劃過(guò)程中，為使有限的人、財(cái)、物力發(fā)揮更好的效果，需要對(duì)保護(hù)的資源進(jìn)行選擇，區(qū)別需要保護(hù)的緊迫性和重要性，客觀地評(píng)價(jià)物種及地理居群受威脅的狀況[1]。受自然環(huán)境變遷及水利工程、過(guò)度捕撈、外來(lái)物種入侵等人類活動(dòng)影響，我國(guó)東北、西北、西南等廣大內(nèi)陸地區(qū)自然水域以魚(yú)類為代表的水產(chǎn)動(dòng)物遺傳資源都趨于衰退[2]。

一般說(shuō)來(lái)，水產(chǎn)動(dòng)物自然群體或小規(guī)模保種群體會(huì)發(fā)生小群體近交、個(gè)體生殖策略改變、個(gè)體小型化等情況，導(dǎo)致基因組某些位點(diǎn)的等位基因頻率改變或者丟失，非等位基因組合發(fā)生顯著變化[3]。這些基因組核苷酸水平的微小變化會(huì)體現(xiàn)在個(gè)體間親緣系數(shù)改變，體質(zhì)量、體型以及適應(yīng)性等重要經(jīng)濟(jì)性狀的加性效應(yīng)和非加性效應(yīng)方差組分變化[4]。有鑒于此，研發(fā)高效易行的水產(chǎn)動(dòng)物基因組近交分析方法并在水生野生動(dòng)物保護(hù)實(shí)踐中應(yīng)用，已成為保護(hù)遺傳學(xué)和水產(chǎn)生物技術(shù)領(lǐng)域共同關(guān)心的重要問(wèn)題之一。

本研究用高通量的SNP 分子標(biāo)記替代傳統(tǒng)家養(yǎng)群體中的系譜記錄信息，擬開(kāi)發(fā)精準(zhǔn)快速的軟件工具包，在無(wú)系譜信息魚(yú)類的野生自然群體或捕獲群體中構(gòu)建實(shí)現(xiàn)分子親緣關(guān)系矩陣，以多種分析方法解析并呈現(xiàn)群體遺傳特征，提供如基因組共祖系數(shù)（Genomic coancestry coefficient）、顯性親緣關(guān)系（Dominance coefficient）、血緣同源（Identity By Descent，IBD）、狀態(tài)同源（Identical By State，IBS）、基因組近交系數(shù)（Genomic inbreeding coefficient）、多維標(biāo)度分析（Multidimensional Scaling，MDS）及群體遺傳結(jié)構(gòu)等統(tǒng)計(jì)基因組分析領(lǐng)域常見(jiàn)統(tǒng)計(jì)參數(shù)的估計(jì)值。本軟件可以應(yīng)用于水產(chǎn)動(dòng)物遺傳資源受威脅程度及保種群體規(guī)模需求量化分析方法研究，并為評(píng)價(jià)內(nèi)陸珍稀水產(chǎn)動(dòng)物小群體保種、親本遺傳管理及種群修復(fù)效果提供新的技術(shù)途徑和科學(xué)可靠的參考依據(jù)。

1 軟件的設(shè)計(jì)開(kāi)發(fā)

1.1 軟件設(shè)計(jì)流程

本軟件以試驗(yàn)群體全部個(gè)體的全基因組范圍內(nèi)SNP 標(biāo)記的分型數(shù)據(jù)為數(shù)據(jù)流起點(diǎn)；經(jīng)過(guò)核心處理過(guò)程對(duì)輸入數(shù)據(jù)進(jìn)行全面的預(yù)處理，如樣本分組、SNP 位點(diǎn)篩選、數(shù)據(jù)格式轉(zhuǎn)換，將預(yù)處理后得到的最終數(shù)據(jù)集作為模型估計(jì)的輸入信息；根據(jù)所框選的統(tǒng)計(jì)模型和分析算法完成基因組SNP 標(biāo)記的各類親緣關(guān)系參數(shù)的估計(jì)，最終以多種分析方法解析并呈現(xiàn)群體遺傳特征，數(shù)據(jù)流和功能結(jié)構(gòu)如圖1所示。

圖1 技術(shù)路線Fig.1 Technical route

1.2 軟件結(jié)構(gòu)設(shè)計(jì)

本軟件采用數(shù)據(jù)處理層、模型工具層和可視化交互層的C/S 三層交互性構(gòu)架結(jié)構(gòu)設(shè)計(jì)（圖2）。數(shù)據(jù)處理層包含SNP 基因分型數(shù)據(jù)集、個(gè)體和分組等數(shù)據(jù)的預(yù)處理；模型工具層包含集成的開(kāi)源分析工具（內(nèi)嵌了R、Perl 等軟件及工作環(huán)境），可作為數(shù)據(jù)處理、模型建立、參數(shù)估計(jì)和計(jì)算結(jié)果統(tǒng)計(jì)分析的工作環(huán)境；可視化層包含工具包可視化交互環(huán)境，以及親緣關(guān)系和群體遺傳結(jié)構(gòu)分析的可視化呈現(xiàn)。

圖2 軟件開(kāi)發(fā)的三層交互結(jié)構(gòu)Fig.2 Three layers interactive structure in software development

1.3 開(kāi)發(fā)工具與運(yùn)行環(huán)境

本軟件是在Windows 操作系統(tǒng)下，基于VB.NET 2012 語(yǔ)言、開(kāi)源工具包（R、PLINK[5]和GVCBLUP[6]、LEA[7]、SnpSift[8]）開(kāi)發(fā)完成，可在安裝有Microsoft.NET Framework 4.0 運(yùn)行環(huán)境的Windows 10 操作系統(tǒng)上運(yùn)行。為使本軟運(yùn)行穩(wěn)定流暢，操作系統(tǒng)的硬件配置不低于酷睿i3 處理器、主頻2.13 GHz、4 GB內(nèi)存和128 GB 硬盤。

2 軟件的功能實(shí)現(xiàn)

本軟件通過(guò)原始數(shù)據(jù)提交、數(shù)據(jù)預(yù)處理、遺傳分析等3 個(gè)功能模塊實(shí)現(xiàn)其功能，具體情況如下：

2.1 原始數(shù)據(jù)提交模塊

工具包使用者通過(guò)原始數(shù)據(jù)提交模塊，選擇需要分析的原始數(shù)據(jù)（VCF 格式文件）。原始數(shù)據(jù)提交后，工具包會(huì)自動(dòng)解析原始數(shù)據(jù)信息，給出所提交數(shù)據(jù)集中樣本數(shù)量和SNP 標(biāo)記數(shù)量信息等相關(guān)信息，如圖3 所示。

圖3 原始數(shù)據(jù)提交模塊Fig.3 VCF file submission module

2.2 數(shù)據(jù)預(yù)處理模塊

該模塊是數(shù)據(jù)分析前的重要預(yù)處理環(huán)節(jié)，用來(lái)對(duì)原始數(shù)據(jù)進(jìn)行數(shù)據(jù)質(zhì)量控制、SNP 標(biāo)記篩選和樣本標(biāo)簽的管理（圖4）。選擇有效的約束條件和控制標(biāo)準(zhǔn)，可以降低數(shù)據(jù)背景噪聲，提高結(jié)果的準(zhǔn)確性和可讀性。使用者可以通過(guò)多種技術(shù)手段高通量獲得全基因組范圍內(nèi)SNP 基因分型信息，如基因組重測(cè)序、簡(jiǎn)化基因組測(cè)序和基因芯片等。由于建庫(kù)方式、測(cè)序量等不同，使全基因組范圍內(nèi)各個(gè)SNP 座位的數(shù)據(jù)完整性不同。該模塊能夠?qū)υ紨?shù)據(jù)進(jìn)行質(zhì)量控制，設(shè)定次要等位基因頻率、哈溫平衡檢測(cè)顯著水平、缺失基因型比例等閾值，使每個(gè)樣本用于進(jìn)一步分析的SNP 標(biāo)記具有較好的一致性。

圖4 數(shù)據(jù)預(yù)處理模塊Fig.4 Data preprocessing module

基因分型文件內(nèi)的樣本標(biāo)簽多為測(cè)序分析公司內(nèi)部編號(hào)。為了滿足分析的實(shí)際要求，需將無(wú)意義的測(cè)序編號(hào)對(duì)應(yīng)有序的一組個(gè)體編號(hào)上。此外，進(jìn)行群體遺傳結(jié)果分析時(shí)，也需要為個(gè)體所在群體指定群體編號(hào)。使用者需要將基因分型的個(gè)體編號(hào)與分析編號(hào)、群體編號(hào)相互對(duì)應(yīng)，保存在文檔中進(jìn)行上傳。

2.3 遺傳分析模塊

該模塊是本軟件的核心模塊，集成了多種開(kāi)源工具包，估計(jì)個(gè)體近交水平和個(gè)體間基因組親緣關(guān)系，分析群體遺傳結(jié)構(gòu)等。使用者只需點(diǎn)選所需分析內(nèi)容（圖5），軟件工具通過(guò)后臺(tái)調(diào)用集成在本軟件內(nèi)的分析工具，即可在后臺(tái)完成全部分析過(guò)程，輸出相應(yīng)結(jié)果。

圖5 遺傳分析模塊Fig.5 Genetic analysis module

2.3.1 個(gè)體近交分析

該子模塊基于全基因組范圍內(nèi)的SNP 標(biāo)記估計(jì)個(gè)體近交水平，使用者通過(guò)個(gè)體分組信息，可進(jìn)一步獲取不同群體的近交水平等信息。該工具包提供兩種個(gè)體基因組近交系數(shù)的估算方法：一是經(jīng)典數(shù)量遺傳學(xué)理論基于個(gè)體加性關(guān)系矩陣對(duì)角線元素計(jì)算個(gè)體近交系數(shù)；二是群體遺傳學(xué)基于觀測(cè)和期望純合子關(guān)系評(píng)估個(gè)體近交系數(shù)。

2.3.2 親緣關(guān)系分析

該子模塊提供多種基因組親緣關(guān)系參數(shù)估計(jì)方法，如共祖系數(shù)、顯性親緣關(guān)系、血緣同源（IBD）、狀態(tài)同源（IBS）等，可從多個(gè)維度解析樣本個(gè)體間的親緣關(guān)系構(gòu)成。使用者可以根據(jù)分析需求，選取關(guān)注的基因組親緣關(guān)系參數(shù)，同時(shí)還能選擇參數(shù)估計(jì)的模型。

2.3.3 群體遺傳分析模塊

該子模塊基于高通量SNPs 標(biāo)記計(jì)算分析常見(jiàn)群體遺傳結(jié)構(gòu)參數(shù)。多維標(biāo)度分析（MDS）是一種類似于主成分分析（Principal components analysis，PCA）的降維分析方法?？稍诙S或三維空間展示個(gè)體間的基因組相似性，呈現(xiàn)出基于基因組相似性的個(gè)體空間散布圖，展示群體結(jié)構(gòu)。本軟件還集成了LEA 工具，利用非連鎖的SNP 標(biāo)記估計(jì)群體遺傳組成單元，推測(cè)現(xiàn)有群體中個(gè)體的原始來(lái)源。

3 實(shí)際應(yīng)用案例

采用本軟件，以利用SLAF-seq 技術(shù)分子標(biāo)記開(kāi)發(fā)的SNP 數(shù)據(jù)集為示例，按如下流程分析了新疆兩個(gè)白斑狗魚(yú)Esox lucius 群體共計(jì)68 尾個(gè)體（圖6）。該數(shù)據(jù)集將酶切片段長(zhǎng)度在414～464 的序列定義為SLAF 標(biāo)簽，預(yù)測(cè)到138 626 個(gè)SLAF 標(biāo)簽。實(shí)驗(yàn)中選擇RsaI+HaeIII 的酶切方案，共得到189.71 Mreads，獲得多態(tài)性的SLAF 標(biāo)簽313 663 個(gè)。

圖6 野生魚(yú)類親緣關(guān)系和近交水平分析流程Fig.6 Analysis flowing chart of genomic relatedness and inbreeding of wild fishes

在數(shù)據(jù)預(yù)處理模塊中，用戶可以根據(jù)自身對(duì)測(cè)序數(shù)據(jù)集質(zhì)量控制的需要設(shè)置相應(yīng)參數(shù)，獲得相應(yīng)的數(shù)據(jù)集用于進(jìn)一步遺傳分析。本研究相關(guān)參數(shù)設(shè)置為測(cè)序質(zhì)量值QUAL 不低于30，次要等位基因頻率不低于0.05，位點(diǎn)哈溫平衡顯著水平設(shè)置為0.01，每個(gè)SNP 位點(diǎn)的完整性為95%，樣本編號(hào)信息采用自定義編號(hào)。原始數(shù)據(jù)進(jìn)行預(yù)處理之后，得到14 124 個(gè)位點(diǎn)的高質(zhì)量數(shù)據(jù)。在遺傳分析模塊中，直接點(diǎn)選全部四種基因組親緣關(guān)系分析方法，并勾選定義6 估計(jì)共祖系數(shù)和顯性親緣關(guān)系；選擇PLINK 估計(jì)方法，估計(jì)個(gè)體的基因組近交系數(shù)。在群體遺傳分析模塊中，選擇默認(rèn)參數(shù)，用LEA 和MDS 兩種分析方法顯示群體遺傳結(jié)構(gòu)。

可以直接下載基因組近交系數(shù)估計(jì)值，也可將個(gè)體的近交系數(shù)估計(jì)值以柱形圖的形式展示（圖7）。

圖7 近交系數(shù)柱形圖Fig.7 Histogram of genomic inbreeding coefficient

共祖系數(shù)、顯性親緣關(guān)系、血緣同源、狀態(tài)同源等基因組親緣關(guān)系參數(shù)均以熱圖方式呈現(xiàn)。本文中以共祖系數(shù)為例，以熱圖形式展示相關(guān)結(jié)果（圖8）。

圖8 共祖系數(shù)熱圖Fig.8 Heatmap of genomic coancestry coefficient

基于LEA 推測(cè)群體構(gòu)成結(jié)構(gòu)（圖9），用MDS方法展示個(gè)體空間散布圖（圖10），兩種群體遺傳結(jié)構(gòu)分析的常用方法也可以直接以圖形的方式從軟件輸出。

圖9 個(gè)體遺傳來(lái)源估計(jì)結(jié)果展示圖Fig.9 Histogram of individual ancestral genetic source estimation

圖10 二維MDS 散點(diǎn)圖Fig.10 Two dimensional plot of MDS

4 分析與討論

目前，增殖放流已成為水生生物資源養(yǎng)護(hù)的一項(xiàng)重要措施之一[9]。但是，野生自然群體沒(méi)有遺傳結(jié)構(gòu)清晰的系譜記錄信息，傳統(tǒng)的標(biāo)志放流方法無(wú)法準(zhǔn)確評(píng)估有效群體規(guī)模、近交率等遺傳滅絕風(fēng)險(xiǎn)指標(biāo)，無(wú)法精確量化評(píng)價(jià)遺傳資源瀕危程度[10]。一些研究者和社會(huì)觀察家認(rèn)為，如果僅靠單純數(shù)量意義上的增殖放流不僅會(huì)在經(jīng)濟(jì)、環(huán)境和增殖放流效果等方面不可持續(xù)，而且會(huì)降低野生群體的遺傳多樣性，使得種群遺傳結(jié)構(gòu)趨于簡(jiǎn)單化[11]。同樣，野生自然群體沒(méi)有遺傳結(jié)構(gòu)清晰的系譜記錄信息，無(wú)法采用傳統(tǒng)家養(yǎng)動(dòng)物的遺傳評(píng)定手段[12]，有時(shí)甚至?xí)霈F(xiàn)放流種群非土著居群后代，占據(jù)了土著種群的生存空間和生態(tài)位，替代本地野生群體；或者放流群體與本地土著群體發(fā)生遺傳交流，改變本地野生群體的遺傳組成，破壞本地魚(yú)類遺傳資源。

現(xiàn)階段，人們可以方便快捷并廉價(jià)地獲得基于二代、三代測(cè)序技術(shù)的高通量SNP 分子遺傳標(biāo)記，應(yīng)用這些標(biāo)記信息就可以準(zhǔn)確地估計(jì)野生水生動(dòng)物的基因組近交水平。在開(kāi)展?jié)O業(yè)資源養(yǎng)護(hù)過(guò)程中，必須考慮水產(chǎn)動(dòng)物基因組近交水平等因素。為此，本軟件為一線野生水生動(dòng)物保護(hù)工作者集成了現(xiàn)有依托于不同環(huán)境如R、Perl 及常用獨(dú)立群體基因組學(xué)軟件包的功能，提供了主流群體基因組和保護(hù)遺傳學(xué)常用基因組參數(shù)估計(jì)，并利用VB.Net 予以界面可視化，實(shí)現(xiàn)全程鼠標(biāo)點(diǎn)選即可實(shí)現(xiàn)高難度的群體基因組學(xué)分析。

遺傳模擬分析中，本軟件集成的各類方法理論上適用各種二倍體有性生殖動(dòng)物。本研究還選取了新疆重要的土著經(jīng)濟(jì)魚(yú)類白斑狗魚(yú)群體的SLAF 簡(jiǎn)化基因組測(cè)序數(shù)據(jù)集作為本研究的實(shí)際案例，用于確認(rèn)軟件功能的有效性，也取得了實(shí)際效果。盡管如此，在樣本量較大或受成本制約等情況下，有時(shí)仍會(huì)考慮使用SSR 等傳統(tǒng)分子標(biāo)記。因此，本研究也拓展了可以兼容SSR 標(biāo)記的分子標(biāo)記共祖分析方法[13，14]，將在下一階段軟件升級(jí)時(shí)集成整合。該軟件的推廣和應(yīng)用有望提升漁業(yè)領(lǐng)域的保種技術(shù)水平,推動(dòng)我國(guó)水產(chǎn)動(dòng)物種質(zhì)資源保存及管理工作。

綜上所述，本軟件包界面簡(jiǎn)約友好，功能較為全面，預(yù)期可以應(yīng)用在水產(chǎn)動(dòng)物種質(zhì)資源保存管理、土著魚(yú)類野生種群修復(fù)效果評(píng)價(jià)等很多領(lǐng)域。本軟件提供的基因組近交分析工具可以協(xié)助建立基于水產(chǎn)動(dòng)物遺傳資源瀕危程度評(píng)價(jià)模型，計(jì)算水產(chǎn)動(dòng)物同一群體內(nèi)或不同群體間個(gè)體親緣關(guān)系、近交系數(shù)、保種所需有效群體規(guī)模等，并在此基礎(chǔ)上，設(shè)計(jì)小群體保種的合理交配策略。