一種基于計(jì)算機(jī)的單核苷酸多態(tài)性位點(diǎn)檢測(cè)方法

馬 磊，賈奇男，張 俊，易青青，賀建峰，張 琪

(昆明理工大學(xué)信息工程與自動(dòng)化學(xué)院生物醫(yī)學(xué)工程研究所，云南昆明650500)

單核苷酸多態(tài)性(single nucleotide polymorphism，SNP)是指在染色體基因組水平上單個(gè)核苷酸的變異引起DNA序列多態(tài)性，它包括單堿基的轉(zhuǎn)換，顛換、插入及缺失等形式，SNP位點(diǎn)的檢出可以在一定程度上預(yù)測(cè)個(gè)體在自然選擇過(guò)程中的演變程度。目前，SNP位點(diǎn)檢測(cè)方法主要有單鏈構(gòu)象多態(tài)性分析(single strand conformation polymorphism，SSCP)［1］，變性梯度凝膠電泳(denatured gradient gel electrophoresis，DGGE)［2］，及其他一些常見方法［3－6］。但是這些方法均具有耗時(shí)長(zhǎng)，過(guò)程繁瑣和技術(shù)難度大，費(fèi)用高等缺點(diǎn)，從而制約了SNP的研究和發(fā)展。近年來(lái)，出現(xiàn)了一些利用生物信息技術(shù)等簡(jiǎn)單易行的方法，進(jìn)行SNP位點(diǎn)篩選分析的研究［7］，但仍存在應(yīng)用范圍較小等問(wèn)題。

本研究提出了一種基于計(jì)算機(jī)的SNP位點(diǎn)檢測(cè)方法，并以HBV的演化為例，應(yīng)用所提出的方法，通過(guò)特征信息的提取，剔除冗余的非疾病基因或非疾病風(fēng)險(xiǎn)基因，采用基于最優(yōu)風(fēng)險(xiǎn)與預(yù)防模式的數(shù)學(xué)算法，研究個(gè)體被HBV感染的風(fēng)險(xiǎn)性及被感染后疾病的演變程度，結(jié)論可為臨床診斷提供參考。

1 方法

1.1 基于信息熵的特征屬性值項(xiàng)的提取

在信息增益的DNA序列中，特征信息選擇的衡量標(biāo)準(zhǔn)是判斷該位置上的特征基因?yàn)槟繕?biāo)屬性帶來(lái)的信息量，信息越多，該特征基因就越重要。對(duì)一個(gè)特征基因而言，特定位置上它存在和不存在時(shí)信息量將發(fā)生明顯變化，而前后信息量的差值就是這個(gè)特征基因給該位置帶來(lái)的信息量，也稱為熵［8］。

假如有變量X，其可能的取值有n種，取到每一種的概率為pi，那么X的熵值計(jì)算公式1為:

X可能的變化越多，X所攜帶的信息量越大，熵也就越大。所以特征信息T給聚類C或者分類C帶來(lái)的信息增益為公式2:

這里包含兩種情況:一種是特征T出現(xiàn)，標(biāo)記為t，一種是特征T不出現(xiàn)，標(biāo)記為t'。所以計(jì)算見公式3所示。

由公式2可知特征與類別的信息增益公式。

該方法可在WEKA工作臺(tái)上執(zhí)行和實(shí)現(xiàn)［9］。

由于基因信息數(shù)據(jù)量大以及生物學(xué)數(shù)據(jù)本身的復(fù)雜性，導(dǎo)致了對(duì)這類數(shù)據(jù)進(jìn)行信息挖掘非常困難。因此在進(jìn)行特征基因選擇之前，首先將數(shù)據(jù)源進(jìn)行格式處理。數(shù)據(jù)源映射圖如圖1所示。sequences是樣本號(hào)對(duì)應(yīng)的每條序列，attribute1至attribute3221表示每條序列對(duì)應(yīng)位置上的堿基，class表示每條序列對(duì)應(yīng)的類屬性。

1.2 最優(yōu)風(fēng)險(xiǎn)與預(yù)防模式

最優(yōu)風(fēng)險(xiǎn)與預(yù)防模式算法MORE［10］的核心思想是采用局部支持度(supp)［11］來(lái)挖掘頻繁模式，其次應(yīng)用流行病理學(xué)中常用的相對(duì)風(fēng)險(xiǎn)值(RR)［12］度量指標(biāo)產(chǎn)生最優(yōu)風(fēng)險(xiǎn)與預(yù)防模式。這里簡(jiǎn)要介紹局部支持度(lsupp)、相對(duì)風(fēng)險(xiǎn)值與比值比的概念及計(jì)算方法。

