馬 揚 劉澤一 梁星星 程光權(quán) 陽方杰 成 清 劉 忠

(國防科技大學(xué)系統(tǒng)工程學(xué)院 長沙 410073)

基因組數(shù)據(jù)是目前最流行的生物信息數(shù)據(jù)類型之一，自人類基因組計劃開始，生物基因序列數(shù)據(jù)開始呈指數(shù)級增長.如何將有效的生物信息特征從海量高維的生物序列數(shù)據(jù)中提取出來，是目前生物信息領(lǐng)域最重要的問題之一.

近年來，基于基因序列的深度學(xué)習(xí)模型研究激增,主要原因之一是該類模型可以直接從原始數(shù)據(jù)進行表示學(xué)習(xí)或特征學(xué)習(xí).深度學(xué)習(xí)可以自動提取特征，省去了大量手工選擇特征的工作量[1].同時，深度學(xué)習(xí)模型也被證明可以有效識別局部基因序列在特定環(huán)境中的作用，并進一步可以用來推斷基因的調(diào)控規(guī)則和表達模式.自編碼器[2]、多層感知機(multilayer perception, MLP)[3-4]、玻爾茲曼機(restricted Boltzmann machines, RBMs)[5]、深度信念網(wǎng)絡(luò)(deep belief networks, DBNs)[6]、卷積神經(jīng)網(wǎng)絡(luò)(convolutional neural network, CNN)、循環(huán)神經(jīng)網(wǎng)絡(luò)(recurrent neural network, RNN)[7]、長短記憶網(wǎng)絡(luò)(long short-term memory, LSTM)[8]等深度學(xué)習(xí)方法與基因組研究的結(jié)合促進了對基因組學(xué)的深刻理解，并為精密醫(yī)學(xué)、藥物設(shè)計[9]等多個領(lǐng)域提供幫助.

Urda等人[10]使用深度學(xué)習(xí)方法在分析RNA序列基因表達數(shù)據(jù)方面優(yōu)于Lasso(least absolute shrinkage and selection operator).Chen等人[11]提出了一種3層前饋神經(jīng)網(wǎng)絡(luò)對基因表達進行預(yù)測，其性能優(yōu)于線性回歸.Xie等人[12]的研究表明，他們基于MLP和SDAs的深度模型在SNP基因型預(yù)測基因表達量化方面優(yōu)于Lasso和隨機森林.曹明宇等人[13]使用了一種新標(biāo)注模式,將藥物實體及關(guān)系的聯(lián)合抽取轉(zhuǎn)化為端對端的序列標(biāo)注任務(wù).Karlic等人[14-17]在實驗中提出了組蛋白修飾與基因調(diào)控之間的相關(guān)性，并且在之前的一些機器學(xué)習(xí)工作中已經(jīng)進行了研究.Singh等人[18]提出了將MLP疊加在CNN之上的統(tǒng)一判別框架DeepChrome，在預(yù)測基因表達水平高低的二元分類任務(wù)中，平均AUC(area under curve)為0.8. DeepTarget和deepMirGene分別使用RNN和LSTM模型使用表達數(shù)據(jù)進行mRNA目標(biāo)預(yù)測，所提方法的一個主要優(yōu)勢是不需要手工制作的特征集[19-20].與此同時，生成模型也被使用來捕捉基因序列中高階、隱藏的相關(guān)性，Way和Greene在癌癥基因組RNA序列數(shù)據(jù)上訓(xùn)練了一個變分自編碼模型來捕獲生物學(xué)相關(guān)特征[21-22].

然而，在實際應(yīng)用中，由于基因組序列之間存在著復(fù)雜的相互依賴和長期的相互作用，直接從基因組序列中學(xué)習(xí)特征非常耗時[23]，且成功的預(yù)測通常在很大程度上依賴于對生物學(xué)知識的正確利用[24].如何在更普遍的無監(jiān)督環(huán)境下構(gòu)建一個基因序列表示模型，使得模型學(xué)習(xí)到的基因表示可以更有效地表示基因特征是目前生物信息領(lǐng)域一個主要的研究問題.且目前基因序列的研究大都集中在基因表達、調(diào)控等領(lǐng)域上[25],對基因數(shù)據(jù)進行低維表示也會結(jié)合基因交互表達網(wǎng)絡(luò)的信息[26]，直接對病毒基因序列進行表示學(xué)習(xí)的研究并不多.

因此，本文針對病毒傳播網(wǎng)絡(luò)構(gòu)建了一種基于傳播網(wǎng)絡(luò)的基因序列表示模型(gene representation based on transmission network, GRTN)，將病毒的基因序列信息與傳播網(wǎng)絡(luò)的結(jié)構(gòu)信息相結(jié)合，進而得到新的基因序列表示.模型首先把輸入序列編碼至特定長度的隱藏層以得到降維嵌入，根據(jù)網(wǎng)絡(luò)傳播的特點使用正向LSTM和反向LSTM分別對節(jié)點的圖上下文節(jié)點信息進行聚合，然后使用注意力機制對其鄰居節(jié)點的序列信息進行聚合，進而得到目標(biāo)節(jié)點基因序列的新的表示.依據(jù)病毒傳播網(wǎng)絡(luò)中鄰居節(jié)點的基因序列相似性高于非鄰居節(jié)點這一特點，對基因序列表示模型進行訓(xùn)練和優(yōu)化.

