基于矩陣的并發(fā)關(guān)系挖掘算法在新冠病毒變異分析中的應(yīng)用

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

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

0 引 言

新冠病毒變異株中最受關(guān)注的是Delta、Omicron、Alpha、Gamma、Beta、Lambda、Mu 變異株。以上變異株或是增加了毒性，或是增強了傳播性，某些變異株還會降低疫苗效用，對公共衛(wèi)生構(gòu)成重大威脅[1]。因此，對于新冠病毒突變的研究至關(guān)重要。

考慮到DNA 序列的特殊性，已經(jīng)有許多專門針對DNA序列挖掘的算法，如ToMMSA、ATRHunter 等[2]。機器學(xué)習(xí)也被應(yīng)用到DNA 序列挖掘[3]。至今已經(jīng)有許多研究人員在針對SARS-CoV-2 的研究中獲得了成果。Qin 等人[4]提出了一種識別共突變模塊的算法，并找到了42 個共突變模塊，基于這些共突變模塊，將SARS-CoV-2 種群分為43 組，并根據(jù)不一致的共突變模塊的數(shù)量確定了各組之間的系統(tǒng)發(fā)育關(guān)系。Liu 等人[5]提出一種計算并發(fā)度的公式。Mishra 等人[6]對SARS-CoV-2 序列進行聚類分析，將突變分為3 類：共同發(fā)生的突變、前導(dǎo)和尾隨突變、相互排斥的突變。

本文以SARS-CoV-2 為例進行生物變異結(jié)構(gòu)關(guān)系挖掘，將共突變模塊的挖掘與結(jié)構(gòu)關(guān)系模式相結(jié)合，提出并發(fā)變異關(guān)系挖掘，并提出更加適合生物變異結(jié)構(gòu)關(guān)系性質(zhì)的公式，嘗試從結(jié)構(gòu)化的角度來挖掘出新冠病毒基因序列的變異信息，并從算法的時間復(fù)雜度與空間復(fù)雜度進行討論。

1 結(jié)構(gòu)關(guān)系模式挖掘

結(jié)構(gòu)關(guān)系模式挖掘是在序列模式挖掘基礎(chǔ)上提出的一種新的數(shù)據(jù)挖掘任務(wù)。結(jié)構(gòu)關(guān)系是一種包括并發(fā)關(guān)系、互斥關(guān)系、關(guān)聯(lián)關(guān)系及這些關(guān)系的復(fù)合關(guān)系的形式，用于表示序列模式之間存在的關(guān)系。并發(fā)關(guān)系[7-9]是指兩個或兩個以上序列大概率同時出現(xiàn)在同一個大的序列中，也就是說它們都是某幾個客戶序列的子序列?；コ怅P(guān)系[10]是指兩個或兩個以上序列大概率不出現(xiàn)在同一個大的序列中。關(guān)聯(lián)關(guān)系是指當(dāng)一個序列出現(xiàn)時另一個序列很大概率也在客戶序列中出現(xiàn)[11]。

2 新冠病毒基因序列以及挖掘算法

2.1 新冠病毒基因序列

新冠病毒基因組具有典型的冠狀病毒結(jié)構(gòu)，大約含有3 萬個核苷酸，核苷酸G 與C 占比約為40%，大約編碼9 860 個氨基酸[12]。截至2022 年1 月NCBI（美國國家生物信息中心）上現(xiàn)存兩百多萬條新冠病毒基因序列。

2.2 變異序列組的構(gòu)建

為了消除臟數(shù)據(jù)，便于構(gòu)建變異堿基矩陣，方便以后的生物變異結(jié)構(gòu)關(guān)系挖掘，我們基于選定的基因數(shù)據(jù)庫VGDB，對原有的病毒基因序列組進行處理與簡化，構(gòu)建出變異序列組。在進行DNA 序列挖掘時大部分算法是對整體DNA 序列進行挖掘，挖掘其序列模式[13-14]。由于DNA 序列很長，所以挖掘時間也較長。在此我們提取變異堿基構(gòu)建變異序列組，只挖掘變異部分，大大縮減了挖掘時間。

