陳秋博，鄒智元，李金龍，張競(jìng)濤，徐青

北京交通大學(xué) a.生命科學(xué)與生物工程研究院；b.軟件學(xué)院；北京 100044

DNA 甲基化是真核細(xì)胞基因組重要的表觀遺傳學(xué)修飾方式之一，能在DNA 序列不發(fā)生變化的前提下調(diào)控基因的表達(dá)[1]。DNA 甲基化修飾主要發(fā)生在基因組CCGG 位點(diǎn)[2-3]，與細(xì)胞分化、染色體失活及基因印記等生物學(xué)過程密切相關(guān)[4]。目前，檢測(cè)DNA 甲基化的方法主要有酶切法、亞硫酸鹽測(cè)序法、甲基化特異性PCR、基因組限制性酶切掃描法等。其中,以酶切和PCR 為基礎(chǔ)的甲基化敏感擴(kuò)增片段多態(tài)性（methylation sensitive amplified polymorphism，MSAP）不但敏感性強(qiáng)，而且只需要常規(guī)的儀器，操作比較簡單，適用于沒有任何信息的全基因組水平的甲基化分析。MSAP 方法由擴(kuò)增片段長度多態(tài)性（AFLP）方法改進(jìn)而來[5]，隨后被廣泛應(yīng)用于動(dòng)植物基因組甲基化相關(guān)研究中。過去，一般采用聚丙烯酰胺凝膠電泳（PAGE）銀染方法檢測(cè)MSAP的選擇性擴(kuò)增產(chǎn)物，但因PAGE 法存在操作復(fù)雜、檢測(cè)通量低和分析時(shí)間長等缺點(diǎn)，逐漸被以熒光標(biāo)記為基礎(chǔ)的測(cè)序膠和毛細(xì)管電泳技術(shù)所替代。

隨著自動(dòng)化程度和靈敏度的增強(qiáng)，毛細(xì)管電泳技術(shù)可以更靈敏、更有效、更高通量地獲得MSAP 的檢測(cè)數(shù)據(jù)。例如，本課題組使用以毛細(xì)管電泳技術(shù)為基礎(chǔ)的ABI 3730xl 遺傳分析儀，平均90 min 內(nèi)可以獲得96 泳道的熒光初始檢測(cè)數(shù)據(jù)，與傳統(tǒng)方法相比，檢測(cè)時(shí)間縮短至約1/12，但產(chǎn)生的數(shù)據(jù)量龐大，約為PAGE-銀染檢測(cè)的8～10 倍。如果每條泳道平均產(chǎn)生200 條數(shù)據(jù)，那么96 泳道可產(chǎn)生近20 000 條數(shù)據(jù)，而往往一個(gè)普通的實(shí)驗(yàn)設(shè)計(jì)需要檢測(cè)約10 個(gè)96 泳道。如此巨大的數(shù)據(jù)量如果依賴傳統(tǒng)的人工手動(dòng)分析，幾乎不可能完成。因此，為了建立能夠?qū)崿F(xiàn)MSAP 熒光檢測(cè)初始數(shù)據(jù)自動(dòng)分析的方法，我們采用毛細(xì)管電泳熒光MSAP 檢測(cè)技術(shù)，分析了北京油雞基因組的甲基化程度，對(duì)產(chǎn)生的MSAP 初始數(shù)據(jù)采用標(biāo)準(zhǔn)的人工轉(zhuǎn)換方法及3 種不同的計(jì)算機(jī)自動(dòng)轉(zhuǎn)換方法進(jìn)行處理，比較和分析了人工轉(zhuǎn)換方法和計(jì)算機(jī)自動(dòng)轉(zhuǎn)換方法獲得的結(jié)果，確定了能夠替代人工轉(zhuǎn)換方法的自動(dòng)轉(zhuǎn)換方法，從而實(shí)現(xiàn)了MSAP 熒光檢測(cè)初始數(shù)據(jù)的自動(dòng)分析，并進(jìn)行編程，創(chuàng)建了適用于MSAP 熒光檢測(cè)初始數(shù)據(jù)自動(dòng)分析的在線軟件。

