MSAP熒光檢測初始數據的自動分析研究

陳秋博，鄒智元，李金龍，張競濤，徐青

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

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

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

1 材料和方法

1.1 材料

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

1.2 方法

MSAP 方法包括酶切反應、連接反應、擴增反應和熒光檢測共4 個主要步驟。本實驗中，酶切反應、連接反應和擴增反應條件與徐青等[6]所用條件相同。

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

2 結果

2.1 MSAP 擴增片段熒光信號的轉換和初始數據的獲得

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

2.2 MSAP熒光標記檢測初始數據的轉換

2.2.1 MSAP 熒光標記檢測初始數據的人工轉換在MSAP 圖譜上，每個樣本基因組對應2 條泳道，其中H 泳道是用HpaⅡ和EcoRⅠ處理的組織DNA 樣品，M 泳道是用MspⅠ和EcoRⅠ處理的組織DNA 樣品。根據擴增產物在2 條泳道出現的情況，個體基因組甲基化帶型可分為3 種：TypeⅠ為非甲基化模式（條帶在H和M 泳道同時出現），TypeⅡ為全甲基化模式（條帶在M 泳道出現而在H 泳道缺失），TypeⅢ為半甲基化模式（條帶在H 泳道出現而在M 泳道缺失）。個體基因組的甲基化水平為TypeⅡ+TypeⅢ與TypeⅠ+TypeⅡ+TypeⅢ的比值。

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

2.2.2 MSAP 熒光標記檢測初始數據的自動轉換由于人工轉換是對樣本毛細管電泳的熒光吸收峰圖譜進行一一比對來確定初始數據Excel 表格中的Size 值的有無，所以是最準確的轉換方法。當樣品數量較少時，可以采用這種方法。但是，在實際研究中，為了獲得較為準確的結果，樣本量需要達到一定規模。大的樣本量獲得的初始數據量非常龐大，人工轉換需要大量時間，很難完成。以人工轉換方法獲得的甲基化水平數據為標準，我們設計了3 種可通過計算機編程實現初始數據自動轉換的方法，用于替代人工轉換方法。這3 種方法可通過計算機編程實現初始數據自動轉換，確定實驗樣品基因組對應位點的甲基化帶型，計算實驗樣品基因組的甲基化水平。

2.2.2.1 直接取整法 將所有Size 值四舍五入為整數值，判斷每一位點H和M 泳道Size 值的有無，確定該位點的甲基化狀態。直接取整法的數據轉換流程為：程序讀取初始數據→實現所有位點Size 值的取整→設置每條泳道中50～800長度為1的區間依次遞增標識→將取整后Size 值填入對應標記的位置中→計算甲基化水平。如圖2DI所示，M 列中172這一值應當填入區間171，但是由于取整丟失了數據精確度，所以存在一定缺陷。

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

2.2.2.3 半區間放置法 與整區間放置法的數據轉換流程相同，但對泳道區間位置標識進行了調整。每條泳道中的位置從49.50 到800.50 用長度為1.00的區間依次標識。如圖2HR 所示，M 列中171.88 應當填入區間171.50～172.50，但171.88 與171.48 的差值小于0.5，經過智能校正，該值被填入170.50～171.50區間。

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

2.3 4種數據轉換方法計算結果比較

圖1 初始數據人工轉換流程圖

圖2 4種方法泳道位置標識及峰位置值填入

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

2.4 用于MSAP熒光檢測初始數據的自動分析的程序編寫

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

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

表1 4種數據轉換方法處理的結果

圖3 MSAP Analyst數據處理軟件的操作界面

3 討論

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

[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]徐青,張沅,孫東曉，等.應用MSAP 方法檢測雞不同組織基因組的甲基化狀態[J].遺傳,2011,33(6):620-626.