定義1 變異序列：在對比病毒基因序列與參考基因序列后，將與參考基因序列不同的堿基（變異點）提取出來，表示為形式，B表示變異后的堿基種類或堿基缺失，S表示變異堿基在病毒基因序列的位點。由變異點組成的序列叫做變異序列，表示為{α1,α2, ...,αn}，其中αi是變異點。

例1：基因序列為AGGAAC…；參考基因序列為AGCATT…；變異序列為{、、...}。

定義2 變異序列組：基因數(shù)據(jù)庫中所有基因序列的變異序列構(gòu)成變異序列組，記為VGDB。

定義3 變異矩陣：儲存病毒基因序列堿基變異情況的N×M的01 矩陣，N行代表所有的病毒基因序列，M列代表此病毒基因序列組所有的變異堿基點，0 代表此基因序列不存在這個變異堿基點，1 代表此基因序列存在這個變異堿基點。變異矩陣的構(gòu)建大大減少了在挖掘中遍歷事務(wù)數(shù)據(jù)庫VGDB 的次數(shù)[15]。

例2：假設(shè)的變異堿基矩陣見表1 所列，有S1、S2、S3三條基因序列，有、、三處變異堿基點。

表1 變異堿基矩陣

由此變異矩陣可得基因序列S1 的變異堿基為、、，基因序列S2 的變異堿基為、，基因序列S3 的變異堿基為、。

變異序列組構(gòu)建步驟：首先從NCBI 數(shù)據(jù)庫獲得新冠病毒基因序列，為了盡可能消除抽樣偏差，我們根據(jù)抽樣日期和國家對序列進行分類。如果基因序列的含糊堿基（未確定是“A”“C”“T”還是“G”，顯示為“N”）過長，則刪除這些基因序列，這樣就獲得了新冠病毒基因序列組。然后對新冠病毒基因序列組進行序列對齊，去除對齊后基因組開頭和結(jié)尾的非編碼區(qū)[4]。設(shè)置一條參考基因序列（參考基因序列應(yīng)選取比新冠病毒基因序列組日期靠前、株系相近或相同、病毒采取地相同的一條新冠病毒基因序列），將新冠病毒基因序列組與參考基因序列對比，獲得變異堿基位點。為了便于表示變異堿基與其堿基位點，提出一種變異堿基點表示形式：，B表示變異后的堿基種類，S表示變異堿基在病毒基因序列的位點。將變異堿基全部轉(zhuǎn)換為此種新形式，由此變異序列組構(gòu)建完成。

2.3 生物變異結(jié)構(gòu)關(guān)系挖掘

生物變異結(jié)構(gòu)關(guān)系由三部分組成：并發(fā)變異關(guān)系、互斥變異關(guān)系、關(guān)聯(lián)變異關(guān)系。在此我們結(jié)合結(jié)構(gòu)關(guān)系模式提出以下概念。

定義4 變異率：存在變異點α的變異序列與變異序列組VGDB 序列數(shù)之比叫做變異率，記作p(α)。

定義5 共變序列：變異點α1、α2、…、αn在變異序列s中同時存在，即A={α1,α2, ...,αn}是s的子序列，稱變異序列A是在s中的共變序列，或稱s支持序列A。

定義6 共變率：包含序列A={α1,α2, ...,αn}的變異序列數(shù)與變異序列組VGDB 序列數(shù)之比，記作p(A)。

定義7 并發(fā)度：包含序列A={α1,α2, ...,αn}的變異序列數(shù)與存在A中任意變異點的變異序列數(shù)之比，稱為序列A的并發(fā)變異度，簡稱并發(fā)度，記作conv(α1,α2, ...,αn)，或conv(A)。

式中，分子是包含序列A={α1,α2, ...,αn}的變異序列數(shù)，分母是至少包含A中任一變異點的變異序列數(shù)。分子分母各除以變異序列組VGDB 序列數(shù)，則分子就是p(A)。

式（3）中分母的計算比較復(fù)雜，考慮簡化計算，將包含k個變異點的序列稱為k-序列。序列A={α1,α2, ...,αn}有n個1-序列、n(n-1)/2 個2-序列、…、n個(n-1)-序列和1 個n-序列，這些統(tǒng)稱為A的k子序列。