1 材料和方法

1.1 材料

選擇10只17周齡北京油雞作為實(shí)驗(yàn)材料，提取基因組DNA。接頭與引物序列設(shè)計(jì)參照徐青等[6]方法，由上海生工生物工程技術(shù)服務(wù)有限公司合成。

1.2 方法

MSAP 方法包括酶切反應(yīng)、連接反應(yīng)、擴(kuò)增反應(yīng)和熒光檢測(cè)共4 個(gè)主要步驟。本實(shí)驗(yàn)中，酶切反應(yīng)、連接反應(yīng)和擴(kuò)增反應(yīng)條件與徐青等[6]所用條件相同。

毛細(xì)管電泳熒光檢測(cè)與數(shù)據(jù)分析方法如下：在96 孔板的各孔中分別加入3μL 選擇性擴(kuò)增產(chǎn)物、6.67μL甲酰胺、0.33μL ROX800分子量內(nèi)標(biāo)，95℃變性5 min，冰上放置10 min，離心1 min，在ABI 3730xl DNA 分析儀上進(jìn)行毛細(xì)管電泳（10 kV 預(yù)電泳1 min，3 kV 進(jìn)樣15 s，10 kV 電泳70 min），收集初始數(shù)據(jù)，系統(tǒng)將各峰值的位置與其泳道中的ROX800 分子量內(nèi)標(biāo)做比較，通過GeneMapper 4.0軟件對(duì)收集的初始數(shù)據(jù)進(jìn)行分析。

2 結(jié)果

2.1 MSAP 擴(kuò)增片段熒光信號(hào)的轉(zhuǎn)換和初始數(shù)據(jù)的獲得

應(yīng)用GeneMapper 4.0 軟件將MSAP 擴(kuò)增片段熒光檢測(cè)信號(hào)轉(zhuǎn)換為MSAP 擴(kuò)增片段的實(shí)體數(shù)據(jù)。首先對(duì)所有泳道的分子量內(nèi)標(biāo)的熒光圖譜進(jìn)行校正，依據(jù)每個(gè)泳道校正后的內(nèi)標(biāo)獲得相應(yīng)泳道的所有MSAP 擴(kuò)增DNA 片段的大小，將結(jié)果導(dǎo)至Excel 表格中。Excel 表格的初始數(shù)據(jù)包括峰位置（Size）、峰高值（Height）、峰面積（Area）和數(shù)據(jù)點(diǎn)（Data Point）等4項(xiàng)有效數(shù)據(jù)。其中，Height是MSAP擴(kuò)增片段的熒光信號(hào)強(qiáng)度，代表擴(kuò)增片段的拷貝數(shù)，在本研究中，Height 的閾值設(shè)為50，熒光強(qiáng)度小于50 的片段不予考慮。Size表示檢測(cè)樣品MSAP擴(kuò)增片段的長度，基本上為非整數(shù)值。在實(shí)驗(yàn)樣品的MSAP 熒光圖譜上，擴(kuò)增片段集中在50～800 bp，所以隨后的數(shù)據(jù)分析只統(tǒng)計(jì)長度為50～800 bp 的片段。在本研究中，定義Size≥50 的值為擴(kuò)增產(chǎn)物在該位點(diǎn)出現(xiàn)，而Size=0 或不出現(xiàn)為擴(kuò)增產(chǎn)物在該位點(diǎn)缺失。這樣，通過定義Height和Size 值的區(qū)間，對(duì)Excel 表格的初始數(shù)據(jù)進(jìn)行過濾，獲得用于進(jìn)行實(shí)驗(yàn)樣品DNA 甲基化分析的初始數(shù)據(jù)，如圖1A所示。

2.2 MSAP熒光標(biāo)記檢測(cè)初始數(shù)據(jù)的轉(zhuǎn)換

