結構關系挖掘及其在生物信息中的應用

陳章昭，陳未如，張 雪，高勝召，韓 靜

（1.沈陽化工大學 計算機科學與技術學院，遼寧 沈陽 110142；2.遼寧省化工過程工業(yè)智能化技術重點實驗室，遼寧 沈陽 110142）

0 引 言

結構關系挖掘是數(shù)據(jù)挖掘領域一個全新的分支，是基于序列模式挖掘提出的一種尋找序列模式之間內在結構關系的挖掘方法。該方法將序列模式之間的關系進一步細分，整合成一種由并發(fā)、互斥、重復及串行關系組成的復合關系[1-3]。

隨著人類基因組計劃的啟動和高通量測序技術的快速發(fā)展，生物醫(yī)學數(shù)據(jù)呈現(xiàn)指數(shù)增長趨勢，面對海量的生物數(shù)據(jù)，生物信息學的重要性日益凸顯[4]。生物信息學的研究內容主要包括發(fā)展新的數(shù)理信息技術以用于開發(fā)生物數(shù)據(jù)的算法和運用數(shù)據(jù)挖掘等計算機技術以用于分析解釋生物基因信息。如今，應用和開發(fā)數(shù)據(jù)挖掘技術來探索生物系統(tǒng)規(guī)律是生物信息學領域最受關注的方向之一，其中包括基因序列分析、蛋白質功能預測、生物進化分析、表達圖譜分析等[5-6]。在生物信息領域，結構關系挖掘方法也有重要應用。王翠青等人提出使用支持向量作為提取蛋白質序列中新模式的算法ConSP[7]并進行了并發(fā)挖掘。Jing Lu 等人使用真實蛋白質數(shù)據(jù)集的實驗突顯了ConSP 方法在蛋白質這種數(shù)據(jù)中的適用性[8]?，F(xiàn)有的結構關系挖掘方法關注的是序列模式間的結構關系，忽略了那些并不是很頻繁但卻可能有意義的序列，而且在挖掘過程中，主要分析的是序列模式之間而非序列內部的關系，這在面向生物信息挖掘時可能會導致分析得到的結果過于冗余，實際意義不大。因此，本文對現(xiàn)有的結構關系挖掘知識體系做了進一步改進，在原有結構關系挖掘的基礎上，改進了并發(fā)度、互斥度以及并發(fā)關系和互斥關系的概念，基于此提出了面向生物基因信息的結構關系挖掘算法框架。改進后的結構關系挖掘方法將序列之間的結構關系進一步細化到項集之間，并且關注了那些并不頻繁但可能存在意義的序列。這樣的改動使得在面向生物基因信息挖掘時能夠得到更加科學客觀的結果，從而確保在研究生物結構、分析生物進化變異等問題時能夠快速準確地挖掘到有效的知識。

1 相關問題描述

1.1 有關序列模式和結構關系模式的知識

I={i1,i2, ...,im}是項目的集合，項集是I的非空子集，記為（x1,x2, ...,xk），其中xj∈I，（1 ≤j≤k≤m）。序列S是項集的有序集合，記為{s1,s2, ...,sn}，其中每個元素si是一個項集。在事務數(shù)據(jù)庫中，包含S的序列數(shù)與事務數(shù)據(jù)庫中的序列總數(shù)之比稱為序列S的支持度，記為sup(S)。用戶指定的最小支持度記為minsup。當序列S的支持度大于等于用戶指定的最小支持度，即sup(S)≥minsup 時，則稱序列S為頻繁序列或序列模式[9]。事務數(shù)據(jù)庫中所有的序列模式構成該數(shù)據(jù)庫的序列模式集，記為SPDB。對于序列S={I1I2...In}和序列S'={I'1I'2...I'n}，m＜n，如果存在m個正整數(shù)1 ≤j1≤j2≤...≤jm≤n，使得Ij1?I'j1，Ij2?I'j2, ...,Ijm?I'jm，則稱序列S包含于序列S'，記為S?S'，也稱S為S'的子序列或S'為S的超序列[10]。

結構關系模式挖掘是一種基于序列模式挖掘所提出的挖掘任務，旨在尋找隱藏在序列模式間的結構關系，如并發(fā)關系模式、互斥關系模式以及重復關系模式等[1-3]?，F(xiàn)有結構關系模式挖掘的研究給出了并發(fā)度、互斥度、并發(fā)序列模式以及互斥序列模式等定義，并在此基礎上提出了幾種結構關系模式挖掘算法。

1.2 結構關系的相關概念