設(shè)序列A={α1,α2, ...,αn}的各k子序列的并發(fā)率之和分別表示為p(Ak)，k=1, 2, ...,n，則有：

根據(jù)容斥原理的性質(zhì)，式（3）可以改寫為：

例3：假設(shè)變異序列組為：

則有：

根據(jù)式（4），可得：

定義8 并發(fā)關(guān)系：序列A={α1,α2, ...,αn}，對于客戶指定的最小并發(fā)度minconv，當(dāng)conv(A)≥minconv 時，稱A存在并發(fā)變異關(guān)系，簡稱并發(fā)關(guān)系，表示為[A]=[ɑ1+ɑ2+...+ɑn]。

在例3 中，設(shè)minconv=0.5，則有并發(fā)關(guān)系[+]、[+]和[+]存在。

并發(fā)關(guān)系有如下性質(zhì)：

（1）性質(zhì)1：并發(fā)關(guān)系具有反單調(diào)性。若序列A={α1,α2, ...,αn}存在并發(fā)關(guān)系，則其任意子序列也一定存在并發(fā)關(guān)系。

證 明：設(shè) 序 列A={ɑ1,ɑ2, ...,ɑn} 存 在 并 發(fā) 關(guān) 系[ɑ1+ɑ2+...+ɑn]，即conv(A)≥minconv，A′為A的一個n-1子序列（k≤n）?？芍赩GDB 中，包含A的所有序列一定包含A′，所以p(A′)≥p(A)。由于少了一個累計變異點，式（3）、（4）的分母值將變?。ㄖ辽俨蛔兇螅?。因此，conv(A′)≥conv(A)≥minconv。即存在并發(fā)關(guān)系的n序列的所有n-1 子序列存在并發(fā)關(guān)系。依此類推，存在并發(fā)關(guān)系的序列的所有子序列存在并發(fā)關(guān)系。反單調(diào)性質(zhì)成立。

為了利用并發(fā)關(guān)系的反單調(diào)性進行并發(fā)關(guān)系挖掘，定義候選并發(fā)關(guān)系。

定義9 候選并發(fā)關(guān)系：若序列的所有子序列都存在并發(fā)變異關(guān)系，則這個變異序列構(gòu)成候選并發(fā)變異關(guān)系，簡稱候選并發(fā)關(guān)系。

根據(jù)并發(fā)關(guān)系的反單調(diào)性質(zhì)，任意一個變異序列存在共變關(guān)系的前提是其所有子系列存在共變關(guān)系，即它首先應(yīng)該是一個候選并發(fā)關(guān)系。通過并發(fā)變異率矩陣很容易得到所有的二元并發(fā)變異關(guān)系，然后再以所有的二元并發(fā)變異關(guān)系為基礎(chǔ)組成三元候選并發(fā)關(guān)系集合，從中篩選出滿足條件的三元并發(fā)關(guān)系。以此類推，可以逐步得到所有并發(fā)關(guān)系。

定義10 單變序列：變異點α1、α2、...、αn在變異序列s中存在且只存在其中一個點，稱A={α1,α2, ...,αn}為在s中的單變序列，該單變序列包含在s中。

定義11 互斥度：包含單變序列A={α1,α2, ...,αn}的變異序列數(shù)與存在A中任意變異點的變異序列數(shù)之比，稱為單變序列A的互斥度，記作 xclv(α1,α2, ...,αn) = xclv(A)。

根據(jù)容斥原理的性質(zhì)，與式（4）類似，將式（5）改寫為：

定義12 互斥關(guān)系：序列A={α1,α2, ...,αn}，對于客戶指定的最小互斥度minxclv，當(dāng)xclv(A)≥minxclv 時，稱A存在互斥變異關(guān)系，簡稱互斥關(guān)系，表示為[A]=[ɑ1⊕ɑ2⊕...⊕ɑn]。

在例3 中，則有：

設(shè)minconv=0.8，則有互斥關(guān)系[⊕]、[⊕]和[⊕]存在。