2.2.1 MSAP 熒光標(biāo)記檢測(cè)初始數(shù)據(jù)的人工轉(zhuǎn)換在MSAP 圖譜上，每個(gè)樣本基因組對(duì)應(yīng)2 條泳道，其中H 泳道是用HpaⅡ和EcoRⅠ處理的組織DNA 樣品，M 泳道是用MspⅠ和EcoRⅠ處理的組織DNA 樣品。根據(jù)擴(kuò)增產(chǎn)物在2 條泳道出現(xiàn)的情況，個(gè)體基因組甲基化帶型可分為3 種：TypeⅠ為非甲基化模式（條帶在H和M 泳道同時(shí)出現(xiàn)），TypeⅡ?yàn)槿谆Ｊ剑l帶在M 泳道出現(xiàn)而在H 泳道缺失），TypeⅢ為半甲基化模式（條帶在H 泳道出現(xiàn)而在M 泳道缺失）。個(gè)體基因組的甲基化水平為TypeⅡ+TypeⅢ與TypeⅠ+TypeⅡ+TypeⅢ的比值。

在MSAP 圖譜中，H和M 泳道的擴(kuò)增片段對(duì)應(yīng)于MSAP 熒光標(biāo)記檢測(cè)初始數(shù)據(jù)表中的H和M 數(shù)據(jù)列的Size 值。如圖1A 所示，初始數(shù)據(jù)中每個(gè)擴(kuò)增片段的Size 值幾乎全為非整數(shù)值，而實(shí)際擴(kuò)增片段應(yīng)為整數(shù)值，所以初始數(shù)據(jù)的非整數(shù)Size 值要轉(zhuǎn)換為對(duì)應(yīng)的整數(shù)Size 值。初始數(shù)據(jù)Size 值的轉(zhuǎn)換處理對(duì)于樣品基因組的甲基化水平估算具有非常大的影響。在本研究中，我們定義人工轉(zhuǎn)換為初始數(shù)據(jù)的標(biāo)準(zhǔn)處理。如圖1C 所示，初始數(shù)據(jù)的標(biāo)準(zhǔn)轉(zhuǎn)換是通過直接比對(duì)樣本毛細(xì)管電泳的熒光吸收峰圖譜和MSAP 熒光標(biāo)記檢測(cè)初始數(shù)據(jù)Excel 表格中的Size值，來確認(rèn)每個(gè)擴(kuò)增位點(diǎn)H和M 數(shù)據(jù)列的Size 值的有或無。MSAP 檢測(cè)中，每個(gè)樣品每對(duì)引物每條泳道平均獲得約200 條目的峰，完成1 個(gè)樣品10 對(duì)引物獲得的MSAP 熒光標(biāo)記檢測(cè)初始數(shù)據(jù)的直接轉(zhuǎn)換需要大約3 h。

2.2.2 MSAP 熒光標(biāo)記檢測(cè)初始數(shù)據(jù)的自動(dòng)轉(zhuǎn)換由于人工轉(zhuǎn)換是對(duì)樣本毛細(xì)管電泳的熒光吸收峰圖譜進(jìn)行一一比對(duì)來確定初始數(shù)據(jù)Excel 表格中的Size 值的有無，所以是最準(zhǔn)確的轉(zhuǎn)換方法。當(dāng)樣品數(shù)量較少時(shí)，可以采用這種方法。但是，在實(shí)際研究中，為了獲得較為準(zhǔn)確的結(jié)果，樣本量需要達(dá)到一定規(guī)模。大的樣本量獲得的初始數(shù)據(jù)量非常龐大，人工轉(zhuǎn)換需要大量時(shí)間，很難完成。以人工轉(zhuǎn)換方法獲得的甲基化水平數(shù)據(jù)為標(biāo)準(zhǔn)，我們?cè)O(shè)計(jì)了3 種可通過計(jì)算機(jī)編程實(shí)現(xiàn)初始數(shù)據(jù)自動(dòng)轉(zhuǎn)換的方法，用于替代人工轉(zhuǎn)換方法。這3 種方法可通過計(jì)算機(jī)編程實(shí)現(xiàn)初始數(shù)據(jù)自動(dòng)轉(zhuǎn)換，確定實(shí)驗(yàn)樣品基因組對(duì)應(yīng)位點(diǎn)的甲基化帶型，計(jì)算實(shí)驗(yàn)樣品基因組的甲基化水平。

