陳梅麗，馬英克，李茹姣，鮑一明

1.中國科學院北京基因組研究所（國家生物信息中心），國家基因組科學數(shù)據(jù)中心和中國科學院基因組科學與信息重點實驗室，北京 100101

2.中國科學院大學，未來技術(shù)學院，北京 100049

引言

隨著人類基因組測序計劃的完成，基因組學的影響力迅速擴大，數(shù)以萬計的動物、植物、微生物基因組被組裝[1-4]，游離DNA 應用于無創(chuàng)產(chǎn)前檢測、通過基因檢測指導靶向藥物治療成為可能[5-6]、單細胞技術(shù)應用于輔助生殖[7]、DNA 編輯技術(shù)廣泛應用[8]、體外提供造血干細胞[9]、長壽基因被找到[10]、抗病蟲害和高產(chǎn)的農(nóng)業(yè)作物新品種被不斷培育[11-12]。這些基因組科學領(lǐng)域的巨大進展，一方面來自于測序和實驗技術(shù)的革新，同時也依賴于為適應實驗和測序技術(shù)進步而不斷發(fā)展的分析手段和方法[4,13-14]。

基因組學測序數(shù)據(jù)增長迅猛，加快數(shù)據(jù)分析速度，提高數(shù)據(jù)處理效率，是對大數(shù)據(jù)整合分析工具和算法開發(fā)的迫切需求。如何用好多源海量基因組學數(shù)據(jù)，去除異質(zhì)性，并對其進行整合分析和深度挖掘，是對數(shù)據(jù)分析工具和算法開發(fā)另一個層面的新要求。科學家們也在不斷地開發(fā)各種算法和工具來提高計算效率，比如第三代測序數(shù)據(jù)組裝算法wtdbg2，將拼接的分析速度提高5倍，少于數(shù)據(jù)產(chǎn)出時間[4]。而在此基礎上，人工智能等先進技術(shù)也被廣泛應用于基因組學大數(shù)據(jù)分析的工具開發(fā)[15]。本文將詳細地闡述隨著測序技術(shù)的發(fā)展，基因組、轉(zhuǎn)錄組、表觀組數(shù)據(jù)的分析算法和工具開發(fā)的現(xiàn)狀，以及大數(shù)據(jù)時代基因組學數(shù)據(jù)算法和工具開發(fā)在未來將面臨的問題和挑戰(zhàn)。

1 基因組測序數(shù)據(jù)分析

隨著測序技術(shù)的發(fā)展，各類基因組相關(guān)研究計劃接踵而來[16-17]，為生物多樣性、物種進化、分子育種、臨床治療等研究提供了寶貴的數(shù)據(jù)資源?；蚪M數(shù)據(jù)分析主要包括基因組組裝、基因組注釋、基因組變異分析等?；蚪M序列和注釋信息包含了生物體的所有遺傳信息和功能信息，是多種組學研究的重要基礎數(shù)據(jù)[18]。通過基因組變異分析可以解析基因組變異對表型、物種進化、疾病等的影響[19]。但是測序數(shù)據(jù)分析時間往往遠大于測序數(shù)據(jù)產(chǎn)出時間，無法匹配數(shù)據(jù)爆發(fā)式增長的趨勢，是基因組數(shù)據(jù)分析面臨的最大挑戰(zhàn)之一。

1.1 基因組組裝

基因組組裝是將測序產(chǎn)生的讀段（read）片段經(jīng)過序列組裝成完整的基因組序列?；蚪M組裝方法主要有兩類：（1）基于參考基因組序列比對的有參組裝，常用于重測序和線粒體/葉綠體等保守細胞器基因組[20-21]；（2）從頭（de novo）組裝，目前主要有純二代測序組裝[22]、二代和三代測序混合組裝[23]、純?nèi)鷾y序組裝[24]等組裝策略。對于動植物等大基因組可以利用遺傳圖譜、Bionano 光學圖譜、Hi-C、10X Genomics 等圖譜信息進行整合輔助組裝，將基因組組裝提升到染色體級別[25-28]。動植物基因組因存在基因組大、雜合度高、GC 含量高、重復序列多和多倍化水平較高等復雜因素，給組裝帶來了很大的挑戰(zhàn)，需要開發(fā)出兼顧效率和計算資源消耗，并且在重復區(qū)域獲得連續(xù)性和完整性都表現(xiàn)很好的高質(zhì)量基因組組裝結(jié)果的算法和工具。目前常用的軟件有SOAPdenovo2[22]、ALLPATHS-LG[29]、Canu[30]、Falcon[31]等。

為了解決大型基因組組裝效率、準確性要求高和計算資源消耗大的問題，阮玨團隊提出了三代數(shù)據(jù)組裝wtdbg2 算法，該算法遵循overlap-layoutconsensus 模式，以快速的讀段“全部對全部”比對實現(xiàn)方式和基于模糊布魯因圖這種新組裝圖理論——一種與稀疏布魯因圖和其變體A-Bruijn 圖有關(guān)的序列組裝的新數(shù)據(jù)結(jié)構(gòu)，來改進現(xiàn)有的組裝程序并提高組裝效率[4]。在一臺計算機上，可在2 天內(nèi)完成4個～30X的人類基因組數(shù)據(jù)集的組裝，極大提高三代測序數(shù)據(jù)的分析效率。與此前提出的Flye 算法[32]相比，wtdbg2的分析速度提升了5倍，內(nèi)存使用僅為Flye的一半，組裝連續(xù)性和精度可與Canu、Falcon、Flye 等其他算法相媲美。wtdbg2 首次將測序數(shù)據(jù)分析時間降低到少于測序數(shù)據(jù)產(chǎn)出時間，是一種兼具高容錯和高準確的高效算法，并可擴展到超大基因組，如32 Gb 大小的蠑螈基因組[4]。

