李鑫，范虹，趙興春，范曉諾，姚若俠

技術(shù)與方法

基于全局最小殘差法快速分析混合STR圖譜

李鑫1，范虹1，趙興春2，范曉諾3，姚若俠1

1. 陜西師范大學(xué)計(jì)算機(jī)科學(xué)學(xué)院，西安 710119 2. 公安部鑒定中心，北京 100038 3. 波士頓大學(xué)生物系，波士頓 02215，美國

在法醫(yī)DNA分析領(lǐng)域，混合短串聯(lián)重復(fù)序列(short tandem repeats，STR)圖譜的分析一直是研究難點(diǎn)。當(dāng)前，國內(nèi)主要依靠法醫(yī)進(jìn)行人工分析，不僅效率低下，分析結(jié)果還存在著主觀性偏好，難以滿足日益增長的STR圖譜分析的需求。本文提出一種新的混合STR圖譜分析方法——全局最小殘差法，不僅可以計(jì)算出分析結(jié)果，還可以預(yù)測(cè)出每個(gè)組分的混合比例。該方法首先給混合比例賦予了新的定義，然后對(duì)等位基因模型進(jìn)行優(yōu)化，進(jìn)而綜合考慮STR圖譜中的所有基因座，將每個(gè)基因座的殘差值進(jìn)行累加求和，選擇累加和最小的混合比例作為推斷結(jié)果，并使用灰狼優(yōu)化算法快速尋找混合比例的最優(yōu)值。對(duì)于二組分STR圖譜，全局最小殘差法能夠兼顧分析的準(zhǔn)確性和分析速度，有利于實(shí)現(xiàn)大量的圖譜分析。本文提出的算法在實(shí)際應(yīng)用中取得了不錯(cuò)的效果，具有較高的應(yīng)用價(jià)值，可為混合STR圖譜分析領(lǐng)域的研究提供新的解決方案。

法醫(yī)遺傳學(xué)；混合STR圖譜；等位基因模型；灰狼優(yōu)化算法

法醫(yī)DNA分析最常用的遺傳標(biāo)記是短串聯(lián)重復(fù)序列(short tandem repeats，STR)[1,2]，當(dāng)DNA樣本來源于兩個(gè)或兩個(gè)以上個(gè)體時(shí)便可得到混合STR圖譜。關(guān)于混合STR圖譜的分析和解釋存在諸多挑戰(zhàn)，如STR分析結(jié)果解釋方法沒有統(tǒng)一的標(biāo)準(zhǔn)、混合樣本的供體個(gè)數(shù)難以確定以及不平衡混合樣本中的次要成分檢測(cè)困難等[3]。由于這些問題始終沒有得到有效解決，這使得分析混合STR圖譜變得困難重重。目前，國內(nèi)大部分實(shí)驗(yàn)室對(duì)于混合STR圖譜的分析都是人工拆分[4]。人工拆分需要分析人員有扎實(shí)的專業(yè)知識(shí)和實(shí)踐經(jīng)驗(yàn)，并且人工拆分時(shí)沒有量化的標(biāo)準(zhǔn)，主觀因素較大，導(dǎo)致分析結(jié)果不確定性較大。而國際上已有專門的拆分軟件，這些軟件拆分的優(yōu)勢(shì)在于：有可量化的標(biāo)準(zhǔn)、可重復(fù)性好、減少主觀性，結(jié)果更加客觀公正[5]。我國已經(jīng)有幾家法醫(yī)DNA實(shí)驗(yàn)室引進(jìn)了國外這些成熟的混合STR圖譜分析系統(tǒng)，如STRmix[6]。STRmix等軟件是基于大量數(shù)據(jù)訓(xùn)練完成的數(shù)學(xué)模型，由于不同國家和地區(qū)人種基因頻率存在差異，且國外和我國的法律條文及法律思想也存在差異，所以STRmix等軟件在國內(nèi)的應(yīng)用性十分有限。而且這些軟件價(jià)格昂貴，這些因素導(dǎo)致STRmix等軟件難以大規(guī)模推廣使用。季現(xiàn)超等[19]自主研發(fā)了一款基于馬爾科夫鏈蒙特卡洛(Markov Chain Monte Carlo，MCMC)算法的混合STR圖譜分析系統(tǒng)——SMART(STR Mixture Analysis and Resolution Tools)，該系統(tǒng)考慮了stutter峰、DNA降解和基因座特異性擴(kuò)增效率等因素，能夠分析現(xiàn)場物證混合STR圖譜，為案件偵破提供線索。但是SMART目前尚未在國內(nèi)推廣使用，也沒有其他適合國內(nèi)應(yīng)用環(huán)境的成熟軟件。當(dāng)前受國際形勢(shì)的影響，STRmix已被列入對(duì)華禁售的高科技產(chǎn)品清單，國內(nèi)需要進(jìn)行分析的混合STR圖譜數(shù)量日益增加，我國亟需自主研發(fā)與STRmix一樣成熟且適用于我國應(yīng)用環(huán)境的分析軟件。當(dāng)檢材供體數(shù)增加時(shí)，STR圖譜分析的難度也會(huì)加大[7]，即使是當(dāng)前國際流行的STR圖譜分析軟件，在分析速度和準(zhǔn)確性上也難以兩面兼顧。所以，有必要建立分析結(jié)果更穩(wěn)健且運(yùn)算速度更快的優(yōu)化模型，并在實(shí)踐中檢驗(yàn)其魯棒性。Gill等[8]建立了利用峰面積分析二組分混合樣本的方法，該方法充分利用STR圖譜的定性信息(等位基因)和定量信息(峰面積)建立等位基因模型(allele model)，使用最小殘差法(the minimum residuals method)推斷混合樣本中各組分的基因型，但是這個(gè)方法的計(jì)算量較大，不能實(shí)現(xiàn)對(duì)STR圖譜的高效分析。鑒于此，本文提出全局最小殘差法，對(duì)等位基因模型進(jìn)行優(yōu)化，并引入灰狼優(yōu)化算法，提高混合STR圖譜的分析效率。

