基于矩陣的并發關系挖掘算法在新冠病毒變異分析中的應用

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

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

0 引 言

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

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

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

1 結構關系模式挖掘

結構關系模式挖掘是在序列模式挖掘基礎上提出的一種新的數據挖掘任務。結構關系是一種包括并發關系、互斥關系、關聯關系及這些關系的復合關系的形式，用于表示序列模式之間存在的關系。并發關系[7-9]是指兩個或兩個以上序列大概率同時出現在同一個大的序列中，也就是說它們都是某幾個客戶序列的子序列。互斥關系[10]是指兩個或兩個以上序列大概率不出現在同一個大的序列中。關聯關系是指當一個序列出現時另一個序列很大概率也在客戶序列中出現[11]。

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

2.1 新冠病毒基因序列

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

2.2 變異序列組的構建

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

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

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

定義2 變異序列組：基因數據庫中所有基因序列的變異序列構成變異序列組，記為VGDB。

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

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

表1 變異堿基矩陣

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

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

2.3 生物變異結構關系挖掘

生物變異結構關系由三部分組成：并發變異關系、互斥變異關系、關聯變異關系。在此我們結合結構關系模式提出以下概念。

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

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

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

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

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

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

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

根據容斥原理的性質，式（3）可以改寫為：

例3：假設變異序列組為：

則有：

根據式（4），可得：

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

在例3 中，設minconv=0.5，則有并發關系[+]、[+]和[+]存在。

并發關系有如下性質：

（1）性質1：并發關系具有反單調性。若序列A={α1,α2, ...,αn}存在并發關系，則其任意子序列也一定存在并發關系。

證 明：設 序 列A={ɑ1,ɑ2, ...,ɑn} 存 在 并 發 關 系[ɑ1+ɑ2+...+ɑn]，即conv(A)≥minconv，A′為A的一個n-1子序列（k≤n）。可知在VGDB 中，包含A的所有序列一定包含A′，所以p(A′)≥p(A)。由于少了一個累計變異點，式（3）、（4）的分母值將變小（至少不變大）。因此，conv(A′)≥conv(A)≥minconv。即存在并發關系的n序列的所有n-1 子序列存在并發關系。依此類推，存在并發關系的序列的所有子序列存在并發關系。反單調性質成立。

為了利用并發關系的反單調性進行并發關系挖掘，定義候選并發關系。

定義9 候選并發關系：若序列的所有子序列都存在并發變異關系，則這個變異序列構成候選并發變異關系，簡稱候選并發關系。

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

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

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

根據容斥原理的性質，與式（4）類似，將式（5）改寫為：

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

在例3 中，則有：

設minconv=0.8，則有互斥關系[⊕]、[⊕]和[⊕]存在。

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

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

在例3 中,則有：

若 設minassv=0.9， 則 關 聯 關 系[{, }→]存在。

3 基于矩陣的并發變異關系挖掘算法步驟

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

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

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

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

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

4 實 驗

4.1 數據選取與處理

從NCBI 數據庫上下載日本地區120 條新冠病毒基因序列，根據日期與株系分為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 組新冠病毒基因序列組與參考序列進行對齊，去除對齊后基因組開頭和結尾的非編碼區。將新冠病毒基因序列組與參考基因序列對比，獲得所有的變異堿基位點，將變異堿基點全部變為形式，3 組變異序列組構建完成。

4.2 挖掘并發變異關系

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

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

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

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

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

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

4.3 實驗結果分析

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

為了避免本次實驗存在偶然性，所以我們又從印度選取120 條新冠病毒序列做一個對照實驗。根據日期與株系分為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組二元并發關系；第二組獲得105 個變異堿基點和359 組二元并發關系；第三組獲得153 個變異堿基點和377 組二元并發關系。通過對挖掘出的變異堿基和并發關系進行分析，三組株系屬于同一個父株系，二、三組屬于同一株系，并且第三組是在第二組的基礎上變異而來，符合B.1.1、B.1.1.306、B.1.1.7 之間的關系。B.1.1.7 曾為值得關注的株系，現為正在監測的株系；B.1.1.306 為普通株系；B.1.1.7 的毒性或傳播性大于B.1.1.306。

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

4.4 算法效率分析

時間復雜度：在求并發變異關系時，時間復雜度主要由構建變異矩陣、構建并發度矩陣、生成候選并發變異關系和對候選并發變異關系計算并發度四部分構成。

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

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

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

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

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

比較基于矩陣的并發變異關系挖掘算法（M-alg）和基于Apriori 的并發變異關系挖掘算法（A-alg）對不同數量的變異序列組進行挖掘實際所需時間。結果如圖1 所示。

圖1 并發關系算法挖掘不同序列條數運行時間比較

由于挖掘對象是變異堿基而不是整條基因序列，并且在最開始將事務數據庫轉變為變異矩陣，挖掘時只需遍歷一遍事務數據庫即可，所以挖掘時間較基于Apriori 的并發變異關系算法大大減少。

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

5 結 語

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