當前Hi-C 技術(shù)越來越多地應用于輔助染色體水平二倍體基因組組裝[25]，Hi-C 技術(shù)是基于將線性距離遠、空間結(jié)構(gòu)近的DNA 片段進行交聯(lián)，并將交聯(lián)的DNA 片段富集后進行高通量測序，對測序數(shù)據(jù)分析揭示染色質(zhì)的遠程相互作用，明確基因組草圖中scaffold和染色體對應關(guān)系、scaffold之間的連接順序和方向，從而將scaffold 掛載到染色體上。但是對于同源多倍體和近期加倍的異源多倍體來說，其同源染色體之間的Hi-C 交聯(lián)信號會將序列相似的等位片段連接在一起，導致同源染色體被錯誤地連接到一起，形成大量嵌合的組裝，給組裝造成了較大困難。針對該問題明瑞光團隊提出了ALLHiC 算法，該算法包括pruning、partition、rescue、optimization、building 5個步驟。ALLHiC 算法一方面通過修剪Hi-C 平行信號和弱信號進行等位基因分型，減少了同源染色體間的嵌合連接；另一方面通過遺傳算法隨機優(yōu)化，極大地提高了contig 序列的排序和定向準確性，成功解決了同源高倍體甘蔗、菠蘿和異源四倍體栽培花生等多倍體組裝難題[33-34]。

而針對于二倍體或者是多倍體，單體型基因組組裝是基因組組裝的最終目標，目前已提出單體型區(qū)塊組裝策略，如FALCON-Unzip[31]、trio binning[35]等方法，嘗試將兩個或多個配子來源的染色體進行分別組裝。但是當前已有的單體型區(qū)塊組裝策略依然無法很好地解決單條染色體長度的單體型基因組組裝問題。針對該問題，Kronenberg 團隊提出了一種單體型基因組組裝的新技術(shù)FALCON-Phase，可把雜合基因組的父本母本分型組裝，并可將其應用于野外采集的樣本或缺乏譜系信息的生物。其主要原理和流程為：將基因組區(qū)域中顯示出高水平雜合性的區(qū)域鑒定為單體型區(qū)塊contig，基于鑒定結(jié)果對所有contig 拆分和打斷，通過將Hi-C 數(shù)據(jù)集成到圖形數(shù)據(jù)結(jié)構(gòu)中構(gòu)建標準化的互作矩陣，最后經(jīng)過Hi-C 輔助組裝掛載獲得單條染色體長度的單體型基因組組裝[36]。

優(yōu)化算法解決非橋式重復等重復序列導致的組裝連續(xù)性問題[32]，通過開發(fā)出對雜合子感知的一致性算法和分型組裝方法，兼具高容錯、高準確、高效、計算資源低耗的優(yōu)點，解決超大基因組組裝和大規(guī)?；蚪M組裝的計算困境，是未來基因組組裝工具開發(fā)需要繼續(xù)改進優(yōu)化的方向。同時高效低耗的組裝算法也能解決真核生物泛基因組（pan-genome）構(gòu)建的困境[37]。在單體型基因組組裝方面，未來的發(fā)展可以優(yōu)化算法將單體型組裝擴大到高雜合性樣品、更高倍性物種的應用上，將分型算法整合到組裝中，并通過基于組裝圖的方法提高單體型區(qū)塊contig 放置的準確率，提高算法性能[36]。

1.2 基因組注釋

基因組注釋是對組裝的基因組進行基因結(jié)構(gòu)和功能注釋?；蝾A測方法依賴于利用兩種類型的基因信息：（1）內(nèi)容——局部位點，如剪接位點、起始密碼子、終止密碼子等，預測編碼和非編碼區(qū)域；（2）信號——蛋白質(zhì)功能位點，檢測功能性位點是否存在。原核生物的基因結(jié)構(gòu)簡單，各種局部位點特異性強、易識別，基因預測方法基本成熟，如GeneMarkS[38]、Glimmer[39]等。而真核生物的基因存在內(nèi)含子結(jié)構(gòu)，要經(jīng)過RNA 剪接步驟才能形成成熟的RNA，使得真核基因的預測更為復雜和困難，需要采用更加復雜的算法來預測真核生物的基因結(jié)構(gòu)[18]。如何準確地識別可變剪切位點是真核基因預測的技術(shù)難點，目前有三種策略來預測真核基因結(jié)構(gòu)：（1）同源預測：有一些基因蛋白在相近物種間是高度保守的，可以使用已有的高質(zhì)量近緣物種注釋信息通過比對的方式確定外顯子邊界和剪切位點，該方法最可靠，但是對數(shù)據(jù)的依賴性很高，所選取物種間的親緣關(guān)系和進化距離對結(jié)果影響較大[40]；（2）從頭預測：通過已有的概率模型（如隱馬爾科夫模型）來預測基因結(jié)構(gòu)，再預測剪切位點和非翻譯區(qū)，這類方法需要預先訓練獲得預測參數(shù)，準確性較低，尤其是對于研究較少的物種，這種情況建議利用近緣物種的同源基因數(shù)據(jù)以訓練基因結(jié)構(gòu)預測參數(shù)[18]；（3）基于轉(zhuǎn)錄組預測：通過物種的RNASeq/EST 等轉(zhuǎn)錄組數(shù)據(jù)輔助注釋，能夠較為準確地確定剪切位點和外顯子區(qū)域[41]，但是低表達基因和轉(zhuǎn)錄組組裝錯誤也是該方法的主要瓶頸。可以開發(fā)出綜合上述三種策略的算法，兼具準確性和特異性地識別可變剪切位點?；诳煽康幕蚪Y(jié)構(gòu)注釋，后續(xù)可以進行基因功能注釋，目前普遍采用比對的方法[40,42]，但是序列相似與實際生物學功能相似無法完全等價，需考慮引入其它方法，如引入基因表達、調(diào)控網(wǎng)絡等信息，開發(fā)算法進行信息整合進一步完善基因功能注釋工作。

1.3 基因組變異分析

基因組變異是指與參考序列相比，基因組中發(fā)生的單堿基變異、DNA 序列片段插入、缺失、擴增和復雜結(jié)構(gòu)變異等。目前基于測序方法進行單核苷酸 多 態(tài) 性（Single Nucleotide Polymorphism,SNP）、短的插入缺失（Insertion-Deletion,InDel）等變異檢測的策略主要有：（1）基于比對結(jié)果直接檢測變異信息；（2）基于從頭組裝結(jié)果比對檢測變異信息。而在結(jié)構(gòu)變異（Structural Variation,SV）檢測方面目前主要有四種策略：（1）基于插入片段長度的異常分布（Read Pair,RP）；（2）基于讀段覆蓋度的異常分布（Read Depth,RD）；（3）基于未比對上的讀段進行分割比對（Split Read,SR）；（4）長讀段或者從頭組裝[43]。不同策略可檢測到的變異類型不同，可將多個策略整合在一起，從而大大降低假陽性率[43]。基于拼接的方法是當前獲得全基因組范圍內(nèi)各類型變異最好的方法，但是該方法的結(jié)果準確性取決于高質(zhì)量組裝，有待高效準確組裝方法的發(fā)展。GATK是當前公認的最主流的基因組變異分析軟件之一，可應用于人和其它物種[44]。但是當前基因組數(shù)據(jù)變異分析在計算處理能力的瓶頸限制了這些技術(shù)的廣泛使用，針對該問題已有用FPGA 進行硬件加速、搭載GPU 加速的Parabricks 軟件等方案來加速GATK的運行效率。測序讀長加長可以提高變異分析的準確度，尤其是大的結(jié)構(gòu)變異分析。未來測序技術(shù)的發(fā)展，兼具測序讀長長和測序準確性，也將極大的改變變異檢測方法的現(xiàn)狀[45]。