1 材料與方法

1.1 STR圖譜

STR圖譜又叫STR分型結(jié)果峰圖，混合STR圖譜就是兩名或兩名以上個(gè)體混合在一起的DNA樣本檢測(cè)生成的STR圖譜。通過STR圖譜可以獲取所檢測(cè)生物樣本的基因信息，STR圖譜上的等位基因的峰面積是對(duì)檢測(cè)樣本中的DNA模板量的映射，法醫(yī)實(shí)驗(yàn)已證實(shí)，DNA的模板量與STR圖譜中的等位基因峰面積存在線性關(guān)系[8]，即檢測(cè)樣本中的DNA模板量越大，該檢測(cè)樣本生成的STR圖譜上的等位基因的峰面積越大。基于此原理，分析者可以通過觀察到的峰面積來推斷可能的基因型組合和各組分的比例。

對(duì)于任意兩個(gè)人的混合樣本，一個(gè)基因座(檢測(cè)位點(diǎn))上的等位基因個(gè)數(shù)最大為4，即單個(gè)基因座上的峰型可以表現(xiàn)為一帶型、二帶型、三帶型或四帶型，它們對(duì)應(yīng)可能存在的基因型組合數(shù)分別為1、7、12、6，詳見表1[9]。極少數(shù)個(gè)體會(huì)發(fā)生染色體變異，細(xì)胞核DNA表現(xiàn)為三倍體[5]，這時(shí)，二組分混合STR圖譜基因座上的等位基因個(gè)數(shù)可能會(huì)超過4個(gè)，這種情況本文不作討論?；蛐徒M合是指兩個(gè)供體基因型的組合，例如基因型組合()表示供體1是基因型為，另一個(gè)供體2的基因型為。二組分混合STR圖譜每個(gè)基因座上的分析結(jié)果在表1中均可以找到，因此表1可視為二組分STR圖譜分析的解空間。

1.2 等位基因模型

分析已知供體數(shù)的混合STR圖譜的第一步是估計(jì)混合比例(mixture proportion，Mx)，混合比例是指混合樣本中某一供體DNA模板量所占的比例。法醫(yī)學(xué)實(shí)驗(yàn)證明混合STR圖譜的所有基因座上的混合比例都是一樣的[8]。所以一旦求出混合比例，就可以建立等位基因模型來估計(jì)表1中基因型組合的期望峰面積。事實(shí)上，當(dāng)基因座上只有一個(gè)等位基因時(shí)，可能的基因型組合只有一個(gè)，即兩個(gè)體的基因型均為純合子，因此一帶型基因座可以直接得到分析結(jié)果，本文對(duì)一帶型基因座不做分析；當(dāng)基因座上有四個(gè)等位基因時(shí)，此時(shí)每一種基因型組合都不存在共享等位基因，估計(jì)出混合比例比較容易做到的。若四帶型基因座的基因型組合為,，即個(gè)體1的基因型為，另一個(gè)供體的基因型為，則可以通過公式(1)計(jì)算出混合比例的一個(gè)估計(jì)值′。

其中，P、P、P、P分別表示四帶型基因座上的等位基因峰面積。

Gill等[8]對(duì)二組分混合樣本建立的等位基因模型如附表1所示?？梢园l(fā)現(xiàn)，如果基因座上有兩個(gè)或三個(gè)的等位基因，這個(gè)過程就不再像公式(1)那樣簡單。二帶型、三帶型和四帶型基因座的基因型組合數(shù)分別為7、12、6，那么后面使用最小殘差法進(jìn)行分析時(shí)必然會(huì)涉及大量計(jì)算。因此，本文對(duì)等位基因模型進(jìn)行優(yōu)化，較原來的等位基因模型可以減少一些計(jì)算量。

在優(yōu)化等位基因模型前，本文重新定義混合比例：混合樣本中DNA模板量最少的供體所占的比例。那么在二組分混合樣本中，假設(shè)個(gè)體1的DNA模板量所占比例為，個(gè)體2的DNA模板量所占比例為1–，此時(shí)的取值范圍由(0, 100%)變?yōu)?0, 50%]。

首先將每個(gè)基因座上的等位基因按照其峰面積大小降序排列，并對(duì)等位基因進(jìn)行標(biāo)記。