該模型學(xué)習(xí)到的基因序列表示可以作為基因特征或先驗信息用于與基因序列相關(guān)的任務(wù)，為了驗證基因序列表示模型的有效性，本文設(shè)計一套在病毒傳播網(wǎng)絡(luò)中發(fā)現(xiàn)未采樣感染者的任務(wù)流程.并且通過在模擬的病毒傳播網(wǎng)絡(luò)、真實的澳大利亞新冠病毒傳播網(wǎng)絡(luò)和真實的艾滋病毒傳播網(wǎng)絡(luò)上進行了發(fā)現(xiàn)未采樣感染者的實驗，驗證了本文提出的基因序列表示模型的有效性.并通過與各類基因表示方法的對比實驗表明模型的高效性.

1 問題描述

病毒傳播網(wǎng)絡(luò)是一種復(fù)雜網(wǎng)絡(luò)，該網(wǎng)絡(luò)可以表示為G=(V,E,A).其中，V表示網(wǎng)絡(luò)中節(jié)點集合，每個節(jié)點表示一個病毒感染者；E表示連邊集合，每條邊(vi,vj)∈E表示節(jié)點vi將疾病傳染給節(jié)點vj，由于疾病傳播是有方向的，故該網(wǎng)絡(luò)為有向圖；A表示節(jié)點屬性集合，在G中，ai∈A表示節(jié)點i的病毒基因序列.

給定一個網(wǎng)絡(luò)G，基于網(wǎng)絡(luò)的基因序列表示的核心任務(wù)是將網(wǎng)絡(luò)的特征結(jié)合到基因序列中，從而為網(wǎng)絡(luò)中每一個節(jié)點的基因序列學(xué)習(xí)一個統(tǒng)一的低維度(β)的表示向量，從而將網(wǎng)絡(luò)中每個節(jié)點的屬性維度由|A|映射到β(β?|A|)，使其能夠捕捉網(wǎng)絡(luò)中的結(jié)構(gòu)和聚合鄰居節(jié)點的特征[27-28].

2 模型描述

本節(jié)介紹基于傳播網(wǎng)絡(luò)的基因序列表示模型，如圖1所示，該模型包括病毒網(wǎng)絡(luò)中節(jié)點病毒序列信息編碼、基于網(wǎng)絡(luò)傳播的結(jié)構(gòu)提取、基于圖上下文信息的特征表示和基于注意力的節(jié)點信息聚合.

Fig. 1 Gene sequence representation model based on transmission network圖1 基于傳播網(wǎng)絡(luò)的基因序列表示模型

作為病毒傳播網(wǎng)絡(luò)的輸入，節(jié)點的基因序列首先經(jīng)過獨熱編碼和維度規(guī)約轉(zhuǎn)化為適合模型輸入的數(shù)據(jù)；然后基于網(wǎng)絡(luò)傳播結(jié)構(gòu)對每個節(jié)點抽取子圖并獲得子圖中的前向節(jié)點信息、后向節(jié)點信息和自身表示信息；基于注意力機制，模型將前向信息、后向信息和節(jié)點本身的表示信息進行聚合，得到節(jié)點的聚合表示；為了能對網(wǎng)絡(luò)節(jié)點的病毒序列信息進行更好的表達，通過設(shè)置基于圖上下文信息的損失函數(shù)，對整個模型經(jīng)過優(yōu)化與訓(xùn)練，最終得到節(jié)點的表示向量.

2.1 病毒基因序列編碼

不同病毒的基因序列長度迥異，但目前已知的大多數(shù)病毒序列中，存在著大量的無效基因序列，而且在對基因序列進行測序時，往往也會存在部分序列信息缺失.同時，病毒在傳播中會不斷發(fā)生變異，不同網(wǎng)絡(luò)節(jié)點的病毒基因序列都存在一定差異.但同種病毒傳播網(wǎng)絡(luò)中，所有節(jié)點的基因序列大部分是相同的.針對病毒傳播網(wǎng)絡(luò)中基因序列的這種特征，本文采用獨熱編碼，該編碼對于不同節(jié)點的相同基因片段可以保證其一致性，同時可以有效地表達病毒基因序列的不同.

病毒在特定的傳播網(wǎng)絡(luò)中，存在基因序列相同的片段，本文通過對這些基因片段的刪減，對輸入數(shù)據(jù)維度進行了維度歸約：若該網(wǎng)絡(luò)中所有節(jié)點的某一位基因值相同，則將該位基因從輸入數(shù)據(jù)中舍棄.通過本節(jié)操作，網(wǎng)絡(luò)中任一節(jié)點v的病毒基因序列g(shù)v在規(guī)約后被編碼為相同長度的0/1編碼序列xv，作為網(wǎng)絡(luò)中節(jié)點的屬性向量.

2.2 網(wǎng)絡(luò)傳播結(jié)構(gòu)抽取

網(wǎng)絡(luò)傳播結(jié)構(gòu)抽取旨在抽取病毒傳播網(wǎng)絡(luò)中與每個節(jié)點關(guān)聯(lián)度高的鄰居節(jié)點集合，從而在保證計算效率的同時，使模型從網(wǎng)絡(luò)中獲取節(jié)點的拓?fù)湫畔?在該抽取的子圖中，目標(biāo)節(jié)點v以外的節(jié)點可以劃分為2類：前向信息節(jié)點集合Γf(v)和后向信息節(jié)點集合Γb(v):