如何解讀這些變異信息，從海量的變異位點中篩選出真正具有生物學意義的基因位點是變異組學運用的主要瓶頸之一。目前可以通過GWAS 分析獲得基因型和表型之間的關(guān)系[46]。在GWAS 分析中需要考慮適合樣本的關(guān)聯(lián)模型。樣本的表型分兩種：（1）數(shù)量性狀的表型，這種關(guān)聯(lián)分析主要采用線性回歸模型，包括一般線性模型和混合線性模型，對于受到多因素共同影響的復雜數(shù)量性狀常采用混合線性模型或者改進后的模型進行分析，仔細選擇模型或者整合互補模型可以對復雜性狀的遺傳基礎有更深入的了解[47-48]；（2）case-control 二元的表型，case和control 比例失衡時，可能會導致較高的假陽性率，針對該問題提出的SAIGE模型顯示出在樣品比例極其失衡時結(jié)果仍較為正常的特點[49]。當前GWAS 分析對于低頻基因的挖掘能力有待提高，可以通過改進開發(fā)新的模型和新的群體設計揭示表型變異的確定原因[47]。而將GWAS 關(guān)聯(lián)結(jié)合轉(zhuǎn)錄組、甲基化組等組學研究的數(shù)據(jù)分析，即eQTL、meQTL 等分析，則可以將表型與基因組變異，再到基因表達、甲基化等之間的關(guān)系聯(lián)系起來。整合多組學數(shù)據(jù)進行聯(lián)合分析揭示復雜性狀的分子機制，這也是目前變異解析的主要挑戰(zhàn)[50]。在識別各類QTLs 位點的時候，需要進行多重假設檢驗和校正，常用的多重假設檢驗校正方法有Bonferroni 校正、BH 校正和隨機打亂后進行Storey 校正等[51]。當前的多重假設檢驗校正方法難以兼具準確性和速度，保證準確度需要耗費大量的內(nèi)存和時間，尤其是對于trans-eQTLs 位點識別時用隨機打亂后進行Storey 校正方法，通常的超級計算機集群是基本不可能完成的，亟需新的模型和算法優(yōu)化來攻克這道難關(guān)[52]。

2 轉(zhuǎn)錄組測序數(shù)據(jù)分析

轉(zhuǎn)錄組是細胞內(nèi)所有RNA 分子的總和。細胞中的DNA 記錄的遺傳信息通過轉(zhuǎn)錄表達，決定生物體的性狀。轉(zhuǎn)錄組中，mRNA是DNA和蛋白質(zhì)之間的中間產(chǎn)物，而非編碼RNA 則有多種多樣的功能。轉(zhuǎn)錄組技術(shù)可以捕捉到細胞內(nèi)所有RNA 分子瞬時的存在情況。轉(zhuǎn)錄組數(shù)據(jù)分析技術(shù)對轉(zhuǎn)錄組技術(shù)產(chǎn)生的數(shù)據(jù)進行分析和挖掘，得到的結(jié)果和發(fā)現(xiàn)對基礎研究、臨床診斷等都有很大的推動。比如Hammond 等[53]對小鼠發(fā)育過程中、老年和大腦受傷后的76 000個小神經(jīng)膠質(zhì)細胞進行了單細胞測序，發(fā)現(xiàn)年輕小鼠的小神經(jīng)膠質(zhì)細胞的異質(zhì)性最高，老年小鼠小神經(jīng)膠質(zhì)細胞中的炎性細胞最多，并鑒定出受傷后發(fā)生反應的小神經(jīng)膠質(zhì)細胞亞型，為今后鑒定和操作健康和疾病狀態(tài)的各種亞型的小神經(jīng)膠質(zhì)細胞奠定了基礎。而在臨床診斷中，轉(zhuǎn)錄組數(shù)據(jù)分析結(jié)果可以得出超出生理范圍的基因表達情況、等位基因的特異表達和異常的選擇性剪切等[54]，對基因組測序數(shù)據(jù)是一個重要的補充。

2.1 RNA-Seq 數(shù)據(jù)分析

RNA-Seq 技術(shù)是目前科研中最常用的轉(zhuǎn)錄組技術(shù)[55]。由于RNA-Seq 讀長短，一般無法測出全長的轉(zhuǎn)錄本，為了檢測哪些轉(zhuǎn)錄本在樣品中表達，在所研究物種有高質(zhì)量的參考基因組的情況下，通常會先把測序短讀段比對到參考基因組，再組裝成全長轉(zhuǎn)錄本[56]。在研究的物種沒有參考基因組或參考基因組質(zhì)量較差的情況下，需要對轉(zhuǎn)錄本進行從頭組裝[41]。如果物種的參考轉(zhuǎn)錄組注釋質(zhì)量高或使用從頭組裝的轉(zhuǎn)錄組，則可以將讀段比對到轉(zhuǎn)錄組上直接定量[57]。對基因表達的定量一般是通過直接計算比對到基因上的讀段數(shù)，然后對基因長度和測序深度做歸一化，得到基因表達的RPKM（Reads Per Kilobase per Million mapped reads），F(xiàn)PKM（Fragments Per Kilobase per Million mapped reads）或TPM（Transcripts Per Million mapped reads）值[56]。也有一些定量方法不需要將讀段比對到基因序列上，比如通過對讀段中k-mer 計數(shù)的方法獲得基因的表達量[58]。在對基因的表達數(shù)據(jù)做樣本間、樣本內(nèi)歸一化、去除批次效應和考慮其它可能產(chǎn)生偏差的因素后，可以基于對表達量統(tǒng)計分布的推斷（一般假設基因的表達量服從負二項分布或泊松分布），通過統(tǒng)計檢驗確定在分組間差異表達的基因[59]。對轉(zhuǎn)錄本水平上的差異表達分析有兩種策略：基于轉(zhuǎn)錄異構(gòu)體（isoform）的策略和基于外顯子的策略。基于轉(zhuǎn)錄異構(gòu)體的策略對轉(zhuǎn)錄異構(gòu)體的表達進行定量，然后檢測差異表達[56]；基于外顯子的策略通過對比對到外顯子上和跨剪切位點的讀段進行計數(shù)來檢測選擇性剪切的信號[60]。融合基因的鑒定類似于新轉(zhuǎn)錄異構(gòu)體的鑒定，但由于融合的兩個基因可能位于不同的染色體上，搜索空間更大，并且要考慮到讀段比對的錯誤和生物學上的重要性對找到的融合基因進行篩選[61]。對小RNA 測序數(shù)據(jù)的分析一般是把讀段直接比對到已知的小RNA 序列上進行定量，或者把讀段比對到參考基因組，提取上下游序列進行莖環(huán)結(jié)構(gòu)預測，預測miRNA 前體序列和可能的新miRNA[62]。