對(duì)于四帶型基因座，等位基因分別標(biāo)記為,,,，即P≥P≥P≥P；

對(duì)于三帶型基因座，等位基因分別標(biāo)記為,,，即P≥P≥P；

對(duì)于二帶型基因座，等位基因分別標(biāo)記為,，即P≥P。

基于以上的內(nèi)容可以構(gòu)造不等式組。

對(duì)于四帶型基因座則有

表1 二組分混合樣本可能的基因型組合

*,,,分別代表基因座上的等位基因。

對(duì)于三帶型基因座則有

對(duì)于二帶型基因座則有

其中P、P、P、P分別為等位基因,,,的期望峰面積。

將表1中每一種基因型組合代入公式(2) (3) (4)中求得基因型組合對(duì)應(yīng)的取值范圍，舍去不等式組無解的基因型組合便可得到優(yōu)化的二組分混合檢材可能的基因型組合，如表2所示。因?yàn)闊o需對(duì)一帶型基因座進(jìn)行分析，所以這里不再將其列入組合表中。對(duì)比表2，可以發(fā)現(xiàn)二帶型、三帶型和四帶型等位基因的基因型組合數(shù)分別減小為4、6、1。根據(jù)表2建立新的等位基因模型，如表3所示。那么最小殘差法基于優(yōu)化的等位基因模型進(jìn)行圖譜分析時(shí)，二組分STR圖譜分析的解空間就變成了表2里的11種基因型組合，相較原來的等位基因模型計(jì)算量減少了近一半。此時(shí)用枚舉法求解混合比例時(shí)，的取值范圍由(0, 100%)縮減為(0, 50%]。在更小的解空間里，分析結(jié)果更不易出錯(cuò)，最小殘差法的本質(zhì)是在解空間里遍歷搜索尋找最優(yōu)的解。本文對(duì)等位基因模型進(jìn)行優(yōu)化相當(dāng)于是對(duì)二組分STR圖譜分析的解空間的優(yōu)化，讓最小殘差法在優(yōu)化的解空間里搜索，一定程度上也實(shí)現(xiàn)了對(duì)最小殘差法的優(yōu)化。

Fan等[20]也對(duì)等位基因模型進(jìn)行了優(yōu)化，但是他們的三帶型等位基因模型只有(,)、(,)、(,)，而本文的三帶型等位基因模型(如表3所示)的基因型組合中有(,)、(,)、(,)、(,)、(,)、(,)，之所以要保留后三者是處于數(shù)學(xué)角度的考量，因?yàn)榛蛐徒M合(,)、(,)、(,)在公式(3)中均可求解到的取值范圍，不能排除不予考慮。

表2 二組分混合檢材可能的基因型組合(優(yōu)化版)

,,,分別代表基因座上的有效等位基因。

表3 等位基因模型(優(yōu)化版)

1.3 全局最小殘差法

最小殘差法(the minimum residuals method)是基于STR圖譜的定性信息(等位基因)和定量信息(峰面積)的一種分析方法。由表3計(jì)算出等位基因峰面積的期望值，如果基因型組合推斷正確，等位基因的峰面積期望值與觀察值應(yīng)該是近似的[8]。最小殘差法使用殘差這個(gè)變量來量化峰面積期望值與觀察值的近似程度，如公式(5)所示。越小，當(dāng)前基因型組合為真實(shí)結(jié)果的概率越大。

其中，residual表示STR圖譜中第個(gè)基因座在基因型組合為時(shí)的殘差值，是基因座的等位基因數(shù)量，P、P分別表示等位基因的峰面積觀察值與期望值。此時(shí)基因座的推測(cè)基因型組合combination可通過公式(6)求得。

其中，是基因型組合變量，它的取值空間為表2中的11種基因型組合，combination表示使得residual取值最小時(shí)的基因型組合，即=combination時(shí)，residual取值最小。

基因座的殘差值的最小值用residual表示，如公式(7)所示。

Gill等[8]提出的分析方法是對(duì)每個(gè)基因座單獨(dú)拆分，在的所有取值中分別計(jì)算所有基因型組合(詳見表1)的，取使得最小時(shí)的基因型組合作為拆分結(jié)果。需要注意的是，這樣的策略忽視了整個(gè)STR圖譜的，因此，本文提出全局最小殘差法，具體方法如下：

在(0, 50%]范圍中取值，尋找使得整個(gè)STR圖譜的所有基因座的最小殘差值之和residual最小的數(shù)值，并把這個(gè)數(shù)值賦值給，如公式(8)和(9)所示。再把代入等位基因模型(表3)中，根據(jù)公式(5)分別計(jì)算每個(gè)基因座的所有基因型組合(表2)的，根據(jù)公式(6)(7)計(jì)算每個(gè)基因座的殘差的最小值residual，取殘差=residual時(shí)的基因型組合combination作為分析結(jié)果。

其中，為STR圖譜中基因座的數(shù)量，residual為整個(gè)STR圖譜中所有基因座的殘差最小值residual之和，公式(9)表示使得residual取值最小時(shí)的混合比例，即=時(shí)，residual取值最小。residual越小，當(dāng)前越接近真實(shí)的混合比例，由公式(6)得出的分析結(jié)果正確的概率越大。