2.2.2.1 直接取整法 將所有Size 值四舍五入為整數(shù)值，判斷每一位點(diǎn)H和M 泳道Size 值的有無，確定該位點(diǎn)的甲基化狀態(tài)。直接取整法的數(shù)據(jù)轉(zhuǎn)換流程為：程序讀取初始數(shù)據(jù)→實(shí)現(xiàn)所有位點(diǎn)Size 值的取整→設(shè)置每條泳道中50～800長度為1的區(qū)間依次遞增標(biāo)識(shí)→將取整后Size 值填入對(duì)應(yīng)標(biāo)記的位置中→計(jì)算甲基化水平。如圖2DI所示，M 列中172這一值應(yīng)當(dāng)填入?yún)^(qū)間171，但是由于取整丟失了數(shù)據(jù)精確度，所以存在一定缺陷。

2.2.2.2 整區(qū)間放置法 將直接取整法的泳道整數(shù)位置標(biāo)識(shí)調(diào)整為區(qū)間位置標(biāo)識(shí)，數(shù)據(jù)轉(zhuǎn)換流程為：程序讀取初始數(shù)據(jù)→設(shè)置每條泳道中50.00～800.99 長度為1 的區(qū)間依次遞增標(biāo)識(shí)→初始Size 值直接填入所屬區(qū)間標(biāo)識(shí)的位置中→程序智能校正填入的初始Size 值，處理重復(fù)的Size 值→計(jì)算甲基化水平。整區(qū)間放置法沒有對(duì)初始Size 值直接進(jìn)行四舍五入，所以保留了數(shù)據(jù)精度，從而可以使用程序?qū)ize 值進(jìn)行智能校正。如圖2WR 所示，M 列中167.91 應(yīng)當(dāng)填入?yún)^(qū)間167.00～167.99，但167.91 與168.07 的差值小于0.5，經(jīng)過智能校正，該值被填入?yún)^(qū)間168.00～168.99。168.9和169.89同樣是這種情況。

2.2.2.3 半?yún)^(qū)間放置法 與整區(qū)間放置法的數(shù)據(jù)轉(zhuǎn)換流程相同，但對(duì)泳道區(qū)間位置標(biāo)識(shí)進(jìn)行了調(diào)整。每條泳道中的位置從49.50 到800.50 用長度為1.00的區(qū)間依次標(biāo)識(shí)。如圖2HR 所示，M 列中171.88 應(yīng)當(dāng)填入?yún)^(qū)間171.50～172.50，但171.88 與171.48 的差值小于0.5，經(jīng)過智能校正，該值被填入170.50～171.50區(qū)間。

在后文中，用4 種處理方法的英文單詞首字母縮寫代表這4 種方法。AA 為人工分析法，DI為直接取整法，WR為整區(qū)間放置法，HR為半?yún)^(qū)間放置法。

2.3 4種數(shù)據(jù)轉(zhuǎn)換方法計(jì)算結(jié)果比較

圖1 初始數(shù)據(jù)人工轉(zhuǎn)換流程圖

圖2 4種方法泳道位置標(biāo)識(shí)及峰位置值填入