由于在實(shí)際的醫(yī)療數(shù)據(jù)集中，數(shù)據(jù)量很大且正反類事例嚴(yán)重不平衡，因此采用局部支持度作為風(fēng)險(xiǎn)模式的支持度，即濾除非患病序列樣本。可以表示為樣本中同時(shí)出現(xiàn)模式P和a的概率與樣本中只出現(xiàn)a的概率的比值。其計(jì)算如公式4所示:

式中supp(p→a)表示模式P的支持度，即同時(shí)出現(xiàn)模式P和a的概率。局部支持度滿足反單調(diào)性，它的含義為:一個(gè)超集的支持度小于或等于它的任一子集的支持度。最優(yōu)風(fēng)險(xiǎn)與預(yù)防模式滿足反單調(diào)性原則就是最優(yōu)風(fēng)險(xiǎn)與預(yù)防模式能夠被挖掘的原因。在本研究中，如果一個(gè)模式的局部支持度大于給定的閾值，則這個(gè)模式就是頻繁的。

相對(duì)風(fēng)險(xiǎn)值表示某種疾病在特定人群中發(fā)生率是在非特定人群中發(fā)生率的多少倍。RR存在閾值的選取，例如閾值為1，則RR=1，表明此模式既不是風(fēng)險(xiǎn)模式也不是預(yù)防模式，此模式中的屬性對(duì)患者患病沒有影響;RR大于1，則說(shuō)明此模式為風(fēng)險(xiǎn)模式，此模式中的屬性為風(fēng)險(xiǎn)因子;RR小于1，表示此模式為預(yù)防模式，此模式中的屬性為預(yù)防因子。相對(duì)風(fēng)險(xiǎn)值(relative risk，RR)的計(jì)算如公式5所示。

表1 模式產(chǎn)生的可能性與輸出結(jié)果Table 1 The results and possibility of patterns

在表1中，－p表示模式P沒有發(fā)生的所有記錄。supp(－p)不在P模式下的所有記錄，－pa是包括a但不在P模式下的所有記錄。supp(p，a)表示supp(p→a)。

1.3 基于最優(yōu)風(fēng)險(xiǎn)和預(yù)防模式集的特征屬性值項(xiàng)差異分析

由于通過(guò)最優(yōu)風(fēng)險(xiǎn)與預(yù)防模式挖掘出的疾病相關(guān)模式或無(wú)關(guān)模式有時(shí)相互矛盾。例如，如果某些模式中的特征屬性值項(xiàng)同時(shí)出現(xiàn)在最優(yōu)風(fēng)險(xiǎn)和最優(yōu)預(yù)防模式中，有可能會(huì)給臨床診斷帶來(lái)偏差。針對(duì)該問(wèn)題，假設(shè)最優(yōu)風(fēng)險(xiǎn)和預(yù)防模式中每個(gè)特征屬性值項(xiàng)是相互獨(dú)立的，然后利用最優(yōu)風(fēng)險(xiǎn)模式中的特征屬性值項(xiàng)形成風(fēng)險(xiǎn)集，利用最優(yōu)預(yù)防模式的特征屬性值項(xiàng)產(chǎn)生預(yù)防集，并統(tǒng)計(jì)風(fēng)險(xiǎn)集和預(yù)防集中的特征屬性值項(xiàng)的頻率，去除相等頻率值的公共特征屬性值項(xiàng)。同時(shí)分別生成最優(yōu)風(fēng)險(xiǎn)頻率項(xiàng)集和最優(yōu)預(yù)防頻率項(xiàng)集。最優(yōu)風(fēng)險(xiǎn)頻率項(xiàng)集中的特征屬性值項(xiàng)為風(fēng)險(xiǎn)因子，表示患者攜帶該基因并有患病風(fēng)險(xiǎn)。最優(yōu)預(yù)防頻率項(xiàng)集中的特征屬性值項(xiàng)為預(yù)防因子，表示患者攜帶該基因不具有患病風(fēng)險(xiǎn)。最后對(duì)每一個(gè)特征屬性值項(xiàng)的頻率歸一化處理，得該特征屬性值項(xiàng)的權(quán)重值，總權(quán)重為100。

2 結(jié)果