1.4 灰狼優(yōu)化算法

灰狼優(yōu)化算法(grey wolf optimizer，GWO)是受灰狼群體捕食行為的啟發(fā)，由Mirjalili等[10]于2014年提出的一種新型群體智能優(yōu)化算法。該算法具有實(shí)現(xiàn)簡單、收斂速度快的優(yōu)勢(shì)，目前被廣泛應(yīng)用于聚類分析[11]、路徑規(guī)劃[12]、圖像分割[13]和參數(shù)優(yōu)化等諸多領(lǐng)域。

研究發(fā)現(xiàn)，在灰狼群體中有著嚴(yán)格的等級(jí)制度，一小部分擁有絕對(duì)權(quán)利的灰狼領(lǐng)導(dǎo)狼群進(jìn)行捕獵活動(dòng)[14]。一般來說，灰狼社會(huì)有4個(gè)等級(jí)，如圖1所示，、、和分別代表4個(gè)社會(huì)等級(jí)。按照這樣的等級(jí)制度，等級(jí)高的灰狼對(duì)等級(jí)低的灰狼有絕對(duì)的支配權(quán)。即，灰狼對(duì)灰狼、和有絕對(duì)的支配權(quán)；灰狼對(duì)灰狼和有絕對(duì)的支配權(quán)；灰狼對(duì)灰狼有絕對(duì)的支配權(quán)。為了模擬灰狼的捕獵行為，假設(shè)、、更接近目標(biāo)獵物的位置(即適應(yīng)度函數(shù)的最優(yōu)解)。另外，灰狼的數(shù)量在灰狼群中的比例最大，同時(shí)灰狼又必須完全服從灰狼、和，所以灰狼群的捕獵行為主要由灰狼、和引領(lǐng)完成。在每次迭代結(jié)束時(shí)，記錄3個(gè)最優(yōu)解作為下一次迭代時(shí)灰狼、、的位置。在灰狼優(yōu)化算法中，為了模擬灰狼種群的等級(jí)制度，同時(shí)又能簡化算法，通常在算法中假設(shè)各有一只灰狼、一只灰狼，一只灰狼和若干只灰狼。

圖1 灰狼社會(huì)等級(jí)

假設(shè)X,X,X分別表示灰狼,,的當(dāng)前位置，X為灰狼的當(dāng)前位置(=,,,)。則灰狼在灰狼(=,,)的引領(lǐng)下的捕獵行為定義如下：

其中，D表示灰狼在灰狼之間的距離，X表示灰狼在灰狼的引領(lǐng)下的更新位置，A和C是系數(shù)向量，是收斂因子，它隨著迭代次數(shù)從2線性遞減到0，1,2為0到1之間的隨機(jī)數(shù)。

那么灰狼在灰狼,,共同引領(lǐng)下的更新位置為

其中，是當(dāng)前的迭代次數(shù)，X(+1)表示灰狼在下次迭代時(shí)的位置。

1.5 GWO-全局最小殘差法

本文在全局最小殘差法中使用灰狼優(yōu)化算法，可以快速找到使得residual取值最小時(shí)的混合比例的估值，即的值。這時(shí)可以認(rèn)為是最接近真實(shí)混合比例的取值，在這個(gè)取值下獲得的拆分結(jié)果正確的概率最大。

對(duì)于公式(9)中的，最簡單的計(jì)算方法是在0到50%的范圍依次取1%，2%，…，49%，50%分別代入公式(8)計(jì)算residual，找到使得residual最小的。如果使用灰狼優(yōu)化算法則能夠快速計(jì)算出這個(gè)結(jié)果。那么，本文把混合比例的值等價(jià)為灰狼優(yōu)化算法中灰狼的位置，公式(8)作為灰狼優(yōu)化算法的適應(yīng)度函數(shù)，如公式(12)所示。

結(jié)合混合STR圖譜分析的具體情況，對(duì)灰狼優(yōu)化算法進(jìn)行如下參數(shù)初始化。

灰狼數(shù)量=4，的取值太大算法運(yùn)行時(shí)間會(huì)變長，取值太小不能保證準(zhǔn)確率，4是在保證準(zhǔn)確率的情況下的最優(yōu)數(shù)值；因?yàn)?<≤50%，所以灰狼的初始位置取值范圍為(0, 0.5]；迭代數(shù)設(shè)置為6，多次實(shí)驗(yàn)發(fā)現(xiàn)灰狼優(yōu)化算法經(jīng)過6次迭代后已收斂至適應(yīng)度函數(shù)的最優(yōu)值。

在全局最小殘差法中引入灰狼優(yōu)化算法(GWO-全局最小殘差法)的算法描述如圖2所示。

輸入：二組分混合STR圖譜(基因座，等位基因，峰面積)過程：初始化灰狼算法：狼群數(shù)量n=4, Mx∈(0,0.5], 迭代次數(shù)iterations=6；初始化灰狼i(i =1,2,…,n)的位置Xi；fort = 1 to iterationsdo 根據(jù)公式(12)計(jì)算每個(gè)灰狼的適應(yīng)度； Xα為當(dāng)前適應(yīng)度函數(shù)F第一最優(yōu)解； Xβ為當(dāng)前適應(yīng)度函數(shù)F第二最優(yōu)解； Xδ為當(dāng)前適應(yīng)度函數(shù)F第三最優(yōu)解； 由公式(10)(11)計(jì)算灰狼群體i(i =1,2,…,n)的更新位置Xi(t+1)； 更新a、A和C；end forMx′=Xα；Mx’代入公式(12)求得residualsum輸出：Mx', residualsum