定義13 關(guān)聯(lián)度：同時包含序列A和B的變異序列數(shù)與包含A的變異序列數(shù)之比，稱為序列A關(guān)聯(lián)B的關(guān)聯(lián)變異度，簡稱關(guān)聯(lián)度，記作 assv(A,B)。

定義14 關(guān)聯(lián)關(guān)系：對于序列A與B，當(dāng)A在某一變異序列中出現(xiàn)時B也有很大概率出現(xiàn)，即A與B的關(guān)聯(lián)度assv(A,B)≥min assv（客戶指定的最小關(guān)聯(lián)變異度），則稱序列A與B滿足關(guān)聯(lián)變異關(guān)系，簡稱關(guān)聯(lián)關(guān)系，表示為[A→B]。

在例3 中,則有：

若 設(shè)minassv=0.9， 則 關(guān) 聯(lián) 關(guān) 系[{, }→]存在。

3 基于矩陣的并發(fā)變異關(guān)系挖掘算法步驟

（1）首先通過變異序列組獲得此新冠病毒基因序列組所有的變異堿基點集合（allVariationBases）。

（2）然后構(gòu)建一個N×M的變異矩陣（variationMatrix），N行代表所有的新冠病毒基因序列VGDB，M列代表此新冠病毒基因序列組所有的變異堿基點集合。矩陣的元素為0 或1，0 代表此基因序列不存在這個變異堿基點，1 代表此基因序列存在這個變異堿基點。

（3）首先進行二元并發(fā)變異關(guān)系挖掘，通過并發(fā)度計算公式計算任意兩個變異堿基的并發(fā)度，將所得結(jié)果構(gòu)建M×M的并發(fā)度矩陣（convMatrix），M行和M列均代表此新冠病毒基因序列組所有的變異堿基點。矩陣的元素是公式計算結(jié)果。

（4） 由 客 戶 指 定 最 小 并 發(fā) 度minconv， 當(dāng)conv(A)≥minconv 時，求得所有的二元并發(fā)變異關(guān)系。

（5）通過步驟（4）得到的二元并發(fā)關(guān)系組成三元候選并發(fā)關(guān)系集合，從中篩選出滿足條件的三元并發(fā)關(guān)系。以此類推，逐步得到所有的并發(fā)關(guān)系。

4 實 驗

4.1 數(shù)據(jù)選取與處理

從NCBI 數(shù)據(jù)庫上下載日本地區(qū)120 條新冠病毒基因序列，根據(jù)日期與株系分為3 組，每組40 條，第一組（日期：2020年7月到2021年11月；株系：B.1.1.214），第二組（日期：2021 年1 月到2021 年5 月；株系：B.1.1.7），第三組（日期：2021 年6 月到2021 年9 月；株系：B.1.1.7）。以日本2020年5 月獲取的一條B.1.1 病毒序列作為參考序列，將3 組新冠病毒基因序列組與參考序列進行對齊，去除對齊后基因組開頭和結(jié)尾的非編碼區(qū)。將新冠病毒基因序列組與參考基因序列對比，獲得所有的變異堿基位點，將變異堿基點全部變?yōu)?B,S>形式，3 組變異序列組構(gòu)建完成。

4.2 挖掘并發(fā)變異關(guān)系

（1）獲取3 組的變異堿基點集合：第一組有93 個變異堿基點；第二組有95 個變異堿基點；第三組有156 個變異堿基點。

（2）分別獲取3 組的變異矩陣：第一組獲得40*93 的變異矩陣；第二組獲得40*95 的變異矩陣；第三組獲得40*156的變異矩陣。

（3）獲取3 組的并發(fā)度矩陣：第一組獲得93*93 的并發(fā)度矩陣；第二組獲得95*95 的并發(fā)度矩陣；第三組獲得156*156 的并發(fā)度矩陣。

（4）設(shè)最小并發(fā)度minconv=0.95，挖掘二元并發(fā)變異關(guān)系：第一組挖掘到10 個二元并發(fā)變異關(guān)系；第二組挖掘到276個二元并發(fā)變異關(guān)系；第三組挖掘到276個二元并發(fā)變異關(guān)系。