Γf(v)=F(v)∪F(i),i∈F(v),

(1)

Γb(v)=S(v)∪S(i),i∈S(v),

(2)

其中,F(v)表示節(jié)點v的父節(jié)點，S(v)表示節(jié)點v的子節(jié)點.前向節(jié)點是位于病毒傳播鏈中目標(biāo)節(jié)點的祖先節(jié)點，即目標(biāo)節(jié)點的父節(jié)點及其祖父節(jié)點；后向節(jié)點是位于病毒傳播鏈中目標(biāo)節(jié)點的子孫節(jié)點，即目標(biāo)節(jié)點的子節(jié)點和孫子節(jié)點.由于病毒基因序列傳播變異較快，間隔多代后基因序列之間存在較大差異.此外，研究發(fā)現(xiàn)[29-31]，網(wǎng)絡(luò)中距離超過3的節(jié)點對之間的互相作用非常小，在本文模型中，選取與目標(biāo)節(jié)點距離小于3的節(jié)點，構(gòu)成目標(biāo)節(jié)點的鄰居集合.如圖1所示，節(jié)點子圖抽取得到的網(wǎng)絡(luò)中，中心節(jié)點到所有節(jié)點的距離均小于3，由于每個節(jié)點只有一個感染源，所以前向信息節(jié)點集合的規(guī)模較小(0～2個)，后向信息節(jié)點集合的規(guī)模相對更多.通過本節(jié)操作，模型中每個節(jié)點v的結(jié)構(gòu)信息只和抽取子圖中節(jié)點Γf(v)和Γb(v)相關(guān).

2.3 基于圖上下文的信息表示

基于圖上下文的信息表示旨在將每個子圖中的信息聚合在一起，子圖中的節(jié)點稱為目標(biāo)節(jié)點的上下文鄰居：節(jié)點的上文鄰居指沿著網(wǎng)絡(luò)中連邊經(jīng)過一定長度隨機游走能到達目標(biāo)節(jié)點的節(jié)點集合，節(jié)點的下文鄰居指從目標(biāo)節(jié)點出發(fā)，經(jīng)過指定長度的隨機游走能到達的節(jié)點集合.由于節(jié)點之間存在傳播關(guān)系，本文采用LSTM將傳播序列信息聚合在一起，對每個節(jié)點v，在計算前向表示向量時，直接采用LSTM將前向信息聚合在一起，并用LSTM最后的隱藏層輸出為前向表示向量；在計算后向表示時，在每個分支上，采用和前向表示類似的LSTM模型，用LSTM最后的隱藏層作為該分支的表示向量，并用所有分支向量的均值作為后向表示向量.對于自身向量，通過多層全連接層將輸入數(shù)據(jù)轉(zhuǎn)換為需要的維度.通過充分利用子圖中所有節(jié)點的信息，可以獲得前向表示Ef、后向表示Eb和自身表示Em三個表示向量:

Ef(v)=LSTMf(BN(FCf(xi))),i∈Γf(v),

(3)

(4)

Em(v)=BN(FC(xv)),

(5)

其中，BN表示批標(biāo)準(zhǔn)化操作，F(xiàn)C表示全連接層，3個表示向量的維度相同，記為do.通過該操作，模型獲得了每個目標(biāo)節(jié)點在網(wǎng)絡(luò)子圖中抽取的結(jié)構(gòu)信息的表示.

2.4 基于注意力的節(jié)點信息聚合

注意力機制已經(jīng)成為許多序列學(xué)習(xí)任務(wù)中的一項標(biāo)準(zhǔn)配置[32]，可以為不同節(jié)點分配不同權(quán)重，訓(xùn)練時依賴于成對的相鄰節(jié)點，而不依賴具體的網(wǎng)絡(luò)結(jié)構(gòu).本文采用注意力機制學(xué)習(xí)節(jié)點前向、自身和后向向量Ek，k∈{f,m,b}的聚合權(quán)重，進而獲得每個節(jié)點的信息聚合E.

E=AGGREGATE(Ef,Em,Eb),

(6)

其中

AGGREGATE(Ef,Em,Eb)=
αfEf+αmEm+αbEb,

(7)

αk,k∈{f,m,b}為注意力權(quán)重系數(shù)，本文通過一個共享的注意力向量a∈2d×1來反映向量Ek對節(jié)點自身向量Em的重要性:

emk=aT(Em⊕Ek),

(8)

其中⊕為拼接操作，a為可學(xué)習(xí)的注意力參數(shù).相應(yīng)的注意力系數(shù)αk可以表示為

(9)

其中σ為激活函數(shù)LeakyRelu.

為了減少數(shù)據(jù)的方差，提升表示結(jié)果的穩(wěn)定性，將該模式推廣到多頭注意力形式：注意力機制被獨立地重復(fù)nh次(每個頭的參數(shù)初始化不同)，并將得到的nh個表示向量拼接作為最終向量:

E=E1⊕E2⊕…⊕Enh,

(10)

2.5 目標(biāo)優(yōu)化與模型訓(xùn)練