北京油雞基因組DNA甲基化程度MSAP熒光檢測(cè)初始數(shù)據(jù)通過人工轉(zhuǎn)換和3 種自動(dòng)轉(zhuǎn)換方法進(jìn)行分析，4 種數(shù)據(jù)轉(zhuǎn)換方法處理結(jié)果見表1。方差分析顯示，對(duì)于非甲基化位點(diǎn)和半甲基化位點(diǎn)的分析，4種轉(zhuǎn)換方法獲得的結(jié)果存在顯著差異（P＜0.05）；對(duì)于全甲基化位點(diǎn)，4 種轉(zhuǎn)換方法獲得的結(jié)果沒有顯著差異。Duncan 多重比較分析顯示，對(duì)于非甲基化位點(diǎn)和全甲基化位點(diǎn)的數(shù)據(jù)轉(zhuǎn)換，DI 與其他3 種方法之間差異顯著（P＜0.05），而AA、WR和HR 三種方法之間沒有顯著差異，HR 與AA 方法之間差異最小。綜合以上結(jié)果，對(duì)于MSAP 熒光檢測(cè)初始數(shù)據(jù)的轉(zhuǎn)換，與AA 標(biāo)準(zhǔn)轉(zhuǎn)換相比，HR 轉(zhuǎn)換結(jié)果與AA 最接近，WR 次之，DI 差異最大達(dá)顯著水平（P＜0.05），所以在處理較大量的數(shù)據(jù)時(shí)，可以采用HR和WR 來替代AA方法，不建議使用DI方法。

2.4 用于MSAP熒光檢測(cè)初始數(shù)據(jù)的自動(dòng)分析的程序編寫

MSAP 方法已被廣泛應(yīng)用于動(dòng)植物基因組的甲基化研究。隨著毛細(xì)管電泳技術(shù)在MSAP 上的應(yīng)用，越來越龐大的數(shù)據(jù)量使研究人員的分析遇到了難題，初始數(shù)據(jù)的計(jì)算機(jī)自動(dòng)化處理已成為高通量MSAP 分析的迫切需要。目前一些公司提供的熒光圖譜收集儀器和配套軟件已經(jīng)初步使初始數(shù)據(jù)的獲得自動(dòng)化，但是在獲得這些數(shù)據(jù)之后，研究人員應(yīng)結(jié)合實(shí)際意義進(jìn)行二次處理，目前還沒有此類成型軟件。為了解決這個(gè)問題，結(jié)合本研究分析結(jié)果，我們?cè)诎雲(yún)^(qū)間放置法的基礎(chǔ)上編寫了一個(gè)數(shù)據(jù)分析軟件——MSAP Analyst。

該軟件采用Microsoft.Net 4.0架構(gòu)，使用C#語言開發(fā)。軟件分為界面模塊、數(shù)據(jù)處理模塊和算法模塊，處理的數(shù)據(jù)文件格式為“.csv”或“.xls”。軟件的使用流程如圖3 所示，主要分為三步：①打開MSAP初始數(shù)據(jù)表格，啟動(dòng)軟件，點(diǎn)擊左上角菜單“文件-打開”，在彈出的窗口中選擇一個(gè)或多個(gè)初始數(shù)據(jù)表格并打開；②點(diǎn)擊左上角菜單“分析-設(shè)置導(dǎo)出選項(xiàng)”，如圖所示，在彈出的新窗口中可以勾選需要統(tǒng)計(jì)的泳道（與96 孔板上樣孔一一對(duì)應(yīng)），點(diǎn)擊“確定”完成；③選擇導(dǎo)出方式并導(dǎo)出，點(diǎn)擊“分析-生成統(tǒng)計(jì)表”可對(duì)一個(gè)原數(shù)據(jù)文件進(jìn)行統(tǒng)計(jì)，點(diǎn)擊“分析-生成所有”則統(tǒng)計(jì)當(dāng)前打開的所有原數(shù)據(jù)文件。在生成統(tǒng)計(jì)表時(shí)，還可以根據(jù)實(shí)驗(yàn)樣本實(shí)際放置的位置，即對(duì)應(yīng)關(guān)系，選擇“橫向”或“縱向”生成統(tǒng)計(jì)表。我們應(yīng)用MSAP Analyst 模擬進(jìn)行了多種類似數(shù)據(jù)的處理，都獲得了理想的分析結(jié)果，因此以半?yún)^(qū)間放置法為基礎(chǔ)的MSAP Analyst 軟件可以應(yīng)用于DNA 甲基化的MSAP 熒光檢測(cè)初始數(shù)據(jù)的自動(dòng)分析。另外，為了方便大家使用，我們?cè)O(shè)計(jì)了MSAP Analyst 軟件的在線分析網(wǎng)站（http://www.fly2leo.cn），支持軟件、測(cè)試數(shù)據(jù)和使用說明書的下載。