(5) 挖掘多元并發(fā)變異關(guān)系：第一組挖掘到一個5 元并發(fā)變異關(guān)系為[++++]；第二組挖掘到一組24元并發(fā)變異關(guān)系為[+++++++++++++++++++++++]；第三組挖掘到一組24 元并發(fā)變異關(guān)系為[+++++++++++++++++++++++]。

三組之間存在共同并發(fā)變異關(guān)系：[+]；第二組與第三組并發(fā)變異關(guān)系相同。

4.3 實驗結(jié)果分析

對本次挖掘結(jié)果進行分析：首先，從變異堿基點集合來看，第一組與第二組變異堿基點數(shù)量相近，第三組的變異堿基點數(shù)量遠多于前兩組，僅從變異堿基點數(shù)量來看，第一組和第二組與參考序列的同源性更高。然后，從并發(fā)變異關(guān)系來看，第一組挖掘到的并發(fā)變異關(guān)系數(shù)量遠遠小于后兩組，93 個變異堿基點只挖掘到10 個二元并發(fā)變異關(guān)系，說明在第一組中堿基較少出現(xiàn)并發(fā)變異，多為獨立變異，變異堿基的普遍度不高。第二組與第三組挖掘到的并發(fā)變異關(guān)系相同，說明第二組與第三組屬于同一株系，再由變異堿基點數(shù)量來看，三組數(shù)量多，說明三組較二組與參考序列的同源性更低，三組是在二組的基礎(chǔ)上變異的。第一組與二、三兩組有共同的并發(fā)變異關(guān)系，說明第一組與第二、三兩組屬于同一個大的株系。分析實驗結(jié)果可知，三組株系關(guān)系與實際情況相同，符合B.1.1、B.1.1.214、B.1.1.7 之間的親緣關(guān)系。B.1.1.7 曾是值得關(guān)注的株系，現(xiàn)為正在監(jiān)測的株系（新冠病毒關(guān)注級別由低至高：正在監(jiān)測的變種、值得留意的變種和值得關(guān)注的變種）；B.1.1.214 為普通株系(株系等級由疾病預(yù)防與控制中心官網(wǎng)得到)；B.1.1.7 的毒性或傳播性大于B.1.1.214。

為了避免本次實驗存在偶然性，所以我們又從印度選取120 條新冠病毒序列做一個對照實驗。根據(jù)日期與株系分為3 組，每組40 條，第一組（日期：2020 年7 月到2021 年12 月；株系：B.1.1.306），第二組（日期：2020 年11 月到2021 年1月；株系：B.1.1.7），第三組（日期：2021 年2 月到2021 年5 月；株系：B.1.1.7）。參考序列選取印度2020 年6 月獲取的一條B.1.1 病毒序列。第一組獲得112 個變異堿基點和6組二元并發(fā)關(guān)系；第二組獲得105 個變異堿基點和359 組二元并發(fā)關(guān)系；第三組獲得153 個變異堿基點和377 組二元并發(fā)關(guān)系。通過對挖掘出的變異堿基和并發(fā)關(guān)系進行分析，三組株系屬于同一個父株系，二、三組屬于同一株系，并且第三組是在第二組的基礎(chǔ)上變異而來，符合B.1.1、B.1.1.306、B.1.1.7 之間的關(guān)系。B.1.1.7 曾為值得關(guān)注的株系，現(xiàn)為正在監(jiān)測的株系；B.1.1.306 為普通株系；B.1.1.7 的毒性或傳播性大于B.1.1.306。

針對日本與印度病毒序列的兩組實驗的結(jié)果大致相同，根據(jù)挖掘到的變異堿基與并發(fā)關(guān)系可以推斷出變異序列之間的同源性，所得結(jié)果與實際大致相符。本文的實驗結(jié)果表明，并發(fā)變異關(guān)系可能是驅(qū)動不同株系致病性或傳播性的主要進化力量，并且可能是決定這些株系毒性或傳播性程度更高和更低的原因。

4.4 算法效率分析

時間復(fù)雜度：在求并發(fā)變異關(guān)系時，時間復(fù)雜度主要由構(gòu)建變異矩陣、構(gòu)建并發(fā)度矩陣、生成候選并發(fā)變異關(guān)系和對候選并發(fā)變異關(guān)系計算并發(fā)度四部分構(gòu)成。