2.2 三代轉(zhuǎn)錄組測序數(shù)據(jù)分析

三代轉(zhuǎn)錄組測序數(shù)據(jù)的特點是讀長長、錯誤率高和通量較低。因其讀段長，大多數(shù)情況下可以直接得到全長的轉(zhuǎn)錄本序列，不需要組裝。但因為其錯誤率較高，需要在轉(zhuǎn)錄本鑒定前對讀段進行錯誤校正[63]。一般錯誤校正的方法分為兩類：（1）基于混合讀段的錯誤校正，即利用對同一個樣品的RNASeq 短讀段序列對易出錯的長讀段進行校正[64]；（2）基于讀段重疊區(qū)的錯誤校正，即通過比對不同讀段間重疊的序列來校正這段被多次測到的序列[65]。最新的三代轉(zhuǎn)錄組的建庫和測序策略（如PacBio 公司的ISO-Seq 技術(shù)[66]和Oxford Nanopore公司的R2C2 方法[67]），充分利用目前測序讀段長度遠長于轉(zhuǎn)錄本長度的特點，在一個長讀段中對同一轉(zhuǎn)錄本連續(xù)進行多次測序，將每次測到的全長轉(zhuǎn)錄本（稱為subread）合在一起比對得到轉(zhuǎn)錄本的一致序列（consensus sequence），有效排除了測序錯誤的影響。由于之前三代測序技術(shù)的通量較低，一般使用二代和三代混合測序數(shù)據(jù)進行轉(zhuǎn)錄本定量。但目前主流三代測序技術(shù)的通量都有很大提高，如PacBio Sequel II 測序儀的一個SMRT Cell的通量已達到8 百萬個長讀段，達到了可以對轉(zhuǎn)錄本進行定量的標準，此時只需要對相應的長讀段進行簡單計數(shù)就可以得到基因或轉(zhuǎn)錄本的TPM值。三代測序數(shù)據(jù)由于讀長長的優(yōu)勢，在轉(zhuǎn)錄異構(gòu)體鑒定[68]、融合基因轉(zhuǎn)錄本鑒定[69]、轉(zhuǎn)錄本單體型鑒定（transcript haplotype phasing）[70]、轉(zhuǎn)錄本重復序列鑒定等應用上無論是準確度還是分析方法的簡便性對比RNASeq 技術(shù)都有很大優(yōu)勢，加上三代的RNA 直測技術(shù)[71]的興起可以去除建庫時反轉(zhuǎn)錄步驟帶來的人為干擾，使轉(zhuǎn)錄本的鑒定和定量更準確，三代轉(zhuǎn)錄組測序技術(shù)是未來轉(zhuǎn)錄組技術(shù)進步的一個方向。

2.3 單細胞RNA-Seq 數(shù)據(jù)分析

單細胞RNA-Seq 技術(shù)由于把瞬時基因表達譜的刻畫提升到了單個細胞的精度，極大地推進了人們對生命系統(tǒng)的認識，促成了很多新的科學發(fā)現(xiàn)，因此是當下最熱門和發(fā)展最快的轉(zhuǎn)錄組技術(shù)。單細胞RNA-Seq 產(chǎn)生的讀段除了轉(zhuǎn)錄本的片段序列外，一般會帶有一段細胞特異的（實際上是井或者液滴特異的）條形碼和一段轉(zhuǎn)錄本特異的唯一分子標識（Unique Molecular Identifier，UMI），首先需要根據(jù)這兩段標記序列給讀段分類，得到一個計數(shù)矩陣（count matrix），每一行代表一個細胞，每一列是一個轉(zhuǎn)錄本在各個細胞中表達的讀段數(shù)[13]。在做基因表達分析之前，需要綜合考慮每個條形碼下的讀段數(shù)（計數(shù)深度，count depth）、每個條形碼下表達的基因數(shù)和每個條形碼下的線粒體基因的讀段數(shù)三個特征來對測序數(shù)據(jù)進行質(zhì)量控制，去除死細胞、細胞膜破裂的細胞和一個井/液滴中含有兩個或多個細胞的情況[72]。由于建庫和測序過程中化學反應的復雜性，即使對于完全相同的細胞，都可能獲得不同的計數(shù)深度，因此在比較細胞間基因表達差異之前，需要在細胞水平上對轉(zhuǎn)錄本表達計數(shù)進行歸一化。單細胞RNA-Seq 計數(shù)矩陣中由于零值較多，需要專門設計歸一化的方法。單細胞RNA-Seq 計數(shù)矩陣歸一化的方法分為線性和非線性兩種：線性方法對一個細胞內(nèi)的所有轉(zhuǎn)錄本使用全局統(tǒng)一的縮放因子[73]；而非線性方法則可以去掉其它的偏差，比如批次因素[74]。在歸一化之后，計數(shù)矩陣會進行l(wèi)og(x+1)變換。