為了客觀衡量病毒傳播網(wǎng)絡(luò)的基因序列表示的效果，本文采用圖上下文損失作為損失函數(shù).圖上下文損失可以看作交叉熵的一種變形：對每個節(jié)點，網(wǎng)絡(luò)中的所有其他節(jié)點，根據(jù)距離該節(jié)點的距離被劃分為正鄰居節(jié)點和負(fù)鄰居節(jié)點，因此其核心是使目標(biāo)節(jié)點的向量表示與采樣正鄰居的相似性更高，并與采樣負(fù)鄰居的差異性更大，基于此的優(yōu)化目標(biāo)被定義為

(11)

其中Cv表示從節(jié)點v開始的隨機游走能到達的節(jié)點集合，p(vc|v;Θ)被定義為softmax形式的概率:

(12)

式(12)通過表示向量內(nèi)積大小衡量節(jié)點對的相似性，通過負(fù)采樣[33]技術(shù)，可以將式(12)中softmax函數(shù)進行對數(shù)近似為

(13)

其中,S表示負(fù)采樣數(shù)目，P(vc)表示負(fù)采樣的概率分布，在本模型中，S=1從而將目標(biāo)函數(shù)簡約為交叉熵的形式[34]

lnθ(EvEvc)+lnθ(-EvEvs),

(14)

從而將損失函數(shù)轉(zhuǎn)化定義為

(15)

算法1.基因序列表示算法.

輸入：病毒傳播網(wǎng)絡(luò)G=(V,E,A)、節(jié)點基因序列g(shù);

輸出：網(wǎng)絡(luò)中每個節(jié)點的基因表示E.

① 通過病毒基因序列編碼將g轉(zhuǎn)化為x；

② for eachv∈V：

③ 分別按照式(1)～(2)抽取節(jié)點v的前向信息節(jié)點集合Γf(v)和后向信息節(jié)點集合Γb(v)；

④ 分別按照式(3)～(5)計算前向表示Ef、后向表示Eb和自身表示Em三個表示向量；

⑤ 按照式(8)計算每個節(jié)點的注意力系數(shù)αk,k∈{f,m,b}，并獲得每個節(jié)點聚合表示E；

⑦ end for

⑧ 根據(jù)式(14)進行模型訓(xùn)練和優(yōu)化；

⑨ 循環(huán)步驟②～⑧直至到達結(jié)束條件.

3 實驗設(shè)計

本節(jié)設(shè)計一種病毒傳播網(wǎng)絡(luò)中發(fā)現(xiàn)未采樣感染者的實驗流程，通過在仿真病毒傳播網(wǎng)絡(luò)、真實的新冠病毒傳播網(wǎng)絡(luò)和艾滋病毒傳播網(wǎng)絡(luò)上發(fā)現(xiàn)未采樣感染者來驗證模型的有效性，并對實驗結(jié)果進行分析.

3.1 未采樣感染者發(fā)現(xiàn)流程

實際情況下，受限于醫(yī)療水平、檢測條件、疾病癥狀的不同，檢測到的病毒傳播網(wǎng)絡(luò)往往會存在缺失節(jié)點(未采樣感染者)，其網(wǎng)絡(luò)規(guī)模小于實際的傳播網(wǎng)絡(luò).如圖2所示，虛線節(jié)點為未采樣感染者，由于這些病患節(jié)點在實際網(wǎng)絡(luò)構(gòu)建時未被發(fā)現(xiàn)，因此網(wǎng)絡(luò)的實際連邊關(guān)系存在異常，原本不應(yīng)有傳播關(guān)系的父節(jié)點和子節(jié)點被認(rèn)為存在直接傳播關(guān)系(存在連邊)，如右圖虛線所示.

Fig. 2 The un-sampled infections detection圖2 未采樣感染者示意圖

基于本文基因序列表示模型，我們設(shè)計2個未采樣感染者發(fā)現(xiàn)流程：

① 通過基于病毒傳播網(wǎng)絡(luò)的基因表示模型，獲取網(wǎng)絡(luò)中每個節(jié)點的基因表示E；

② 將發(fā)現(xiàn)傳播網(wǎng)絡(luò)中未采樣節(jié)點任務(wù)轉(zhuǎn)化為連邊得分排序任務(wù)，即通過計算兩節(jié)點的表示的相似度，評價這條連邊的穩(wěn)固性.對于節(jié)點對i，j，其連邊相似度為

simij=EiEj,

(16)

在所有存在連邊的排序中，得分最低的k條連邊所對應(yīng)的節(jié)點對被認(rèn)為是存在未采樣感染者的節(jié)點組合，進而在實際的醫(yī)學(xué)調(diào)查診斷相應(yīng)環(huán)節(jié)進行排查.

3.2 數(shù)據(jù)描述

本文采用仿真和真實的病毒傳播網(wǎng)絡(luò)數(shù)據(jù)集進行實驗.

仿真數(shù)據(jù)集是根據(jù)病毒傳播和變異的規(guī)律，生成傳播網(wǎng)絡(luò)和網(wǎng)絡(luò)中每個節(jié)點的基因序列，其中，每次仿真基因序列的變異是基于其父節(jié)點的基因序列，該網(wǎng)絡(luò)的生成規(guī)則為：

1) 給定病毒基因序列長度L，每位基因取值數(shù)量為n，假定每位基因的變異是互相獨立的，病毒基因在每一次傳播時序列發(fā)生變異的概率為p；