構(gòu)建變異矩陣：假設(shè)變異序列組存在m條變異序列，這些變異序列分別與數(shù)量為n的變異堿基集合進行比對，時間復(fù)雜度為O(mn2)。

構(gòu)建并發(fā)度矩陣：計算conv 需要O(m)，需要計算n個變異堿基，由于并發(fā)度矩陣是中心對稱的，只需計算矩陣的上三角形部分，所以時間復(fù)雜度為O(mn)。

生成候選并發(fā)變異關(guān)系：時間復(fù)雜度為O(mn2)。

對候選并發(fā)變異關(guān)系計算并發(fā)度：時間復(fù)雜度為O(n2)，所以算法總的時間復(fù)雜度為O(mn2)。

基于Apriori 的并發(fā)變異關(guān)系挖掘算法在生成每個候選并發(fā)變異關(guān)系時，要通過遍歷事務(wù)數(shù)據(jù)庫來獲得并發(fā)度，時間復(fù)雜度為O(mn2)。因此，當(dāng)生成ɑ個候選并發(fā)變異關(guān)系時，基于Apriori 的并發(fā)變異關(guān)系挖掘算法的要比基于矩陣的并發(fā)變異關(guān)系挖掘算法多O(ɑmn2)。

比較基于矩陣的并發(fā)變異關(guān)系挖掘算法（M-alg）和基于Apriori 的并發(fā)變異關(guān)系挖掘算法（A-alg）對不同數(shù)量的變異序列組進行挖掘?qū)嶋H所需時間。結(jié)果如圖1 所示。

圖1 并發(fā)關(guān)系算法挖掘不同序列條數(shù)運行時間比較

由于挖掘?qū)ο笫亲儺悏A基而不是整條基因序列，并且在最開始將事務(wù)數(shù)據(jù)庫轉(zhuǎn)變?yōu)樽儺惥仃嚕诰驎r只需遍歷一遍事務(wù)數(shù)據(jù)庫即可，所以挖掘時間較基于Apriori 的并發(fā)變異關(guān)系算法大大減少。

空間復(fù)雜度由兩部分構(gòu)成：（1）變異矩陣：需要一個n*m的二維數(shù)組，空間復(fù)雜度為O(nm)；（2）并發(fā)度矩陣：需要一個n*n的二維數(shù)組，空間復(fù)雜度為O(n2)。因此，算法總的空間復(fù)雜度為O(n2)（目前實驗數(shù)據(jù)量：變異堿基數(shù)量n與新冠病毒序列數(shù)量m的范圍均為0 ～200）。就目前新冠病毒變異堿基數(shù)量與我們所需要挖掘的新冠病毒序列條數(shù)來看用內(nèi)存來存儲完全夠用，所以在追求較快運行速度的前提下優(yōu)先使用內(nèi)存來存儲，在以后挖掘時遇到內(nèi)存不夠時也會選用外存來存儲數(shù)據(jù)。

5 結(jié) 語

由于環(huán)境和宿主的選擇壓力，一部分基因一起發(fā)生突變，這樣就形成了共突變，就相當(dāng)于本文的并發(fā)變異關(guān)系。在癌癥細胞中存在互斥變異關(guān)系，在同一致病通路上的基因只需要變異一個即可以導(dǎo)致功能異常，因此多個基因的變異會使得功能冗余，不具有選擇優(yōu)勢[16]。在生物變異的過程中一個基因的變異會導(dǎo)致另一個基因的變異，這就是本文中的關(guān)聯(lián)變異關(guān)系。因此，對生物變異的并發(fā)變異關(guān)系、互斥變異關(guān)系、關(guān)聯(lián)變異關(guān)系的挖掘是非常有必要的。本文通過挖掘新冠病毒基因序列的并發(fā)變異關(guān)系，可以判斷新冠病毒序列之間的同源性，證明結(jié)構(gòu)關(guān)系模式挖掘是可以應(yīng)用到生物變異信息挖掘的；此外就目前數(shù)據(jù)規(guī)模而言，算法效率是可以接受的。