為了減少下游分析的復雜度和減少數(shù)據(jù)中的噪聲，避免“維數(shù)災難”，需要對計數(shù)矩陣進行特征選擇和降維。特征選擇步驟選出對于下游分析信息量最大的基因，通常會選擇表達量波動最大的1 000-5 000個基因（Highly Variable Genes,HVGs），作為降維算法的輸入[75]。降維算法在盡量保持數(shù)據(jù)結(jié)構(gòu)的基礎上，進一步將高維空間的數(shù)據(jù)映射到低維空間。最常用的降維方法分為兩類：（1）線性降維方法，主要關(guān)注在映射后保持表達差異點的距離，如主成分分析法（Principal Component Analysis,PCA）和多維標度法（Multidimensional Scaling,MDS）；（2）非線性降維方法，主要關(guān)注在映射后保持表達模式相似的點足夠接近、保持數(shù)據(jù)的局部結(jié)構(gòu)，同時也兼顧保持全局結(jié)構(gòu)，比如曲線成分分析（Curvilinear Components Analysis,CCA）、t-SNE（t-Distributed Stochastic Neighbor Embedding）、UMAP（Uniform Manifold Approximation and Projection）、擴散映射（Diffusion maps）等。在數(shù)據(jù)降維后，可以實現(xiàn)對單細胞RNA-Seq 數(shù)據(jù)的可視化。

對經(jīng)上述步驟處理好的表達數(shù)據(jù)進行模式識別是單細胞RNA-Seq 數(shù)據(jù)分析中的核心步驟。單細胞RNA-Seq 常見的下游分析有聚類分析、聚類注釋、組成分析、軌跡推斷和差異表達分析等。聚類分析根據(jù)基因表達譜定義的細胞間的距離對細胞聚類。細胞間距離的度量方式包括歐幾里得距離、余弦相似度、基于相關(guān)性的距離度量和通過對每個數(shù)據(jù)集用高斯核學習得到的距離度量等，后三種方式由于具有標度不變性，考慮值之間的相對差異，對于建庫大小、細胞大小不同的細胞之間的比較更健壯。傳統(tǒng)上，聚類分析方法主要有：（1）k-均值法，其中最常用的是RaceID和SIMLR，這兩種方法針對識別罕見細胞類型做了特殊處理；（2）層次聚類法，常用的有CIDR和PCAReduce，分別針對測序深度低導致的dropout 事件和罕見細胞類型做了優(yōu)化；（3）社區(qū)發(fā)現(xiàn)法（community detection），該類方法的思路是把密集連接的點，而非距離近的點歸為一類，這類方法被認為針對大數(shù)據(jù)集有更快的速度和更好的聚類效果，其中基于k-近鄰圖[76]的Louvain算法[77]在被Seurat和SCANPY 工具包采用后，在單細胞RNA-Seq 數(shù)據(jù)聚類分析中得到了最廣泛的應用。近期針對單細胞RNA-Seq 數(shù)據(jù)，出現(xiàn)了一些新的聚類分析思路，比如基于概率模型的CellAssign、BAMM-SC、基于深度神經(jīng)網(wǎng)絡的GOAE、GONN、ACTINN、scScope和scDeepCluster，這些方法在一些數(shù)據(jù)集上顯示出具有更高的聚類準確度，未來有望得到更廣泛的應用。對聚類的注釋需要依靠Mouse Brain Atlas和Human Cell Atlas[78]這樣的參考數(shù)據(jù)庫，或者先找出研究的一類細胞與所有其它細胞差異表達的基因（稱為標記基因）與數(shù)據(jù)庫中細胞的標記基因?qū)φ眨蛘哂眠@類細胞的整個基因表達譜與數(shù)據(jù)庫中的細胞表達譜對照，推斷這類細胞的類型。對聚類進行注釋之后，可以對細胞組成進行統(tǒng)計建模，用統(tǒng)計檢驗對不同數(shù)據(jù)集間的細胞組成差異進行分析。測序樣品中的細胞如果處于連續(xù)變化的過程中，比如樣品中同時存在處于發(fā)育、疾病進展、對處理條件響應等過程中不同階段的細胞時，旨在尋找離散的分組的聚類分析技術(shù)無法準確刻畫樣品中細胞的組成結(jié)構(gòu)，而軌跡推斷或偽時序分析技術(shù)則可以對細胞變化的路徑進行解析，把每類細胞在進程中的先后順序表示為線性的、二叉的、復雜圖的、樹的或多叉的拓撲結(jié)構(gòu)[79]。在基因?qū)用嫔?，可以對細胞間差異表達的基因進行分析。根據(jù)測試，為普通RNA-Seq 設計的差異表達分析方法加上對基因進行加權(quán)的方法的效果要好于目前專門為單細胞RNA-Seq 設計的差異表達分析方法[80]。

轉(zhuǎn)錄組分析技術(shù)已成為科學研究、臨床應用中不可缺少的關(guān)鍵技術(shù)環(huán)節(jié)。一方面，現(xiàn)有轉(zhuǎn)錄組技術(shù)的不斷完善和新興的單細胞三代測序技術(shù)、空間轉(zhuǎn)錄組測序技術(shù)等新轉(zhuǎn)錄組技術(shù)的不斷涌現(xiàn)和發(fā)展為轉(zhuǎn)錄組數(shù)據(jù)分析技術(shù)的發(fā)展提出了更高的要求；另一方面，人們對科學問題理解的逐步加深又要求轉(zhuǎn)錄組分析技術(shù)不斷創(chuàng)新，能更充分更深入地挖掘基因表達數(shù)據(jù)中蘊含的現(xiàn)象和規(guī)律。可以預見，轉(zhuǎn)錄組分析技術(shù)在未來將發(fā)揮更大的作用。

3 表觀組測序數(shù)據(jù)分析

隨著測序技術(shù)的發(fā)展，表觀組學相關(guān)研究已經(jīng)成為生命科學領(lǐng)域炙手可熱的研究方向?？茖W家們發(fā)現(xiàn)，表觀組學對許多重要的細胞過程的調(diào)控起著關(guān)鍵作用，在精準醫(yī)學、生長發(fā)育、分子育種、公共安全等領(lǐng)域研究中占有不可忽視的地位。表觀組研究主要包括組蛋白修飾、DNA 甲基化、染色質(zhì)可及性，以及染色質(zhì)三維結(jié)構(gòu)等。表觀組測序?qū)嶒灱夹g(shù)的特殊性決定了用于數(shù)據(jù)分析的各種計算方法的特征和特異性。隨著表觀組研究技術(shù)發(fā)展的突飛猛進，必將開發(fā)出更多新的算法，產(chǎn)生更多解決問題的新工具用于數(shù)據(jù)整合和深度挖掘，也必將使得更多的生命科學秘密被發(fā)現(xiàn)，更多的改變生命過程的表觀遺傳標記被發(fā)現(xiàn)并應用于提高國計民生。