基因項：基因項ij=＜loc, base＞由兩部分組成, 其中base為生物堿基{A, G, C, T}或蛋白質{G, A, V, L, I, F, W, Y, D, N,E, K, Q, M, S, T, C, P, H, R}，loc 為該堿基或蛋白質在所對應序列的絕對位置。

基因序列：由若干個基因項構成的集合稱為基因序列，記為{i1,i2, ...,in}，其中每個元素ij為一個基因項。

如序列S：{＜1, A＞, ＜2, T＞, ＜3, T＞, ..., ＜210, T＞, ＜211, T＞,＜212, G＞, ..., ＜29 561, T＞, ＜29 562, G＞}是由29 562 個基因項構成的新冠病毒基因序列。

并發(fā)度：對于序列A={α1,α2, ...,αn}，序列數(shù)據(jù)庫SDB中包含A的序列個數(shù)與包含A中任意項集的序列個數(shù)之比，稱為序列A的并發(fā)度，記作con(α1,α2, ...,αn)，或con(A)。

并發(fā)關系：對于序列A={α1,α2, ...,αn}，給定客戶指定的最小并發(fā)度mincon，當con(A)≥mincon 時，稱A存在并發(fā)關系，表示為[A]=[α1+α2+...+αn]。α1,α2, ...,αn構成一組并發(fā)集。特別的，若A為基因序列且并發(fā)集中包含n個基因項，則稱該并發(fā)集為n-基因并發(fā)集。

表1 為包含了4 條新冠序列的基因序列數(shù)據(jù)庫GSDB。

表1 基因序列數(shù)據(jù)庫GSDB

若給定最小并發(fā)度mincon = 70%, 根據(jù)并發(fā)度的定義可以得出序列S={＜1, A＞, ＜210, T＞, ＜211, T＞}的并發(fā)度con(＜1,A＞＜210, T＞, ＜211, T＞)=3/4 ≥mincon, 則稱序列S存在并發(fā)關系。表示為[S]=[＜1, A＞+＜210, T＞+＜211, T＞]。＜1, A＞, ＜210,T＞, ＜211, T＞構成一組3-基因并發(fā)集。

并發(fā)關系具有反單調性：對于給定的序列數(shù)據(jù)庫GSDB，如果序列A={α1,α2, ...,αn}存在并發(fā)關系[α1+α2+...+αn]，則A的任意一個子序列也存在并發(fā)關系。