2.1 實(shí)驗(yàn)數(shù)據(jù)源

本研究使用的數(shù)據(jù)是從 Genbank［13］下載。樣本集包括FJ349205-FJ349241，EU859930-EU859937和EU859939-EU859956共40條序列，包括9例急性乙型肝炎患者，22例慢性乙型肝炎患者，9例肝硬化患者。在乙型肝炎病毒序列中，實(shí)驗(yàn)的目的是挖掘HBV的SNP位點(diǎn)即單突變位點(diǎn)，因此這里將所有HBV序列的每一個(gè)垂直列映射為特征屬性進(jìn)行數(shù)據(jù)處理。HBV序列數(shù)據(jù)類型分為3類，即急性(acute)、慢性(chronic carrier)和肝硬化(cirrhosis)。

2.2 最優(yōu)風(fēng)險(xiǎn)和預(yù)防模式的生成

如前所述，將其應(yīng)用到乙型肝炎轉(zhuǎn)化關(guān)系研究中。由于實(shí)驗(yàn)數(shù)據(jù)源HBV序列經(jīng)過(guò)多條序列比對(duì)后長(zhǎng)為3 221 bp，則映射為3 221個(gè)垂直列。為了盡可能獲得更多的最優(yōu)風(fēng)險(xiǎn)與預(yù)防模式，對(duì)模式長(zhǎng)度和相對(duì)風(fēng)險(xiǎn)閾值進(jìn)行了多次選取實(shí)驗(yàn)，最后選定了一個(gè)最佳方案即在乙型肝炎轉(zhuǎn)化關(guān)系研究中，設(shè)置模式長(zhǎng)度為6，特征屬性的閾值為0.20，相對(duì)風(fēng)險(xiǎn)閾值為1。在此條件下，實(shí)驗(yàn)共返回400個(gè)最優(yōu)風(fēng)險(xiǎn)與預(yù)防模式，分別為36個(gè)最優(yōu)風(fēng)險(xiǎn)模式和364個(gè)最優(yōu)預(yù)防模式。限于篇幅，這里只列舉了部分具有代表性的最優(yōu)風(fēng)險(xiǎn)模式(表2)和最優(yōu)預(yù)防模式(表3)。

表2 在局部支持度為0.05，模式長(zhǎng)度為6，特征屬性選取閾值為0.20，相對(duì)風(fēng)險(xiǎn)閾值為1的條件下由HBV序列生成的部分最優(yōu)風(fēng)險(xiǎn)模式Table 2 lsupp=0.05，pattern length=6，threshold of feature property=0.20，RR threshold=1，selected optimal risk patterns are generated from HBV sequence

針對(duì)表2和表3的部分實(shí)驗(yàn)結(jié)果，以最優(yōu)風(fēng)險(xiǎn)模式中的Pattern 1為例解釋說(shuō)明。模式中l(wèi)ength=1表示模式長(zhǎng)度為1，說(shuō)明此模式包括1個(gè)特征屬性值項(xiàng)，RR=3.6667表示相對(duì)風(fēng)險(xiǎn)值為3.6667。

2.3 乙肝病毒序列的差異分析

在此實(shí)驗(yàn)中，假設(shè)模式中每一個(gè)特征屬性值項(xiàng)是相互獨(dú)立的，HBV序列的特征屬性值項(xiàng)的最優(yōu)風(fēng)險(xiǎn)與預(yù)防權(quán)重如表4和表5所示，表示為最優(yōu)風(fēng)險(xiǎn)集和最優(yōu)預(yù)防集中各特征屬性值項(xiàng)的權(quán)重。每個(gè)特征屬性值項(xiàng)的權(quán)重來(lái)自它們?cè)谧顑?yōu)風(fēng)險(xiǎn)與預(yù)防集中的百分比。它可以了解某個(gè)特征屬性值項(xiàng)對(duì)患者患某種疾病的風(fēng)險(xiǎn)性與預(yù)防性。