3.1 組蛋白修飾

目前，組蛋白修飾研究主要通過免疫共沉淀技術(shù)與深度測序技術(shù)相結(jié)合的辦法（ChIP-Seq）獲取實驗數(shù)據(jù)。ChIP-Seq 數(shù)據(jù)分析中的核心步驟是富集峰識別（peak calling）。首先是通過整合比對到特定基因組區(qū)域的讀段數(shù)來生成信號圖譜。然后，依靠滑動窗口方法將讀段計數(shù)的離散分布平滑為連續(xù)的信號輪廓分布，估計超出預定義的峰的讀段數(shù)量[81]，進一步考慮正鏈和負鏈中讀段計數(shù)的對應關(guān)系[82-83]，以提高峰的分離度。還有其他一些工具使用更復雜的方法在序列窗口中整合信號，例如使用局部泊松模型來識別基因組位置的局部偏差[84]、依賴核密度估計[85]、使用貝葉斯分層t-混合模型用于平滑基因組信號圖譜中的讀段計數(shù)[86]。在與樣本的比較分析中，背景分布的選擇也是識別富集峰中必不可少的步驟，多使用非結(jié)合蛋白質(zhì)的對照抗體實驗的ChIPSeq 數(shù)據(jù)。大多數(shù)峰識別算法通過估計其相應的p值和q值來對更重要的峰進行排序和選擇[87]。通常，不同的峰識別工具產(chǎn)生的結(jié)果會有所不同，因此，使用者需根據(jù)研究的靶向蛋白和實驗的情況，比如陰性對照樣本的選擇、是否有生物學重復等來選擇峰識別的工具。

3.2 DNA 甲基化

通常情況下，DNA 甲基化是指胞嘧啶甲基化。檢測DNA 甲基化狀態(tài)的實驗主要包括重亞硫酸鹽芯片、MeDIP-Seq、重亞硫酸鹽測序BS-Seq 等等。Illumina 公司的Infinium Methylation Assay 能夠?qū)θ蚪M范圍內(nèi)的單堿基分辨度的甲基化水平進行定量測量。在處理Infinium 芯片測定的數(shù)據(jù)時，要進行背景校正，以消除非特異性信號和重復樣品之間的差異[88]，完成歸一化[88]和進行批次效應校正[89]等。MeDIP-Seq是利用5-甲基胞嘧啶特異性抗體識別甲基化DNA的免疫沉淀技術(shù)結(jié)合高通量測序的方法。MeDIP-Seq 數(shù)據(jù)可以采用與ChIP-Seq 數(shù)據(jù)相同的生物信息學方法識別富集區(qū)域，也可以使用專門為甲基化富集數(shù)據(jù)量身定制的方法[89]。在重亞硫酸鹽測序數(shù)據(jù)中，甲基化的胞嘧啶不與重亞硫酸鹽發(fā)生化學反應，而未甲基化的胞嘧啶則在重亞硫酸鹽處理和測序后變成胸腺嘧啶。對于BS-Seq 數(shù)據(jù)，將已經(jīng)發(fā)生堿基變化的讀段比對到基因組，是甲基化數(shù)據(jù)分析的第一個難點。目前，已開發(fā)出多種比對策略來解決這個問題。一類是利用已有通用比對軟件[90-91]，采用三堿基（即T，G和A 或者A，T和C）方式進行序列比對[92-94]。由于三堿基比對降低了序列的特異性，在比對到大的參考基因組序列時，一方面增加了比對時間，另一方面減少了唯一比對的讀段數(shù)量，造成測序讀段利用率變低。一旦重亞硫酸鹽測序讀段與參考基因組唯一比對，就可以估計特定基因組區(qū)域的甲基化水平，定量Cs和Ts的頻率。另一種比對策略是在比對過程中，不斷調(diào)整標記比對評分的矩陣，以適應堿基不匹配的情況[94]。這類比對工具往往高估了高甲基化的區(qū)域，但比對效率較高。在為BS-Seq 數(shù)據(jù)選擇比對工具時，要考慮測序數(shù)據(jù)是來自于定向還是非定向建庫測序。在BSSeq 建庫時往往會加入非甲基化的Lambda DNA，用于比對后評估重亞硫酸鹽反應的轉(zhuǎn)化效率。近年來，單細胞實驗技術(shù)和測序技術(shù)迅猛發(fā)展，單細胞甲基化研究多采用重亞硫酸鹽建庫方式，測序數(shù)據(jù)的比對率一直都比較低。2019年一項研究發(fā)現(xiàn)[95]，在單細胞BS-Seq 文庫制備時，通過重亞硫酸鹽轉(zhuǎn)化后基因組近端序列與微同源區(qū)（MR）重組產(chǎn)生嵌合分子，而且MR 內(nèi)的DNA 甲基化高度可變，必須刪除這些區(qū)域以準確估算DNA 甲基化水平。而常規(guī)比對工具采用end-to-end 比對，難以達到較高的比對率。于是科學家們開發(fā)出支持單細胞甲基化測序數(shù)據(jù)局部比對的比對工具，在執(zhí)行質(zhì)控和重亞硫酸鹽測序數(shù)據(jù)局部比對的同時，進行嵌合分子測定和MR 去除，大大提高了單細胞重亞硫酸鹽測序數(shù)據(jù)的比對率和功能元件的回收率。