2 結(jié)果與分析

2.1 混合STR圖譜的全局最小殘差值

實(shí)驗(yàn)數(shù)據(jù)來自文獻(xiàn)[15](詳見附表 2)，由R. Wickenheiser的Crime實(shí)驗(yàn)室提供，本文把該數(shù)據(jù)標(biāo)記為Data 1。該STR圖譜包含有13個(gè)基因座，同時(shí)還給出兩個(gè)供體的真實(shí)基因型和每個(gè)基因座的等位基因?qū)?yīng)的峰面積。因此可以使用該數(shù)據(jù)驗(yàn)證全局最小殘差法的分析效果。

全局最小殘差法對(duì)混合比例估計(jì)值的計(jì)算可以使用枚舉法，即依次取值1%，2%，…，49%，50%，分別計(jì)算residual，取使residual最小值時(shí)的作為最終的估計(jì)值，然后在混合比例為時(shí)的拆分結(jié)果作為最終的分析結(jié)果。

首先，將Data 1中的每個(gè)基因座的等位基因按照峰面積的大小降序排序，然后進(jìn)行歸一化。這一步會(huì)把低于主峰峰面積15%的峰形視為stutter峰，并從圖譜中刪除不參與后面的計(jì)算。實(shí)驗(yàn)數(shù)據(jù)的歸一化結(jié)果如表4所示，歸一化結(jié)果對(duì)應(yīng)為公式(5)中的P。

然后，將依次賦值0.01，0.02，…，0.49，0.50，計(jì)算residual的值。這里以=0.20為例(即個(gè)體1的DNA模板量在混合樣本中所占比例為20%)。參考表3可以得到=20%時(shí)的三等位基因模型如表5所示。

表4 Data 1的歸一化結(jié)果

“等位基因”列上的數(shù)字代表基因座上短串聯(lián)重復(fù)序列的重復(fù)數(shù)。

因此，當(dāng)=0.2時(shí)，基因座的推斷分析結(jié)果為(,)，即(16/16，15/18)。

同理可以得到其他三帶型等位座在=0.20時(shí)的最小殘差值及分析結(jié)果，如表7所示。

同理可以計(jì)算出四帶型基因座和二帶型基因座在=0.20時(shí)的最小殘差值及分析結(jié)果，如此一來，便可以得到整個(gè)圖譜在=0.20時(shí)的分析結(jié)果，如表8所示。

根據(jù)上述步驟，可以得到取其他值時(shí)的residual。根據(jù)公式(9)，存在一個(gè)取值,使得residual最小。本文把=時(shí)的分析結(jié)果作為最終的分析結(jié)果。

依次取0.01，0.02，…，0.49，0.50進(jìn)行上述計(jì)算，可以繪制出與residual關(guān)系圖像(如圖3所示)。圖3中的黑色圓點(diǎn)是residual最小值的坐標(biāo)位置(0.3, 0.030011)，該坐標(biāo)的縱坐標(biāo)0.030011即為全局最小殘差值；當(dāng)取值0.3時(shí)，residual最小，那么的估計(jì)值=30%，將=0.3時(shí)全局最小殘差法得出的分析結(jié)果為最終分析結(jié)果。

表5 Mx=20%時(shí)的三等位基因模型

表7 三帶型基因座在Mx=0.20時(shí)的最小殘差值及分析結(jié)果

*表示基因座的最小殘差值。

表8 STR圖譜在Mx=0.20時(shí)的最小殘差值及分析結(jié)果

*表示所有基因座的最小殘差值residual之和，即當(dāng)=0.2時(shí)，residual=0.133323。

2.2 對(duì)比實(shí)驗(yàn)

mixsep (Forensic Genetics DNA Mixture Sepa-rator)是由Tvedebrink等[16]用R語言編寫的用于法醫(yī)混合DNA分析的經(jīng)典軟件。為了對(duì)比全局最小殘差法的分析效果，本文使用最小殘差法和mixsep對(duì)Data 1進(jìn)行分析，三種方法的分析結(jié)果如表9所示。根據(jù)附表2中供體的真實(shí)基因型，在表9中用加粗字體標(biāo)注出分析有誤的結(jié)果。由表9可知全局最小殘差法與mixsep的分析結(jié)果全部正確，而最小殘差法在基因座和上卻給出了兩種分析結(jié)果，還需要人工進(jìn)一步判斷。

圖3 residualsum與Mx的關(guān)系圖

黑色圓點(diǎn)為該曲線的最低點(diǎn)。

為了驗(yàn)證灰狼優(yōu)化算法的優(yōu)化效果，本文使用python語言實(shí)現(xiàn)灰狼優(yōu)化算法，將其融入全局最小殘差法中尋找混合比例的最優(yōu)值。這里同樣使用Data 1進(jìn)行對(duì)比實(shí)驗(yàn)，表10從混合比例、residual最小值、分析準(zhǔn)確率以及分析耗時(shí)四個(gè)方面比較全局最小殘差法和GWO-全局最小殘差法。