證明：假設序列A={α1,α2, ...,αn} 且存在并發(fā)關系[α1+α2+...+αn]，即con(α1,α2, ...,αn)≥mincon，A'為序列A的一個n-1 子序列。在序列庫SDB 中，包含A的序列肯定也包含A'，即con(A')的分子要大于等于con(A)的分子；由于序列A'較A相比少了一個元素，因此，con(A')的分母要小于等于con(A)的分母。綜上可得，con(A')≥con(A)≥mincon，即序列A的任意一個n-1 序列也存在并發(fā)關系。以此類推，序列A的任意一個子序列都存在并發(fā)關系。

完全并發(fā)集：對于并發(fā)關系[C1]=[α1+α2+...+αm]和[C2]=[β1+β2+...+βn]，m＜n。若對?i（1 ≤I≤m）都存在αi?βj（1 ≤j≤n），則稱并發(fā)關系[C2]包含并發(fā)關系[C1]。若基因序列S存在并發(fā)關系且不被任意一個并發(fā)關系所包含，則稱并發(fā)關系[S]為完全并發(fā)關系，該并發(fā)關系的所有基因項構成一組完全并發(fā)集。

互斥度：對于序列A={α1,α2, ...,αn}，序列數(shù)據(jù)庫SDB中包含且僅包含A中一個項集的序列個數(shù)與包含A中任意項集的序列個數(shù)之比稱為序列A的互斥度，記作xcl(α1,α2, ...,αn)或xcl(A)。

互斥關系：對于序列A={α1,α2, ...,αn}，給定客戶指定的最小并發(fā)度minxcl，當xcl(A)≥minxcl 時，稱A存在互斥關系，表示為[A]=[α1⊕α2⊕...⊕αn]。α1,α2, ...,αn構成一組互斥集。特別的，若A為基因序列且互斥集中包含n個基因項，則稱該互斥集為n-基因互斥集。

對于給出的GSDB，若給定最小互斥度minxcl =60%，根據(jù)互斥度的定義可以得出基因序列S={＜29 655,C＞, ＜29 656, A＞} 的互斥度為：xcl(＜29 655, C＞, ＜29 656,A＞)=2/3 ≥minxcl，稱序列S存在互斥關系，表示為[S]=[＜29 655, C＞⊕＜29 656, A＞]。

根據(jù)互斥度與并發(fā)度關系，以及并發(fā)關系的反單調性質可知，任何一個存在互斥關系的序列A的超序列（包含該A的序列）很容易滿足互斥關系，滿足這一條件的互斥關系稱為平凡互斥關系，這樣的互斥關系不是我們關心的，只有那些任意子序列間都存在互斥關系的序列才有意義。

非凡互斥關系：序列A={α1,α2, ...,αn}存在非凡互斥關系，當且僅當A及其所有子序列都滿足互斥關系。顯然，非凡互斥關系滿足反單調性。

關聯(lián)度：同時包含序列A和B的序列占包含序列A的序列的比例，稱為序列A關聯(lián)B的關聯(lián)度，記作association(A,B)。

關聯(lián)關系：對于序列A與B，當A在某一序列中出現(xiàn)時B也有很大概率出現(xiàn)，即A與B的關聯(lián)度ass(A,B)≥minass（minass 為客戶指定的最小關聯(lián)度），則稱序列A與B存在關聯(lián)關系，表示為[A→B]。

對于給出的GSDB，若給定最小關聯(lián)度minass =90%，根據(jù)關聯(lián)度的定義可以得出基因序列A={＜1, A＞, ＜2, T＞, ＜2,C＞}與B={＜210, T＞, ＜212, G＞}的關聯(lián)度為：association(A,B)=1.0 ≥minass， 稱存在關聯(lián)關系[{＜1, A＞, ＜2, T＞, ＜2,C＞}→{＜210, T＞, ＜212, G＞}]。

2 結構關系挖掘具體過程描述

獲取數(shù)據(jù)集：數(shù)據(jù)集一般是通過訪問資料庫、網頁抓取和問卷調查手動收集等方式獲得。特別的，對于生物基因數(shù)據(jù)而言，可以從生物基因數(shù)據(jù)庫下載進行研究，如NCBI（https：//www.ncbi.nlm.nih.gov/）為美國國家生物技術信息中心，該數(shù)據(jù)庫包含人類基因組、病毒、微生物和新冠病毒等生物基因信息；GISAID（https：//db.cngb.org/gisaid/）是全球最大的流感及新型冠狀病毒數(shù)據(jù)平臺，該數(shù)據(jù)庫不僅具有最完整的新冠病毒基因組序列數(shù)據(jù)以及相關臨床和流行病學數(shù)據(jù)，更匯聚了全球諸多科研團隊對COVID-19 的研究成果。

預處理：數(shù)據(jù)預處理是為了提高數(shù)據(jù)的質量，保證數(shù)據(jù)的準確性、完整性和一致性。對于生物基因序列而言，從基因數(shù)據(jù)庫下載的基因序列可能存在基因缺失或未知堿基數(shù)過多的情況，所以需要將下載的序列進行預處理，去除重復和低質量的序列。其次，雖然同類生物基因序列相似度很高，但序列長度會存在略微偏差，在預處理階段還需要進行序列對齊操作。

獲取變異基因序列組：生物進化的實質是遺傳物質的變異。在面向生物信息的研究過程中，由于生物序列和事務序列的特征存在很大差異，如生物序列是由有限個體（堿基或蛋白質）組成的超長序列，且同類生物基因序列的相似性很高[11]，因此可以對基因序列的變異點進行針對性研究，這不僅可以很大程度提高挖掘效率，更使得分析生物變異進化過程更加科學、客觀。已有研究表明，病毒基因組間的共突變是研究病毒進化的重要標志。例如，Deng Lizong 等人利用氨基酸序列的共突變網絡來預測埃博拉病毒的致命性[12]。Olabode E. Omotoso 等人分析發(fā)現(xiàn)，新冠病毒序列S蛋白D614G 與其他復發(fā)性蛋白共突變對病毒ACE-2 宿主進入產生了影響[13]。Qin Luyao 等人根據(jù)SARS-CoV-2 基因組發(fā)現(xiàn)了一些共突變模塊來推測病毒的進化傳播過程[14]。經過預處理后的基因序列長度一致，我們首先選取一條序列作為參考序列，然后將基因組中的序列和選取的參考序列進行序列比對，去除序列中具有一致核苷酸的保守位點，剩余部分則構成了變異基因序列。所有變異基因序列構成變異序列組vGSDB。

挖掘序列間的結構關系：通過結構關系挖掘算法挖掘序列間的結構關系。本文給出了面向生物基因信息的結構關系挖掘算法。

結果可視化表達：由挖掘序列間的結構關系步驟可以得到SDB 結構關系，根據(jù)這些結構關系通過相關可視化方法可進一步分析序列庫中的信息。如本文根據(jù)挖掘得到的新冠序列結構關系生成了GSDB 系統(tǒng)發(fā)育框架，從而更好捕獲生物基因的進化變異情況。

3 結構關系挖掘算法

由于實驗對象選取的是新冠病毒Sars-Cov-2 序列，因此本文在Apriori、Prefixspan 等序列模式挖掘算法的基礎上，結合結構關系定義給出了適用于生物信息領域的挖掘算法。

3.1 基于Apriori 的并發(fā)關系挖掘算法conApriori

輸入：基因序列數(shù)據(jù)庫GSDB，最小并發(fā)度mincon

輸出：所有的并發(fā)集allConcurrentItemSets

算法：

（1）獲取GSDB 中所有的基因項T，生成初始候選并發(fā)集Lk = T?T（k 為2）；

（2）令Ck = null； allConcurrentItemSets = null；

（3）do

for each s of Lk

if（con（s）≥mincon）

將s 存入k-并發(fā)集Ck 中；

將Ck 存入allConcurrentItemSets；

Lk+1=Ck?Ck；

for each c of Lk+1

if（c 存在k-子序列不被Ck 所包含）

將c 從Lk+1 中刪除；

while（Lk+1 is not null） ；

3.2 基于PrefixSpan 的并發(fā)關系挖掘算法conPrefix

輸入：基因序列數(shù)據(jù)庫GSDB，最小并發(fā)度mincon

輸出：所有的并發(fā)集allConcurrentItemSets

算法：

（1）獲取GSDB 中所有的基因項T；

（2）令pre = null，preDB = null，prefDBItem = null；

（3）for each s of T

（4）pre = s；

（5）獲取前綴pre 對應的投影數(shù)據(jù)庫preDB；

if（preDB is not null）

獲取preDB 中的所有基因項preDBItem；

for each c of preDBItem

pre = pre + c；

if （con（pre）≥mincon）

將pre 存入allConcurrentItemSets；

執(zhí)行步驟（5）；

3.3 基于Apriori 的互斥關系挖掘算法excApriori

輸入：基因序列數(shù)據(jù)庫GSDB，最小互斥度minxcl

輸出：所有非凡互斥集allExclusionaryItemSets

算法：

（1）獲取GSDB 中所有的基因項T，生成初始候選互斥集Lk = T?T（k 為2）；

（2）令Ek = null； allExclusionaryItemSets= null；

（3）do

for each s of Lk

if（xcl（s）≥minxcl）

將s 存入k-非凡互斥集Ek 中；

將Ek 存入allExclusionaryItemSets；

Lk+1=Ek?Ek；

for each c of Lk+1

if（c 存在k-子序列不被Ek 所包含）

將c 從Lk+1 中刪除；

4 實驗

4.1 數(shù)據(jù)準備和處理

GISAID 數(shù)據(jù)庫按照Pangolin 分類方法將SARS-CoV-2基因組序列劃分成了若干個分支。本實驗首先將GenBank 中的SARS-CoV-2 參考基因組作為參考序列，再從B 分支，B.1分支以及B.1.126 分支分別選取15 條SARS-CoV-2 基因組序列作為數(shù)據(jù)集。數(shù)據(jù)集經數(shù)據(jù)清洗和序列對齊等預處理操作生成新的基因序列集GSDB。將基因序列集GSDB 中的每條序列和選取的參考序列進行序列比對，去除序列中具有一致核苷酸的保守位點，剩余部分組成變異基因序列。處理得到45 條變異基因序列，構成變異基因序列集vGSDB。

4.2 實驗過程及分析

通過conApriori，conPrefix 以及excApriori 算法對變異基因序列集vGSDB 進行結構關系挖掘。

表2 列舉了變異基因序列集vGSDB 在最小并發(fā)度mincon=0.9 下挖掘得到的部分并發(fā)集。由表2 可知，多處核苷酸或氨基酸位點存在高并發(fā)變異的情況，如ORF1ab 蛋白的4517H 和S 蛋白的614G 之間滿足高并發(fā)關系。已有研究表明，SARS-CoV-2 序列產生并發(fā)突變可能會對病毒的傳播和感染能力產生一定影響，如S 蛋白的L452R 和T478K 并發(fā)突變會使刺突蛋白以更高的親和力附著在ACE2 受體上，從而影響病毒的傳 播感染性[15]。

表2 vGSDB 在mincon=0.9 下挖掘得到的部分并發(fā)集

圖1 表示數(shù)據(jù)集GSDB 的系統(tǒng)發(fā)育情況。在mincon=0.6時，vGSDB 通過算法挖掘得到了15 個完全并發(fā)集。構建一個n×m的矩陣，其中n代表選取的45 條樣本基因序列，m代表挖掘得到的15 個完全并發(fā)集。若基因序列中存在某完全并發(fā)集，則將該完全并發(fā)集對應位置的值記為1，否則記為0。然后，先將具有相同完全并發(fā)集的基因序列放在同一組，再根據(jù)組中包含完全并發(fā)集的個數(shù)將基因序列分配到不同水平。若一個組包含i個完全并發(fā)集，則將該組分配到第i級。同時，如果第i級中的完全并發(fā)集都被第i+1 級的完全并發(fā)集包含，則將第i級的組認定為第i+1 級組的父級。根據(jù)該劃分規(guī)則，樣本基因序列集最終被劃分為15 個分支。現(xiàn)有一些對SARS-CoV-2 進行分類的方法，如Nextstrain 和Pangoli[16]分別根據(jù)分離株和突變的數(shù)量以及時空分布的變化來構建系統(tǒng)發(fā)育樹，從而實現(xiàn)對SARS-CoV-2 的劃分。本文提出的根據(jù)完全并發(fā)集的劃分方法可以看作是一種更加精細的劃分規(guī)則，屬于同一分支下的序列親緣性較高，該方法能夠準確有效地識別SARSCoV-2 群體之間的層次關系，這為病毒追根溯源提供了相關依據(jù)。

圖1 數(shù)據(jù)集GSDB 的系統(tǒng)發(fā)育情況

圖2 表示conPrefix 和conApriori 兩種挖掘算法在不同相關度下的運行時間效率情況。由圖2 可知，隨著最小并發(fā)度的不斷增大，兩種挖掘算法所消耗的時間也不斷減少。此外，在同一最小并發(fā)度下，conPrefix 算法消耗的時間相對較短，效率優(yōu)于conApriori 算法。

圖2 conPrefix 和conApriori 算法在不同并發(fā)度下的運行時間變化曲線

圖3 表示變異基因序列集vGSDB 通過算法在不同minxcl 下挖掘得到的非凡互斥集個數(shù)變化情況。從圖3 可得知，在最小互斥度為0.9 時，依然存在著大量的非凡互斥集，如{＜1822, I＞, ＜5951, Y＞, ＜69, C＞}、{＜2611, I＞, ＜5951, Y＞,＜69, C＞, ＜167, F＞}均滿足非凡互斥關系，這表明病毒發(fā)生了進一步突變，部分序列出現(xiàn)了較新突變位點。

圖3 vGSDB 在不同最小互斥度下的非凡互斥集變化曲線

圖4 表示變異基因序列集vGSDB 通過算法在不同minass 下挖掘得到的關聯(lián)關系變化情況。從圖4 可得知，在不同的最小關聯(lián)度下，挖掘得到的關聯(lián)關系個數(shù)相差不大，說明變異序列間的關聯(lián)性很強，如關聯(lián)關系{＜5951, Y＞ ＜73,F＞} →{＜4715, H＞ ＜614, G＞}，當ORF1ab 蛋白的5951 號位氨基酸突變?yōu)閅，E 蛋白的73 號位氨基酸突變?yōu)镕 時，ORF1ab 蛋白的4715 號位氨基酸突變?yōu)镠，S 蛋白的614 號位氨基酸突變?yōu)镚 的可能性很大。

圖4 vGSDB 在不同最小關聯(lián)度下的關聯(lián)關系變化曲線

5 結 語

在生物信息領域，結構關系是一種普遍存在的關系，如新冠病毒序列S 蛋白D614G 與其他復發(fā)性蛋白并發(fā)變異會對病毒ACE-2 宿主進入產生影響；癌癥患者的基因中，具有互斥關系的基因集合在患者群體中會表現(xiàn)出有且只有一個基因變異的現(xiàn)象；妊娠乳腺癌患者的基因存在關聯(lián)關系，檢測到其中一種修飾，則其它三種也極有可能存在。本文在現(xiàn)有結構關系挖掘方法的基礎上做了進一步改進，給出了面向生物信息的結構關系挖掘算法，并將其應用到SARS-CoV-2研究中，通過實驗挖掘得到了隱藏在序列集中的并發(fā)、互斥和關聯(lián)關系，由并發(fā)關系生成的系統(tǒng)進化樹可用于新冠病毒序列間的進化傳播關系研究，也進一步驗證了挖掘方法的正確有效性。此外，結構關系挖掘還可以應用于諸多場景，合理的應用結構關系挖掘方法對生物基因信息研究具有一定意義。

注：本文通訊作者為陳章昭。