表4中，特征屬性值項(xiàng)attribute2531=T出現(xiàn)在最優(yōu)風(fēng)險(xiǎn)頻率集中，attribute2531對(duì)應(yīng)在HBV序列中的位置為第2531位堿基位點(diǎn)。風(fēng)險(xiǎn)權(quán)重為10.4427，是最優(yōu)風(fēng)險(xiǎn)集中最大的風(fēng)險(xiǎn)權(quán)重，表明attribute2531在HBV序列第2531位堿基為T時(shí)發(fā)生急性乙肝轉(zhuǎn)慢性乙肝的可能性在所有特征屬性值項(xiàng)中最大，這些是導(dǎo)致此處發(fā)生堿基突變的決定因素。attribute2567=T出現(xiàn)在最優(yōu)預(yù)防集中，且預(yù)防權(quán)重為8.9241，表明attribute2567在 HBV序列第2567位堿基為T時(shí)不發(fā)生急性乙肝轉(zhuǎn)慢性乙肝的可能性很大，是此處防止堿基突變的決定因素。如果某一屬性同時(shí)出現(xiàn)在最優(yōu)風(fēng)險(xiǎn)集與最優(yōu)預(yù)防集中，說(shuō)明此處為決定急性轉(zhuǎn)慢性的候選SNPs位點(diǎn)。

根據(jù)上述表述，基于表4中最優(yōu)風(fēng)險(xiǎn)權(quán)重集，實(shí)驗(yàn)共檢測(cè)出18處急性轉(zhuǎn)慢性候選SNPs位點(diǎn)，其中4 處屬于堿基替換突變(nt1655，nt1753，nt1762，nt1764)，1處屬于堿基缺失突變(nt1896)，這5處均已在前期的文獻(xiàn)中報(bào)道過(guò)［14－15］。其余13處是新發(fā)現(xiàn)的候選SNP位點(diǎn)(其中nt2576為替換突變，其他12處為缺失突變)?；诒?中最優(yōu)風(fēng)險(xiǎn)權(quán)重集，在實(shí)驗(yàn)中共檢測(cè)出3處慢性乙型肝炎轉(zhuǎn)肝硬化候選SNPs位點(diǎn)，其中2處屬于堿基替換突變(nt2159、nt902)，1處屬于堿基缺失突變(nt507)。

3 討論

近年來(lái)SNP的檢測(cè)方法已被廣泛研究，國(guó)內(nèi)外專家學(xué)者也相應(yīng)提出了多種方法檢測(cè)SNP，但是大多需要依賴昂貴的儀器或?qū)I(yè)人員的技術(shù)支持，且操作相對(duì)復(fù)雜。本研究針對(duì)40條HBV病毒序列數(shù)據(jù)(9條急性乙型肝炎序列，22條慢性乙型肝炎序列，9條肝硬化序列)提出了一種基于最優(yōu)風(fēng)險(xiǎn)與預(yù)防模式算法來(lái)研究HBV病毒序列的SNP位點(diǎn)檢測(cè)問(wèn)題。實(shí)驗(yàn)共檢測(cè)出18處急性乙型肝炎轉(zhuǎn)慢性乙型肝炎候選SNPs位點(diǎn)，其中5處已有報(bào)道，另外13處是新發(fā)現(xiàn)的候選SNP位點(diǎn)。同時(shí)檢測(cè)出3處慢性乙型肝炎轉(zhuǎn)肝硬化候選SNPs位點(diǎn)。該方法成本低廉，操作簡(jiǎn)便，并能在龐大的基因數(shù)據(jù)中選出SNP位點(diǎn)，可以對(duì)乙型肝炎的臨床診斷和生物醫(yī)學(xué)研究起到有益的參考和借鑒作用，有可能成為適用于臨床針對(duì)肝炎及其他疾病的SNPs檢測(cè)方法。

optimal risk and preventive pattern(supp=0.43，info gain＞0.15，pattern length=7)property in risk pattern property in preventi ve pattern risk factor loci weight preventive factor loci weight attribute2531=T nt2531 10.4427 attribute2567=T nt2567 8.9241 attribute1896=G nt1896 9.2698 attribute1773=T nt1773 8.0144 attribute1588=A nt1588 7.3401 attribute3090=A nt3090 7.0107 attribute2689=T nt2689 6.8104 attribute1421=G nt1421 5.756 attribute1320=C nt1320 6.5834 attribute1657=C nt1657 5.458 attribute1655=T nt1655 6.3942 attribute49=G nt49 5.3795 attribute1762=T nt1762 6.3564 attribute1764=G nt1764 5.0502 attribute681=C nt681 5.1457 attribute2576=T nt2576 4.5483 attribute3073=T nt3073 4.7673 attribute1999=A nt1999 4.3444 attribute3094=C nt3094 4.5781 attribute1334=C nt1334 4.313 attribute1764=A nt1764 4.0484 attribute917=G nt917 4.313 attribute1753=T nt1753 4.0484 attribute3129=A nt3129 4.2033 attribute3094=C nt3094 3.3674 attribute1762=T nt1762 4.1248 attribute2819=A nt2819 3.1025 attribute600=C nt600 4.1248 attribute2687=G nt2687 2.9512 attribute2260=A nt2260 3.9837 attribute2576=G nt2576 2.6485 attribute8=A nt8 3.3407 attribute2588=C nt2588 2.4593 attribute2717=T nt2717 3.2465 attribute681=C nt681 2.1188 attribute9=A nt9 3.2152 attribute1655=C nt1655 3.2152 attribute513=C nt513 3.074 attribute1753=A nt1753 2.3683 attribute896=T nt896 1.9918