對于DNA 甲基化測序數(shù)據(jù)的分析，如何在全基因組水平上準確識別DNA 甲基化差異區(qū)域（Differentially Methylated Region,DMR）仍然是一個挑戰(zhàn)。DMR的識別通常采用以下幾種方法，第一種使用二項分布、負二項分布或具有過度分散參數(shù)的離散分布對甲基化/未甲基化讀段的數(shù)量進行建模[96-98]。一種是應用平滑算子（smoothing operator），考慮相鄰CpG 位點之間甲基化模式的相關(guān)性[98]。Metilene[99]采用分割算法來檢測單個或者一組重復實驗之間的DMR，不需要對數(shù)據(jù)生成機制進行任何模型假設。雖然識別DMR的方法很多，但是往往存在參數(shù)的依賴性，選擇不同的參數(shù)，會得到不同的DMR 識別結(jié)果。為了解決參數(shù)依賴和缺乏通用性的問題，科學家們采用完全貝葉斯方法開發(fā)工具ABBA[14]自動平滑甲基化圖譜，并可靠地識別DMR。ABBA 基于潛在的高斯模型（Latent Gaussian Model,LGM）和集成嵌套拉普拉斯近似（Integrated Nested Laplace Approximation,INLA），對WGBS 實驗的隨機采樣過程（甲基化和非甲基化的讀段數(shù)量分布作為非高斯響應變量）進行建模，采用概率分布指定所有未知量。ABBA 通過對平滑的未觀測到的甲基化圖譜的后驗概率的評估來計算每個CpG 位點的后驗概率，并對兩組數(shù)據(jù)之間的全基因組后驗甲基化概率進行差異比較，進而識別DMRs。ABBA 考慮了WGBS 數(shù)據(jù)的一些內(nèi)在特征，比如通過具有特定組內(nèi)差異的隨機效應進行建模以消除組內(nèi)實驗重復之間的DNA 甲基化差異；通過潛在的高斯場方程反映模型的鄰域結(jié)構(gòu)，并自動適應數(shù)據(jù)的變化，進而明確DNA 甲基化模式的相關(guān)性。ABBA為潛在的高斯場方程的參數(shù)分配先驗分布，從而充分考慮了這些量的不確定性。所有這些功能使ABBA的模型能夠根據(jù)實際情況進行調(diào)整，而無需任何用戶定義的參數(shù)。盡管此方法實現(xiàn)了DMR 識別的無參數(shù)，但是算法繁瑣復雜，并且對于大樣本間的差異甲基化識別仍有很大的困難和挑戰(zhàn)。

3.3 染色質(zhì)可及性

目前，研究染色質(zhì)可及性的實驗手段主要是通過酶解或者物理化學手段對開放區(qū)域的DNA 進行片段化處理，進一步對DNA 片段進行基因組定位，明確富集區(qū)域，并對富集區(qū)域進行功能分析。目前研究染色質(zhì)可及性的方法主要包括DNase-Seq、FAIRE-Seq、MNase-Seq和ATAC-Seq。MNase-Seq是研究核小體定位的方法，通過對核小體保護區(qū)域的DNA 序列進行定位，進而反映出失去核小體保護的DNA的區(qū)域，其他三種方法則是通過對染色質(zhì)開放區(qū)域的直接定位，反映染色質(zhì)開放性。分析這幾種實驗類型測序數(shù)據(jù)的思路與ChIP-Seq 測序數(shù)據(jù)的分析思路基本上是一致的，即識別富集區(qū)域，并對富集區(qū)域進行功能分析。但是染色質(zhì)開放性峰信號通常是寬序列讀取峰，需要區(qū)分真實峰和假象峰，也要對峰進行分類[100-102]。科學家們不斷開發(fā)和優(yōu)化算法，使染色質(zhì)可及性測序數(shù)據(jù)的分析越來越準確。近年來，單細胞染色質(zhì)可及性研究越來越多[103]，單細胞染色質(zhì)可及性對于在單細胞水平了解細胞核內(nèi)染色質(zhì)動態(tài)變化、揭示染色質(zhì)可及性在細胞間的異質(zhì)性、鑒定細胞特異性的染色質(zhì)活性區(qū)域等具有重要意義。

3.4 染色質(zhì)三維結(jié)構(gòu)

三維基因組結(jié)構(gòu)研究利用的染色質(zhì)構(gòu)象技術(shù)在經(jīng)歷了3C、4C和5C 技術(shù)后，又開發(fā)出了Hi-C 技術(shù)。通過Hi-C 數(shù)據(jù)分析揭示染色質(zhì)的遠程相互作用，推導出基因組的三維空間結(jié)構(gòu)和可能的基因之間的調(diào)控關(guān)系。Hi-C 數(shù)據(jù)的上游分析包括：（1）比對基因組，Hi-C 數(shù)據(jù)需要用特定的策略來比對基因組[104-105]；（2）基因組的一定區(qū)間范圍（bin）的選擇，bin的大小直接決定了最終分析結(jié)果的分辨度。由此，也產(chǎn)生了一些確定最佳bin 大小的方法，使選擇的bin 盡可能的小[106]；（3）歸一化，如迭代校正和特征向量分解歸一化，順序分量歸一化等[104-105,107-108]，并生成數(shù)據(jù)矩陣。Hi-C 數(shù)據(jù)的下游分析是從Hi-C 數(shù)據(jù)矩陣中提取有生物學意義的結(jié)果，包括：（1）隔室（A/B Compartments）識別，隔室識別多通過主成分分析獲得，由第一主成分值的正負來表示A Compartment或者B Compartment，也有依靠一種更快且內(nèi)存效率更高的方法來定義隔室得分，以反映出任何給定的bin 進入“A”隔室的可能性[105,109]；（2）拓撲相關(guān)結(jié)構(gòu)域（TAD）識別，即識別沿著接觸矩陣的對角線顯示為高度自關(guān)聯(lián)的連續(xù)區(qū)域[106,110-112]；（3）相互作用點識別，也就是識別距離較遠的染色質(zhì)區(qū)域之間相互聯(lián)系的特定位點。相互作用的計算識別需要定義背景模型，以便于識別相互作用頻率高于預期的相互作用聯(lián)系[113-114]。另外，數(shù)據(jù)的可視化在Hi-C 數(shù)據(jù)的分析中顯得尤為重要，可視化可以直觀地幫助研究者描述TAD 模式，也可以揭示分析產(chǎn)生的基因組結(jié)構(gòu)特征性信息。因此，科學家們不斷地開發(fā)一些可視化工具用于Hi-C 數(shù)據(jù)的分析，比如