表1 4種數(shù)據(jù)轉(zhuǎn)換方法處理的結(jié)果

圖3 MSAP Analyst數(shù)據(jù)處理軟件的操作界面

3 討論

人工分析法是直接將2 種酶切反應(yīng)的熒光圖譜進(jìn)行重合比較，根據(jù)熒光圖譜的差異確定甲基化類型，計(jì)算甲基化程度。數(shù)據(jù)處理過程中，可以及時(shí)地進(jìn)行人工調(diào)整，所以能夠最大程度地減少誤差，得到的分析結(jié)果最準(zhǔn)確，因而把這個(gè)結(jié)果定義為標(biāo)準(zhǔn)結(jié)果，把人工分析法定義為標(biāo)準(zhǔn)方法。而對(duì)于3 種自動(dòng)轉(zhuǎn)換方法來說，由于堿基的最小值是1 bp，所以設(shè)置的區(qū)間跨度只能為1，因此3種分析方法都會(huì)出現(xiàn)2 個(gè)峰位置值在同一泳道位置（距離較短，近似重合）重復(fù)出現(xiàn)的情況，即重復(fù)位置值，重復(fù)位置值的出現(xiàn)會(huì)導(dǎo)致甲基化水平的計(jì)算誤差。直接取整法中，對(duì)峰位置值直接進(jìn)行簡單的四舍五入，忽略了初始數(shù)據(jù)的小數(shù)部分，產(chǎn)生了沒有進(jìn)行任何糾正的重復(fù)位置值，導(dǎo)致了較大的計(jì)算誤差。為了減小重復(fù)位置值的出現(xiàn)，引入?yún)^(qū)間放置的概念。整區(qū)間放置法通過程序智能校正峰位置值，減少了部分重復(fù)位置值，所以甲基化程度的計(jì)算誤差比直接取整法小。根據(jù)人工分析法獲得的經(jīng)驗(yàn)，對(duì)整區(qū)間放置法進(jìn)行了改進(jìn)，得到半?yún)^(qū)間放置法。半?yún)^(qū)間放置法本身與樣品峰值圖吻合度最高，而且經(jīng)過程序智能校正，重復(fù)位置值進(jìn)一步減少，雖然不能完全消除重復(fù)位置值的產(chǎn)生，但可以保證甲基化程度的計(jì)算誤差最小。為了最大程度地減少重復(fù)位置值的出現(xiàn)，我們?cè)O(shè)計(jì)的軟件中可以將每個(gè)泳道的重復(fù)值進(jìn)行標(biāo)識(shí)，使用人員可根據(jù)原始的峰位置圖進(jìn)行手動(dòng)糾正這些重復(fù)位置值，以進(jìn)一步減小自動(dòng)分析與標(biāo)準(zhǔn)分析之間的誤差。

[1]Jones P A.Functions of DNA methylation:islands,start sites,gene bodies and beyond[J].Nat Rev Genet,2012,13(7):484-492.

[2]Smith Z D,Meissner A.DNA methylation:roles in mammalian development[J].Nat Rev Genet,2013,14(3):204-220.

[3]Moore L D,Le T,Fan G.DNA methylation and its basic function[J].Neuropsychopharmacology,2013,38(7):23-38.

[4]Gupta R,Nagarajan A,Wajapeyee N.Advances in genomewide DNA methylation analysis[J].Biotechniques,2010,49(4):iii-xi.

[5]Xu M,Li X,Korban S S.AFLP-based detection of DNA methylation[J].Plant Mol Biol Rep,2000,18(4):361-368.

[6]徐青,張沅,孫東曉，等.應(yīng)用MSAP 方法檢測(cè)雞不同組織基因組的甲基化狀態(tài)[J].遺傳,2011,33(6):620-626.