由表10可知，全局最小殘差法在引入灰狼優(yōu)化算法后，得到的混合比例估計(jì)值和residual最小值，與全局最小殘差法的分析結(jié)果很接近，分析結(jié)果的準(zhǔn)確率依然是100%，而且減少了約60%的分析耗時(shí)。該對(duì)比實(shí)驗(yàn)證明，在全局最小殘差法引入灰狼優(yōu)化算法，可以在不影響分析結(jié)果的情況下，提高分析效率。圖4為灰狼優(yōu)化算法求解混合比例最優(yōu)值時(shí)適應(yīng)度函數(shù)的收斂曲線，由該圖可以觀察到灰狼優(yōu)化算法每次迭代所求得的residual的最小值，經(jīng)過6次迭代得到了residual的最小值0.0298612。

2.3 實(shí)踐應(yīng)用

為了驗(yàn)證全局最小殘差法在實(shí)踐應(yīng)用中的表現(xiàn)，本文使用一組來自國內(nèi)某法醫(yī)實(shí)驗(yàn)室的STR圖譜數(shù)據(jù)進(jìn)行實(shí)驗(yàn)(見附表3)，本文把該數(shù)據(jù)標(biāo)記為Data 2。Data 2是一份有21個(gè)基因座(檢測(cè)位點(diǎn))的二組分混合STR圖譜數(shù)據(jù)，由于數(shù)據(jù)保密性，本文只展示了部分的STR圖譜信息和供體基因型。

表9 最小殘差法、全局最小殘差法和mixsep分析結(jié)果比較

加粗字體表示該分析結(jié)果與供體的真實(shí)基因型不一致。

表10 全局最小殘差法與GWO-全局最小殘差法的對(duì)比實(shí)驗(yàn)

準(zhǔn)確率=分析正確的位點(diǎn)數(shù)/STR圖譜的位點(diǎn)數(shù)。

圖4 灰狼優(yōu)化算法收斂曲線

該數(shù)據(jù)經(jīng)過全局最小殘差法分析后求得混合比例的估值為32%。這里同樣使用最小殘差法和mixsep對(duì)Data 2進(jìn)行分析，這三種方法的分析結(jié)果如表11所示。可以發(fā)現(xiàn)，最小殘差法在基因座和上分析有誤，而全局最小殘差法和mixsep，除了個(gè)體1在基因座上的分析結(jié)果有誤，其余分析結(jié)果均與真實(shí)基因型一致。

表12列舉了基因座在=32%時(shí)四種基因型組合的期望峰形，可以發(fā)現(xiàn)全局最小殘差法選擇(7/9, 9/9)作為分析結(jié)果，而該基因座的真實(shí)基因型組合為(7/7, 9/9)，比較這個(gè)兩種結(jié)果的殘差值，可以發(fā)現(xiàn)兩者的殘差值很懸殊。圖5是基因座的四種基因型組合的殘差值與的關(guān)系圖，可以發(fā)現(xiàn)(7/9, 9/9)的曲線(紅)和(7/7, 9/9)的曲線(綠)的相交點(diǎn)橫坐標(biāo)在0.3附近，而全局最小殘差對(duì)混合比例的估值=32%也正好在臨界區(qū)域附近。

表11 實(shí)踐應(yīng)用的分析結(jié)果

加粗字體表示該分析結(jié)果與真實(shí)結(jié)果不一致。準(zhǔn)確率=分析正確的位點(diǎn)數(shù)/STR圖譜的位點(diǎn)數(shù)。

表12 Mx =32%時(shí)基因座TH01的期望峰形

生物檢材在檢測(cè)生成STR圖譜的過程中會(huì)受外界干擾從而產(chǎn)生峰形波動(dòng)，例如出現(xiàn)雜合不平衡、stutter峰、drop-in等[17]，這都會(huì)影響等位基因模型的擬合效果。相比最小殘差法，全局最小殘差法對(duì)的分析結(jié)果無誤，說明面對(duì)不可避免的干擾，全局最小殘差法穩(wěn)定性更好，且效果并不差于mixsep。由附表3可知，基因座和均為二帶型基因座，它們的真實(shí)基因型組合均為兩個(gè)純合子，雜合均衡比分別為0.4215、0.2892，那么說明兩個(gè)峰并不是完全來自一個(gè)供體。換而言之，這兩個(gè)基因座的峰形差異性較大，基因座受到的干擾使得它的峰形嚴(yán)重偏離等位基因模型上的期望峰形，錯(cuò)把其他期望峰形匹配為最相似的峰形。由圖5可以發(fā)現(xiàn)，如果可以降低這個(gè)干擾，那么紅綠兩條曲線的交點(diǎn)會(huì)向右偏移，這樣基因型組合(7/7, 9/9)將會(huì)是相似度最大的期望峰形，全局最小殘差法就能夠更正這個(gè)錯(cuò)誤。

圖5 基因座TH01的殘差值residual與混合比例Mx關(guān)系圖

3 討論