表5 慢性乙型肝炎與肝硬化數(shù)據(jù)生成的最優(yōu)風(fēng)險(xiǎn)集與預(yù)防集中的每個(gè)特征屬性的頻率并按照降序排列，每一個(gè)位置對(duì)應(yīng)的是HBV序列中的堿基位置Table 5 The descending ranking of frequency of each feature p?roperty in the optimal risk and preventive sets generated from the data of chronic HBV carrier and cirrhosis

［1］Kubo KS，Stuart RM，F(xiàn)reitas-Astúa J，et al．Evaluation of the genetic variability of orchid fleck virus by single-strand conformational polymorphism analysis and nucleotide sequencing of a fragment from the nucleocapsid gene［J］．Arch Virol，2009，154:1009 －1014．

［2］Drabovich AP，Krylov SN．Identification of base pairs in single-nucleotide polymorphisms by MutS protein-mediated capillary electrophoresis［J］．Anal Chem，2006，78:2035－2038．

［3］Jiang HH，Huang SY，Zhou DH，et al．Genetic characterization of Toxoplasma gondii from pigs from different localities in China by PCR-RFLP［J］．Parasites ＆ Vectors，2013，6:227．

［4］Dunnen JT，Antonarakis SE．Mutation nomenclature extensions and suggestions to describe complex mutations:A discussion［J］．Human Mutat，2000，15:7 －12．

［5］ Sapolskya RJ，Hsie L，Berno A，et al．High-throughput polymorphism screening and genotyping with high-density oligonucleotide arrays［J］．Elsevier，1999，14:187 － 192．

［6］ Gross E，Arnold N，Goette J，et al．A comparison of BRCAI mutation analysis by direct sequencing，SSCP and DHPLC［J］．Human Genet，1999，105:72 －78．

［7］陳客宏，曾靈，顧瑋等．利用生物信息學(xué)技術(shù)挑選TLR2基因標(biāo)簽單核苷酸多態(tài)性位點(diǎn)［J］．基礎(chǔ)醫(yī)學(xué)與臨床，2010，30:242－245．

［8］ Cover TM，Thomas JA．Elements of information theory［M］．Wiley-interscience，2006．

［9］Hall M，F(xiàn)rank E，Holmes G，et al．The WEKA data mining software:an update［J］．SIGKDD Explor Newsl，2009，11:10－18．

［10］Li JY，F(xiàn)u AW，He HX，et al．Efficient discovery of risk patterns in medical data［J］．Artif Intell Med，2009，45:77－89．

［11］Ohsaki M，Kitaguchi S，Okamoto K，et al．Evaluation of rule interestingness measures with a clinical dataset on hepatitis．In:Boulicaut J-FF．，Esposito F，Giannotti F，Pedreschi D，editors．Proceedings of the eighth European conference on principles and practice of knowledge discovery in databases［M］，New York:Springer-Verlag，2004:362－373．

［12］Triola MM，Triola MF．Biostatistics for the biological and health sciences［M］，Boston:Addison-Wesley，2006．

［13］Benson DA，Cavanaugh M，Clark K，et al．GenBank［J］．Nucleic Acids Res，2013，41:36 －42．

［14］Yuan JM，Ambinder A，F(xiàn)an Y，et al．Prospective evaluation of hepatitis B 1762T/1764A mutations on hepatocellular carcinoma development in Shanghai，China［J］．Cancer Epidemiol Biomarkers＆ Prevention，2009，18:590－594．

［15］Wang Y，Liu H，Zhou Q，et al．Analysis of point mutation in site 1896 of HBV precore and its detection in the tissues and serum of HCC patients［J］．World JGastroenterol，2000，6:395－397．