2) 給定病毒傳播能力值R，每一代傳播人數(shù)在以R為中心服從泊松分布(為計算方便，模擬生成網(wǎng)絡(luò)中傳染人數(shù)在0～6之間對稱分布).假定網(wǎng)絡(luò)中所有節(jié)點起始于“零號病人”，即有一個同一父節(jié)點，則根據(jù)R值，零號病人生成子節(jié)點，每個子節(jié)點的基因序列的變異是獨立的；

3) 重復(fù)步驟2)，直至傳播網(wǎng)絡(luò)節(jié)點數(shù)達到給定網(wǎng)絡(luò)規(guī)模N.

基于上述規(guī)則，本文生成病毒傳播仿真網(wǎng)絡(luò)，相關(guān)參數(shù)如表1所示:

Table 1 Parameters of Simulating Network表1 仿真網(wǎng)絡(luò)參數(shù)設(shè)置

真實數(shù)據(jù)集包括Aussie數(shù)據(jù)集和HIV數(shù)據(jù)集.Aussie數(shù)據(jù)集是部分澳大利亞新型冠狀病毒患者地基因序列，病人信息已經(jīng)匿名化處理，該數(shù)據(jù)集基因序列數(shù)據(jù)和標(biāo)準(zhǔn)參考序列比對出現(xiàn)多處缺失，序列質(zhì)量不高.我們使用TransPhylo對該基因序列集進行分析，TransPhylo是一個可以通過基因數(shù)據(jù)構(gòu)建病毒傳播網(wǎng)絡(luò)的工具[35].分析得到病毒傳播網(wǎng)絡(luò)如圖3所示，該網(wǎng)絡(luò)中共有1 032個節(jié)點，且約有10%的未采樣節(jié)點.網(wǎng)絡(luò)的每個節(jié)點代表一名新冠病毒患者，節(jié)點的有向連邊表示新冠病毒傳播關(guān)系，每個節(jié)點的屬性為其感染的新冠病毒的基因序列，基因長度為29 915.基因的每個位置的取值范圍為{A,T,C,G,-}，其中-表示該位基因缺失.如圖3所示，其中節(jié)點顏色越深，表示該節(jié)點距離零號病人越近.

Fig. 3 Aussie: 2019-nCoV transmission network圖3 Aussie新冠病毒傳播網(wǎng)絡(luò)

HIV數(shù)據(jù)集是加拿大阿爾伯塔省南部部分HIV患者的病毒序列信息，并且包含經(jīng)TransPhylo分析且經(jīng)醫(yī)院驗證的病毒傳播網(wǎng)絡(luò)信息.該網(wǎng)絡(luò)節(jié)點數(shù)量是140，病毒基因序列長度是1 522，該網(wǎng)絡(luò)中未采樣節(jié)點數(shù)約為41，即該傳播網(wǎng)絡(luò)包含較多未采樣節(jié)點，網(wǎng)絡(luò)如圖4所示.病毒傳播信息采集中，往往無法獲取0號病人的信息，該HIV數(shù)據(jù)集中同樣缺失0號病人基因序列信息.本文參考Lauren等人[36]研究推斷與0號病人傳播關(guān)系最接近的病毒序列，隨機改變約30%的序列信息后作為0號病人的病毒序列.

Fig. 4 HIV: HIV transmission network圖4 HIV病毒傳播網(wǎng)絡(luò)

3.3 驗證集數(shù)據(jù)生成

為了驗證基因表示模型的有效性，我們會在仿真和真實傳播網(wǎng)絡(luò)中刪除一些節(jié)點，這些刪除的節(jié)點即未采樣感染者，然后生成一個存在缺失節(jié)點的傳播網(wǎng)絡(luò)，存在缺失節(jié)點的連邊將作為驗證集，之后使用模型在該網(wǎng)絡(luò)上進行缺失節(jié)點預(yù)測.通過模型預(yù)測的結(jié)果和驗證集的比較，進而計算模型在缺失節(jié)點還原結(jié)果的準(zhǔn)確度.其中缺失網(wǎng)絡(luò)的生成規(guī)則為：

1) 生成指定數(shù)目缺失節(jié)點，檢查缺失節(jié)點合理性：若生成節(jié)點處于網(wǎng)絡(luò)兩端(沒有父節(jié)點或沒有子節(jié)點)，移除該節(jié)點；

2) 將缺失節(jié)點的父節(jié)點與子節(jié)點連接，生成新的連邊；

3) 重復(fù)步驟1)～2)，直至缺失節(jié)點數(shù)目達到給定值.

3.4 參數(shù)設(shè)置

對于基因編碼序列屬性，模型形式采用獨熱編碼；模型訓(xùn)練代數(shù)為1 000，初始學(xué)習(xí)速率lr=0.01，每進行100代訓(xùn)練，將學(xué)習(xí)率減半.表示向量維度為do=128，隱藏層維度為dh=1024，L2正則化參數(shù)β=0.001，注意力頭數(shù)nh=4，隨機數(shù)種子為1，隨機游走長度Lw=10，隨機游走重啟概率pr=20%.

3.5 對比算法