Hi-Browse[115]、HicPlotter[116]、3D-GNOME[117]、HiC-3Dviewer[118]和3D genome browser[119]等以熱圖的方式可視化染色質(zhì)相互作用矩陣數(shù)據(jù)，并大多采用web 形式實現(xiàn)人機交互。近兩年開發(fā)的，如Delta[120]通過使用4個可視化視圖模塊，即Virtual 4C 繪圖、線性基因組視圖、基因組環(huán)形視圖和物理視圖，實現(xiàn)對不同類型的特征數(shù)據(jù)進行全面的可視化展示；GITAR[121]除了可實現(xiàn)Hi-C 數(shù)據(jù)的全面分析和可視化，還提供大量人類和小鼠的數(shù)據(jù)集，便于用戶進行數(shù)據(jù)集之間的比較；HiCcompare[122]能夠可視化分析兩個Hi-C 數(shù)據(jù)集之間的差異；GenomeFlow[123]在實時可視化Hi-C 數(shù)據(jù)3D 基因組模型的同時，還允許用戶在3D 基因組模型上附上基因注釋、基因表達數(shù)據(jù)和基因組甲基化數(shù)據(jù)；SilkDB 3.0[124]專門用于家蠶數(shù)據(jù)的交互式可視化分析?？傮w而言，Hi-C 技術(shù)的迅速發(fā)展，帶來了染色質(zhì)三維結(jié)構(gòu)數(shù)據(jù)集的快速增加，這也加快了數(shù)據(jù)分析方法的迅速興起，然而，雖然可以使用大量的方法學解決方案來識別TAD，但尚未完全了解TAD的結(jié)構(gòu)和功能，尤其是TAD的內(nèi)部結(jié)構(gòu)仍然難以捉摸[125]。因此，全面了解染色質(zhì)三維結(jié)構(gòu)的功能和作用，面臨著生物學和技術(shù)上的雙重挑戰(zhàn)。

4 挑戰(zhàn)和展望

從本世紀初高通量組學測序技術(shù)開始得到廣泛應用至今，隨著測序花費的下降和測序通量的不斷提高，各國政府和企業(yè)資助的大型基因組、轉(zhuǎn)錄組、表觀組計劃相繼開展[16-17,78,126]，對基因組學數(shù)據(jù)分析技術(shù)處理海量數(shù)據(jù)的速度提出了極高的要求。為了應對基因組學數(shù)據(jù)的飛速增長，一方面，算法技術(shù)不斷改進，如數(shù)據(jù)分析中耗時步驟是將讀段比對到基因組的算法，從線性比對[127]，到對基因組建立散列索引[128]，到建立后綴數(shù)組/Burrows-Wheeler 變換索引，到建立更加復雜的FM 索引[129]，到更高效的層次化圖結(jié)構(gòu)的FM 索引（Hierarchical Graph FM index，HGFM）[130]，一系列改進使得序列比對速度有了上千倍的提高；另一方面，專門針對基因組數(shù)據(jù)分析設計的并行計算和高性能計算方案也開始涌現(xiàn)，比如使用多核共享內(nèi)存的超級計算機、特殊的硬件設備（FPGA、GPU、TPU 等）、多節(jié)點的高性能計算機、云計算、容器部署等[45]。已有文章針對上述硬件開發(fā)了應用于多種基因組學數(shù)據(jù)分析的并行算法，比如基于Apache Spark的GATK RNA-Seq流程、基于迭代隨機森林的基因表達網(wǎng)絡構(gòu)建、基于MapReduce的應用于年齡預測的CpG 島篩選、基于Apache Spark的對BWA-MEM的加速等，都達到了很好的效果。盡管有上述這些進展，現(xiàn)有的算法和硬件技術(shù)發(fā)展對處理數(shù)據(jù)效率的提升仍然遠遠趕不上組學數(shù)據(jù)增長的速度，很多分析流程的復雜性決定了計算過程必須是并行和串行混合的，并且有很難突破的限速步驟，并行計算和分布式計算本身也存在理論上的極限[131]。因此，如何開發(fā)新的方法，提高組學大數(shù)據(jù)的處理效率，仍是目前亟待解決的問題。另外，實踐“功能即服務”的理念——把復雜的數(shù)據(jù)分析流程分割成微服務，部署在商業(yè)云或?qū)I(yè)云上，可以增加每部分服務的可伸縮性，對于世界各地需要進行組學數(shù)據(jù)分析的中小實驗室，可以減去其購買高性能服務器等基礎設施建設和搭建軟件環(huán)境的壓力，提高其數(shù)據(jù)分析的效率，也是未來的發(fā)展趨勢之一。

海量組學數(shù)據(jù)分析除了數(shù)據(jù)量大帶來的計算問題，數(shù)據(jù)來源廣產(chǎn)生的異構(gòu)性、不完整、尺度大、質(zhì)量參差不齊，也給組學大數(shù)據(jù)分析帶來了極大的挑戰(zhàn)。目前已有很多去除組學大數(shù)據(jù)批次效應的方法和模型[132-133]，但是統(tǒng)計建模過程從特定數(shù)據(jù)類型，轉(zhuǎn)變?yōu)槎嘟M學整合交叉研究的自動化建模是未來發(fā)展的方向。針對異構(gòu)性和數(shù)據(jù)缺失問題，雖然已提出了多核學習（multiple kernel learning）算法、DARTS 計算框架、有參模型、無參模型等大數(shù)據(jù)機器學習、人工智能、統(tǒng)計建模技術(shù)，為多源異構(gòu)組學數(shù)據(jù)構(gòu)建預測模型，但是算法的性能、精度和可靠性仍有很大的提升空間[15,134-135]。測序技術(shù)的發(fā)展也將推動組學數(shù)據(jù)分析方法的發(fā)展，但是如何整合分析多組學數(shù)據(jù)全面解讀生物系統(tǒng)也成了一項挑戰(zhàn)[13]。多組學數(shù)據(jù)整合分析將帶來維度災難、擴展性、歸一化等問題，有待優(yōu)化和開發(fā)一些深度學習方法促進多維組學大數(shù)據(jù)整合分析[136-137]。此外，臨床上將多組學數(shù)據(jù)與環(huán)境暴露數(shù)據(jù)、健康醫(yī)療檔案信息、臨床影像信息聯(lián)合分析將能更好的服務于精準醫(yī)學的發(fā)展，雖然已經(jīng)有了一些研究范式[17,135]，但是制定數(shù)據(jù)標準化體系、建設生物醫(yī)學大數(shù)據(jù)共享平臺、開發(fā)多維生物醫(yī)學數(shù)據(jù)分析的新算法、建立多維知識圖譜等都是亟待解決的問題[138]。我國于2019年成立國家基因組科學數(shù)據(jù)中心來存儲、質(zhì)控、管理、發(fā)布、共享這些多維、異構(gòu)基因組學數(shù)據(jù)[139]，期待該中心未來也能開發(fā)出能解決數(shù)據(jù)多維、異構(gòu)問題的算法和工具。

利益沖突聲明

所有作者聲明不存在利益沖突關(guān)系。