利用STR圖譜的定量信息峰面積對(duì)混合DNA樣本進(jìn)行分析，能夠?qū)崿F(xiàn)對(duì)混合STR圖譜的自動(dòng)化分析。本文提出的全局最小殘差法在引入灰狼優(yōu)化算法后，可以快速分析二組分混合DNA樣本生成的STR圖譜，該分析方法可以推斷出混合樣本中各組分的混合比例和基因型，具有較高的準(zhǔn)確率和較快的分析速度，對(duì)實(shí)現(xiàn)對(duì)混合STR圖譜的批量化分析具有重要意義，可以更好地協(xié)助法醫(yī)工作者進(jìn)行圖譜分析，提高STR圖譜的利用率。

混合STR圖譜受多種因素的影響，如DNA降解、峰飽和、等位基因共享、雜合子不平衡、基因座的特異性擴(kuò)增效率以及平行擴(kuò)增效率等[17]，這些因素都會(huì)造成等位基因模型難以很好的擬合峰面積觀察值。后續(xù)的研究工作可以針對(duì)影響峰面積的因素進(jìn)行數(shù)學(xué)建模，在進(jìn)行圖譜分析前，可以對(duì)圖譜數(shù)據(jù)進(jìn)行降噪處理，這樣可以最大程度還原數(shù)據(jù)。但是數(shù)據(jù)降噪處理的過程必然會(huì)增加計(jì)算量。如何對(duì)數(shù)據(jù)進(jìn)行降噪處理，增強(qiáng)等位基因模型的擬合性，以及在保證準(zhǔn)確分析的同時(shí)兼顧分析速度是后續(xù)研究的重點(diǎn)。另外，本文提出的方法是否可以有效地分析來自三個(gè)及以上的混合STR圖譜，仍有待實(shí)驗(yàn)檢驗(yàn)。混合STR圖譜分析本質(zhì)上是一個(gè)數(shù)學(xué)組合問題，由表1可知，對(duì)于任意兩個(gè)人的混合STR圖譜，基因座上等位基因的可能組合數(shù)最大為12，但是隨著供體數(shù)的增加，最大基因型組合數(shù)將指數(shù)式增長。比如供體數(shù)為3時(shí)，最大基因型組合數(shù)為294；供體數(shù)為4時(shí)，最大基因型組合數(shù)為16020。同時(shí)，多組分STR圖譜的等位基因模型推導(dǎo)將變得困難，這時(shí)引入智能優(yōu)化算法優(yōu)化分析過程依然有重要意義。而且供體數(shù)超過2時(shí)，STR圖譜發(fā)生stutter峰和drop-in的情況會(huì)增加[17]，需要人工完成數(shù)據(jù)預(yù)處理，或者直接丟棄該數(shù)據(jù)。

在法醫(yī)應(yīng)用中，往往期望從STR序列信息中獲得更多有效等位基因，這樣可以提高個(gè)體識(shí)別能力。目前二代測(cè)序技術(shù)趨近成熟[18]，能夠挖掘出混合DNA中更多的信息。相比毛細(xì)管電泳STR分型方法，二代測(cè)序技術(shù)的檢測(cè)靈敏度更高，對(duì)于一些陳舊疑難樣本有較高的檢出率。對(duì)于基因序列變異，二代測(cè)序技術(shù)能夠利用等位基因序列的差異，對(duì)主要組分等位基因、次要組分等位基因以及stutter峰加以區(qū)分，這更有利于解釋混合STR圖譜。因此，后續(xù)可以收集二代測(cè)序STR數(shù)據(jù)驗(yàn)證模型的穩(wěn)定性，再結(jié)合人工模擬的混合樣本驗(yàn)證推斷的混合比例是否存在較大誤差，從而設(shè)計(jì)更有效的混合STR圖譜分析方法，為法醫(yī)工作者提供更可靠的參考。

附加材料見文章電子版www.chinagene.cn。

[1] Cowell RG, Lauritzen SL, Mortera J. Identification and separation of DNA mixtures using peak area information., 2007, 166(1): 28–34.

[2] Gill P, Sparkes R, Kimpton C. Development of guidelines to designate alleles using an STR multiplex system., 1997, 89(3): 185–197.

[3] Bleka ?, Benschop CCG, Storvik G, Gill P. A comparative study of qualitative and quantitative models used to interpret complex STR DNA profiles., 2016, 25: 85–96.

[4] Budowle B, Onorato AJ, Callaghan TF, Manna AD, Gross AM, Guerrieri RA, Luttman JC, McClure DLMixture interpretation: defining the relevant features for guidelines for the assessment of mixed DNA profiles in forensic casework., 2009, 54(4): 810–821.

[5] Taylor D, Bright JA, Buckleton J. The interpretation of single source and mixed DNA profiles., 2013, 7(5): 516–528.

[6] Swaminathan H, Grgicak CM, Medard M, Lun DC. NOCIt: a computational method to infer the number of contributors to DNA samples analyzed by STR genotyping., 2015, 16: 172–180.

[7] Curran JM. A MCMC method for resolving two person mixtures., 2008, 48(4): 168–177.