為了驗證GRGC模型對基因序列表示的高效性，本文采用5種算法作為基準(zhǔn)對比算法，包括不對基因序列降維處理直接用于計算的Direct方法與可以對高維基因序列進行降維及特征提取的PCA和Auto Encoder(AE)方法.此外，由于本文的方法是基于病毒傳播網(wǎng)絡(luò)構(gòu)建的，因此也選用了GCN和SDNE兩種網(wǎng)絡(luò)表示方法作為對比.

1) Direct方法.直接使用基因序列進行計算，每個節(jié)點的基因序列編碼后，不對其進行壓縮和特征提取，直接將其用于發(fā)現(xiàn)缺失節(jié)點任務(wù).

2) AE方法.自編碼器[37]由編碼器和解碼器2部分組成，其中編碼器能將輸入屬性壓縮成潛在的空間表征，解碼器能從潛在空間表征重構(gòu)數(shù)據(jù)的輸入.自編碼器常被用于降維或特征學(xué)習(xí)，本文首先使用AE模型對基因序列進行訓(xùn)練，然后使用訓(xùn)練好的AE模型的中間層即降維后的序列表示，作為發(fā)現(xiàn)缺失節(jié)點任務(wù)的輸入.

3) 主成分分析方法(principal component analysis， PCA)[38].是一種使用最廣泛的數(shù)據(jù)降維算法.PCA的主要思想是通過正交變化將可能存在相關(guān)性的N維特征數(shù)據(jù)映射到線性不相關(guān)的k維上，在原有N維特征的基礎(chǔ)上重新構(gòu)造k維特征，進而實現(xiàn)對原始數(shù)據(jù)的高效和精確的表示，轉(zhuǎn)換后的k維特征則被稱為主成分.本文使用基因序列的主成分作為發(fā)現(xiàn)缺失節(jié)點任務(wù)的輸入.

4) 圖卷積神經(jīng)網(wǎng)絡(luò)方法(graph convolutional neural network, GCN)[39].通過低通濾波器產(chǎn)生節(jié)點嵌入，將每個節(jié)點的信息傳遞到圖中的鄰居節(jié)點上去，通過卷積運算使得圖能夠逐層傳遞信息.本文構(gòu)建了一個2層GCN模型(隱藏層維度為256，輸出層維度為128)，并用上下文損失函數(shù)訓(xùn)練一個無監(jiān)督GCN模型，通過輸出層的基因表示，計算節(jié)點對相似度并發(fā)現(xiàn)缺失節(jié)點.

5) 結(jié)構(gòu)網(wǎng)絡(luò)嵌入方法(structural deep network embedding， SDNE)[40].使用深層神經(jīng)網(wǎng)絡(luò)對節(jié)點表示間的非線性進行建模，整個模型可以被分為2個部分：一個是由Laplace矩陣監(jiān)督的建模第1級相似度的模塊，另一個是由無監(jiān)督的深層自編碼器對第2級相似度關(guān)系進行建模.本文采用3層神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)(2層隱藏層維度分別為500和300，輸出層維度為128).

3.6 評價指標(biāo)

對所有對比算法和本文方法，本文通過AUC,RS和k-RS三個指標(biāo)評價其效果.

鏈路預(yù)測中常用AUC指標(biāo)[41]衡量算法的效果：每次隨機地從驗證集中選取1條存在的邊和1條不存在的邊并比較這2條邊的指標(biāo)值，如果存在的邊的分值大于不存在的邊的分值，則記為1分；如果2個邊的分?jǐn)?shù)值相等，則記為0.5分；若存在邊的分值小于不存在邊的分值，則記為0分.獨立重復(fù)地比較N次，如果有n′次驗證集中的邊的分?jǐn)?shù)值大于不存在的邊，有n″次驗證集中邊分?jǐn)?shù)值等于不存在邊的分?jǐn)?shù)，則AUC指標(biāo)計算結(jié)果為

(17)

AUC的指標(biāo)值越接近1，說明獲得的表示向量在連邊預(yù)測上的能力越強.

RS即排序分(ranking score)，考慮了驗證集中的邊在模型預(yù)測結(jié)果中排序的位置.令H=U-ET為未知邊的集合(相當(dāng)于驗證集中的邊和不存在的邊的集合)，re表示e∈EP在排序中的排名.那么這條驗證邊的排序分為

(18)

該數(shù)據(jù)集上的總的排序分為

(19)

根據(jù)目標(biāo)，本文希望存在未采樣感染者的連邊對的排名更靠后，因此，RS指標(biāo)越大，表明該表示向量在預(yù)測任務(wù)中表現(xiàn)越好.

k-RS統(tǒng)計RS得分前k條邊中驗證集邊所占的比例.如果模型預(yù)測結(jié)果中有m條邊在驗證集中，k-RS定義為

(20)

k-RS指標(biāo)側(cè)重于評價驗證集中特征最明顯的那部分邊，在實際應(yīng)用時，我們只會關(guān)心排名特征最明顯的部分(排名最靠前的部分)，從而在準(zhǔn)確率和覆蓋率中進行權(quán)衡，k-RS指標(biāo)不再關(guān)注本身特征就不明顯的邊，從而更好地反映和評價算法在真實使用時的效用.

4 實驗分析

本文分別在仿真病毒傳播網(wǎng)絡(luò)和Aussie新冠病毒傳播網(wǎng)絡(luò)、HIV病毒傳播網(wǎng)絡(luò)進行實驗.為了使實驗更符合真實場景，即現(xiàn)實中獲取的病毒傳播網(wǎng)絡(luò)往往存在未采樣感染者，因此我們首先在仿真或者真實數(shù)據(jù)集的基礎(chǔ)上生成存在一定節(jié)點缺失的網(wǎng)絡(luò)，然后在該存在缺失節(jié)點的網(wǎng)絡(luò)上訓(xùn)練基因表示模型.

在仿真數(shù)據(jù)實驗中，我們首先將模型在生成的病毒傳播網(wǎng)絡(luò)上進行訓(xùn)練，然后使用節(jié)點缺失比例為10%和20%的傳播網(wǎng)絡(luò)對模型進行測試.節(jié)點缺失網(wǎng)絡(luò)按照3.3節(jié)缺失網(wǎng)絡(luò)生成規(guī)則生成，即選取未采樣節(jié)點并將其在數(shù)據(jù)集中刪除.然后通過分析模型成功找出缺失節(jié)點的正確率來評價算法的表現(xiàn).在Aussie數(shù)據(jù)上，測試網(wǎng)絡(luò)的缺失節(jié)點比例分別為10%和20%.在HIV傳播網(wǎng)絡(luò)實驗中，由于病毒傳播網(wǎng)絡(luò)規(guī)模較小，10%的節(jié)點缺失總數(shù)是14，較小的驗證集會導(dǎo)致測試結(jié)果波動較大，因此在HIV實驗中我們使用全網(wǎng)絡(luò)進行模型訓(xùn)練，在30%節(jié)點缺失的網(wǎng)絡(luò)上進行測試.

4.1 結(jié)果分析

本節(jié)實驗提到的k-RS指標(biāo)中，百分?jǐn)?shù)值表示排名靠前相應(yīng)比例的連邊數(shù)(25%表示RS得分排名前25%的連邊)，Top 20表示最可能包含缺失節(jié)點的20條連邊.表中實驗結(jié)果均是10次生成缺失節(jié)點實驗的平均值，分值越高代表算法預(yù)測存在缺失節(jié)點的正確率越高.

從表2和表3可以看出，在仿真數(shù)據(jù)中，采用Direct方法的基因表示在進行缺失節(jié)點發(fā)現(xiàn)任務(wù)中的表現(xiàn)最好，這是因為仿真數(shù)據(jù)的基因序列生成及變異是嚴(yán)格按照規(guī)則生成，因此仿真基因序列中無效序列較少，對基因序列進行降維特征提取反而丟失了重要的特征.但由實驗結(jié)果可見，GRTN方法在仿真數(shù)據(jù)集上的表現(xiàn)依然優(yōu)于Direct以外的對比方法.

Table 2 Results on Simulating Network with 10% Nodes Missing

Table 3 Results on Simulating Network with 20% Nodes Missing

為探究本文方法在實際應(yīng)用中的價值，我們也在真實新冠病毒網(wǎng)絡(luò)中進行了實驗.分析表4和表5可得，在真實的新冠病毒病毒基因序列數(shù)據(jù)集中，采用Direct進行基因表示的實驗結(jié)果要遠(yuǎn)低于本文提出的GRTN方法.而GRTN方法通過聚合節(jié)點本身的基因序列信息和傳播網(wǎng)絡(luò)的特征信息，在多數(shù)驗證指標(biāo)上取得了最優(yōu)的成績，在Top-5和RS等指標(biāo)上結(jié)果略差于GCN方法，但其綜合表現(xiàn)在所有方法中最好.

Table 4 Results on Aussie with 10% Nodes Missing表4 缺失10%節(jié)點Aussie中各指標(biāo)表現(xiàn) %

Table 5 Results on Aussie with 20% Nodes Missing表5 缺失20%節(jié)點Aussie中各指標(biāo)表現(xiàn) %

為探究本文模型的適用性，本文選取HIV病毒傳播網(wǎng)絡(luò)進行試驗.因為HIV病毒基因網(wǎng)絡(luò)節(jié)點數(shù)目相對較少，只關(guān)注排名靠前的節(jié)點的指標(biāo)時(比如k-RS指標(biāo)的前3%)，實驗結(jié)果波動極大，由表6中結(jié)果可得，雖然HIV的病毒序列特征與COVID-19病毒的極為不同，即HIV的序列長度要遠(yuǎn)遠(yuǎn)小于COVID-19，且該網(wǎng)絡(luò)本身即包含大量缺失節(jié)點，但從實驗結(jié)果可以看出，基于病毒傳播網(wǎng)絡(luò)的基因表示模型(GRTN)在處理缺失節(jié)點發(fā)現(xiàn)的任務(wù)上，表現(xiàn)要遠(yuǎn)遠(yuǎn)好于其他對比方法.

Table 6 Results on HIV with 30% Nodes Missing表6 缺失30%節(jié)點HIV數(shù)據(jù)集中各指標(biāo)表現(xiàn) %

4.2 模型分析

4.1節(jié)實驗驗證了GRTN模型的有效性，本節(jié)從表示聚合分模塊能力和注意力使用等多個方面分析模型的表現(xiàn)效果.

4.2.1 表示聚合分析