[8] Gill P, Sparkes R, Pinchin R, Clayton T, Whitaker J, Buckleton J. Interpreting simple STR mixtures using allele peak areas., 1998, 91(1): 41–53.

[9] Tvedebrink T, Eriksen PS, Mogensen HS, Morling N. Identifying contributors of DNA mixtures by means of quantitative information of STR typing., 2012, 19(7): 887–902.

[10] Mirjalili S, Mirjalili SM, Lewis A. Grey wolf optimizer., 2014, 69(3): 46–61.

[11] Zhang S, Zhou YQ. Grey wolf optimizer based on powell local optimization method for clustering analysis., 2015, 2015: 1–17.

[12] Zhang S, Zhou YQ, Li ZM, Pan W. Grey wolf optimizer for unmanned combat aerial vehicle path planning., 2016, 99: 121–136.

[13] Khairuzzaman AKM, Chaudhury S. Multilevel threshol-ding using grey wolf optimizer for image segmentation., 2017, 86: 64–76.

[14] Zhang XF, Wang XY. Comprehensive review of grey wolf optimization algorithm., 2019, 46(3): 30–38.張曉鳳, 王秀英. 灰狼優(yōu)化算法研究綜述. 計(jì)算機(jī)科學(xué), 2019, 46(3): 30–38.

[15] Wang T, Xue N, Birdwell JD. Least-square deconvolution: a framework for interpreting short tandem repeat mixtures., 2006, 51(6): 1284–1297.

[16] Tvedebrink T. mixsep: An R-package for DNA mixture separation., 2011, 3(1): e486–e488.

[17] Clayton TM, Whitaker JP, Sparkes R, Gill P. Analysis and interpretation of mixed forensic stains using DNA STR profiling., 1998, 91(1): 55–70.

[18] B?rsting C, Morling N. Next generation sequencing and its applications in forensic genetics., 2015, 18: 78–89.

[19] Ji XC, Chi LJ, Xu Z, Peng Z, Ye J, Tu Z, Chen H. SMART: an analysis system for mixed str profiles., 2022, 47(1): 1–9.季現(xiàn)超, 池連江, 徐珍, 彭柱, 葉健, 凃政, 陳華. SMART：自主研發(fā)的混合STR圖譜分析系統(tǒng). 刑事技術(shù), 2022, 47(1): 1–9.

[20] Fan HL, Xie QQ, Wang LX, Ru K, Tan XH, Ding JY, Wang X, Huang J, Wang Z, Li YN, Wang XH, He YT, Gu CH, Liu M, Ma SW, Wen SQ, Qiu PM. Microhaplotype and Y-SNP/STR (MY): A novel MPS-based system for genotype pattern recognition in two-person DNA mixtures., 2022, 59: 102705.

Rapid analyzing mixed STR profiles based on the global minimum residual method

Xin Li1, Hong Fan1, Xingchun Zhao2, Xiaonuo Fan3, Ruoxia Yao1

The analysis of mixed short tandem repeat (STR) profiles has been long considered as a difficult challenge in the forensic DNA analysis. In the context of China, the current approach to analyze mixed STR profiles depends mostly on forensic manual method. However, besides the inefficiency, this technique is also susceptible to subjective biases in interpreting analysis results, which can hardly meet up with the growing demand for STR profiles analysis. In response, this study introduces an innovative method known as the global minimum residual method, which not only predicts the proportion of each contributor within a mixture, but also delivers accurate analysis results. The global minimum residual method first gives new definitions to the mixture proportion, then optimizes the allele model. After that, it comprehensively considers all loci present in the STR profile, accumulates and sums the residual values of each locus and selects the mixture proportion with the minimum accumulative sum as the inference result. Furthermore, the grey wolf optimizer is also employed to expedite the search for the optimal value. Notably, for two-person STR profiles, the high accuracy and remarkable efficiency of the global minimum residual method can bring convenience to realize extensive STR profile analysis. The optimization scheme established in this research has exhibited exceptional outcomes in practical applications, boasting significant utility and offering an innovative avenue in the realm of mixed STR profile analysis.

forensic genetics; mixed STR profile; allele model; grey wolf optimizer

2023-06-15;

2023-08-14;

2023-08-30

國家自然科學(xué)基金項(xiàng)目(編號(hào)：12271324)，陜西省自然科學(xué)基金重點(diǎn)項(xiàng)目(編號(hào)：2022ZJ-39)和法醫(yī)遺傳學(xué)公安部重點(diǎn)實(shí)驗(yàn)室開放課題(編號(hào)：2021FGKFKT07) 資助[Supported by the National Natural Science Foundation of China (No. 12271324), the Key Program of Shaanxi Natural Science Foundation (No. 2022ZJ-39), and Open Projects of the Key Laboratory of Forensic Genetics of the Ministry of Public Security (No. 2021FGKFKT07)]

李鑫，碩士研究生，專業(yè)方向：生物信息。E-mail: 15514872144@163.com

姚若俠，博士，教授，研究方向：符號(hào)計(jì)算。E-mail: rxyao@snnu.edu.cn

范虹，博士，教授，研究方向：生物信息。E-mail: fanhong@snnu.edu.cn

10.16288/j.yczz.23-101

(責(zé)任編委: 朱波峰)