在GRTN模型中，節(jié)點的基因表示由Ef,Eb和Em經(jīng)過注意力聚合得到，在本節(jié)，用GRTN-f表示去掉目標(biāo)節(jié)點上文節(jié)點信息聚合的模型，GRTN-b表示去掉下文節(jié)點信息聚合的模型，用GRTN-cut表示對輸入基因序列不進行維度約減的操作，用GRTN-2n表示不考慮目標(biāo)節(jié)點距離為2的鄰居信息的模型.綜合上述多個角度，對方法進行消融實驗.

從表7結(jié)果可知，本文模型通過結(jié)合前向、后向和自身的表示信息表現(xiàn)效果最佳，缺失前向表示的信息對模型的效果影響更大，后向節(jié)點表示信息對模型表現(xiàn)也有較大影響，單純使用節(jié)點自身表示的退化模型，效果較差.當(dāng)只聚合目標(biāo)節(jié)點的子節(jié)點和父節(jié)點信息時，模型效果較差，這說明只使用一階鄰居的信息進行基因表示是不充分的.當(dāng)不對輸入數(shù)據(jù)維度進行約減時，模型效果有明顯的下降，這說明模型對基因序列降維后進行了有效特征提取，相比不降維的全序列數(shù)據(jù)，可以去掉基因序列數(shù)據(jù)中無用的干擾信息.而對比同樣是對數(shù)據(jù)進行降維提取特征的PCA和AE，說明基于傳播網(wǎng)絡(luò)上下文的基因表示模型效果更優(yōu).

Table 7 Ablationstudy on Aussie with 20% Nodes Missing表7 缺失20%節(jié)點Aussie中各模塊消融實驗 %

4.2.2 注意力權(quán)重分析

在Aussie中，本文提取節(jié)點的Ef、Eb和Em向量的注意力權(quán)重，統(tǒng)計如表8所示.對于處于網(wǎng)絡(luò)不同位置、屬性信息不同的節(jié)點，其注意力系數(shù)有明顯差異.從整體影響力來說，前向注意力的系數(shù)權(quán)重比較高，這和4.2.1節(jié)中前向節(jié)點影響大的結(jié)論一致.圖5展示了網(wǎng)絡(luò)中節(jié)點7和節(jié)點621的局部結(jié)構(gòu)，以節(jié)點7為例，該節(jié)點處于網(wǎng)絡(luò)中的末端位置，沒有后向節(jié)點，其屬性結(jié)果很大程度上和其父節(jié)點相關(guān)，故αf高，雖然模型給出了αb的系數(shù)，但是該節(jié)點的Eb為0；節(jié)點621同時存在子節(jié)點和父節(jié)點，注意力系數(shù)和三者皆有關(guān)系.

Table 8 Values of Forward, Backward and Self Attention表8 前向、后向和自身注意力系數(shù)

Fig. 5 Nodes attention analysis圖5 節(jié)點注意力分析

4.2.3 多頭注意力分析

為了探討不同注意力機制使用對模型效果的影響，列出采用不同注意力機制的模型在Aussie中的表現(xiàn)結(jié)果，如圖6所示.本節(jié)使用的不同注意力機制有：

1) average-4.采用4頭注意力機制，每個注意力的表示維度為128，將4個注意力的表示向量輸出在每個維度上取平均值.

2) concat-2.采用2頭注意力機制，每個注意力的表示維度為128，將2個注意力的表示向量輸出拼接在一起.

3) concat-4.采用4頭注意力機制，每個注意力的表示維度為128，將4個注意力的表示向量輸出拼接在一起.

4) single.單頭注意力表示，輸出的表示向量維度為128.

Fig. 6 Model performances with different attention mechanisms on Aussie with 20% nodes missing圖6 采用不同注意力機制的算法在Aussie上預(yù)測準(zhǔn)確率表現(xiàn)(20%節(jié)點缺失)

從圖6分析可知，采用單頭注意力的模型表示結(jié)果的隨機性較大，采用拼接操作和取均值操作均可以提高模型的效果，其中采用拼接操作的效果更優(yōu)，并且隨著拼接的注意力模塊的增加，模型表現(xiàn)有一定提高.

5 結(jié) 論

本文基于病毒傳播網(wǎng)絡(luò)構(gòu)建了一種新的基因序列表示方法，在對節(jié)點病毒序列進行編碼后，使用注意力機制對節(jié)點及其鄰居節(jié)點的序列信息進行聚合，進而得到目標(biāo)節(jié)點病毒基因序列的表示.為了驗證基因序列表示模型的有效性，本文設(shè)計一套在病毒傳播網(wǎng)絡(luò)中發(fā)現(xiàn)未采樣感染者的任務(wù)流程，可以在實際的醫(yī)學(xué)調(diào)查診斷的相應(yīng)環(huán)節(jié)進行輔助排查.本文用網(wǎng)絡(luò)缺失點發(fā)現(xiàn)任務(wù)對得到的基因序列表示進行評估，模型在模擬病毒傳播網(wǎng)絡(luò)、真實的澳大利亞新冠病毒傳播網(wǎng)絡(luò)和艾滋病毒網(wǎng)絡(luò)上的實驗結(jié)果驗證了模型的有效性.

作者貢獻聲明：馬揚和劉澤一提出論文理論方法、思路，進行實驗及論文撰寫，為本文同等貢獻第一作者；程光權(quán)和劉忠指導(dǎo)論文思路與撰寫，為本文共同通信作者；梁星星、陽方杰、成清對論文提出修改建議.