原貼地址：https://blog.csdn.net/joshua_hit/article/details/54982018
?????????????? https://bioconductor.org/packages/release/bioc/vignettes/ChAMP/inst/doc/ChAMP.html
我目前的研究領(lǐng)域是DNA甲基化琅捏，這個領(lǐng)域并沒有很好的分析軟件柄延，頂多一些R包吧搜吧，而且很多學(xué)生物的人滤奈，想要分析DNA甲基化數(shù)據(jù)并不容易蜒程，雖然這樣一些包已經(jīng)存在的昭躺，但是首先领炫，他們不一定易用帝洪，例如minfi（我領(lǐng)域第一號分析包）碟狞，開發(fā)人在寫的時候族沃，就已經(jīng)默認使用者有相當(dāng)?shù)腞語言知識和統(tǒng)計知識泌参。其次是沽一，包的功能分散铣缠，并 不集中蝗蛙，經(jīng)常需要從一個軟件跑出來的東西捡硅，弄來弄去導(dǎo)入到其他軟件中去進行下一步分析壮韭。
我在英國留學(xué)的時候喷屋，遇到了ChAMP包的原作者屯曹，她已經(jīng)不想再管這個軟件了是牢。所以我接手了驳棱，在用了一段時間之后社搅，我實在感覺這個包沒法用……所以將其升級了一下形葬。用到了Shiny和Plotly兩個我超級喜歡的R工具笙以。
目前最新的軟件已經(jīng)發(fā)布在了Bioconductor上猖腕。我也負責(zé)回答有關(guān)于這個軟件的所有問題倘感，處理bug老玛，幫助全球科學(xué)家做分析等等蜡豹。不過個人覺得镜廉，寫這個軟件也算是一個小工程問題桨吊。
先上一張圖：
這就是最后的R包運行結(jié)果视乐。其中的GUI展示窗口佑淀。Shiny并不是什么很新的技術(shù)了伸刃，我很喜歡捧颅。只需要幾行代碼碉哑，就可以形成一個上述的網(wǎng)站扣典，非常有用贮尖，其中還可以有動態(tài)控件讓用戶調(diào)控參數(shù)湿硝。后臺的接入其實就是普通的R程序图柏，你要是有txt文件蚤吹，可以直接讀取進來裁着；你要是有數(shù)據(jù)庫二驰，可以用RSQLite等包從數(shù)據(jù)庫里直接調(diào)用……完全不受限制桶雀。另一個很好用的工具是Plotly全肮。這個工具有很多語言的接口辜腺，我純粹使用了它的R語言接口而已评疗，效果是非常好的百匆，直接生成了可交互的圖片胧华。對于R語言熟悉的同學(xué)矩动，用會這兩個包不會太久悲没。
我之所有用這個辦法寫了這個包有幾個原因：
首先這個包畢竟還是用來幫助科學(xué)家研究DNA甲基化數(shù)據(jù)的示姿，并不是什么用于Visualization的花殼子栈戳，后邊的統(tǒng)計分析嚴謹認真子檀，也是我研究領(lǐng)域一個小有名氣的軟件了亩进。
其次归薛，科學(xué)家對于可視化絕對是有要求的主籍，且不說圖片能不能達到文章要求崇猫，但是可交互的圖片，可調(diào)整的參數(shù)屋厘，可以節(jié)省人們大量的時間汗洒。
第三溢谤，有人問我為什么不像很多人一樣，做網(wǎng)頁數(shù)據(jù)庫瞻坝，讓用戶上傳數(shù)據(jù)所刀，然后再后臺運算分析（這是科研屆太常見的一種模式）浮创？我的回答是斩披，生物數(shù)據(jù)太大雏掠，一批數(shù)據(jù)如果有500個以上樣本摧玫，上傳都要好久诬像，我也不知道怎么搞到那么大內(nèi)存的服務(wù)器坏挠？也不知道如果很多人同時上傳會不會有問題？而且程序參數(shù)眾多（整個包參數(shù)加起來幾百個）邪乍，有人要調(diào)怎么辦降狠？所以我最后想到了這個辦法——在每個人自己的服務(wù)器上實現(xiàn)可視化，本地就可以基于數(shù)據(jù)生成可交互的網(wǎng)頁庇楞。這樣我不用出錢租服務(wù)器榜配，不用考慮數(shù)據(jù)大小問題吕晌，又成功把用戶體驗做出來了蛋褥。
鑒于這是CSDN博客，我之前寫了更多R語言技術(shù)睛驳，現(xiàn)在還是介紹一些R包背景烙心，也許對于國內(nèi)做計算生物分析DNA甲基化的同學(xué)有一些幫助。乏沸、
DNA甲基化指的是包裹在DNA上CG堿基（又叫CpG）的外周的一些蛋白發(fā)生變化淫茵，進而導(dǎo)致基因或者一些調(diào)控因子的表達出現(xiàn)變化，進而導(dǎo)致那些基因或者調(diào)控因子控制的表型出現(xiàn)變化蹬跃，進而導(dǎo)致疾病或者差異的發(fā)生痘昌。（沒錯，生物研究就這樣，一層關(guān)聯(lián)一層辆苔，一層決定一層算灸，非常難研究，比起寫一個網(wǎng)站直接用鼠標(biāo)點點看好不好用難太多了……）驻啤。我研究的部分就是上述鏈條的第一步：DNA最初的甲基化是如何導(dǎo)致基因或者調(diào)控因子的表達出現(xiàn)變化的菲驴。
正如大家所聽聞，測序是目前生物研究的“標(biāo)配”了骑冗。同樣科學(xué)家也可以把DNA甲基化測出來赊瞬，然后比較兩群人的差異，簡而言之：100個癌癥病人贼涩，100個正常人巧涧，把它們每一個人的DNA甲基化測出來，然后比較一下遥倦，看看人體那幾百萬個CpG位點差異如何谤绳，哪些有差異。這就是大部分生信人的工作袒哥。但是這并不容易缩筛，因為統(tǒng)計方法需要學(xué)習(xí)，數(shù)據(jù)處理需要學(xué)習(xí)堡称，如果有雜音需要消除瞎抛，如果有結(jié)果要看相關(guān)，如果有相關(guān)要看因果……所以却紧，個人覺得桐臊，計算生物學(xué)（或者生信）做得好的人，絕對在IT和統(tǒng)計界也算高手了晓殊。
目前主流的測序方法大概有全基因組測序（whole genome bisulfite sequencing）豪硅、我最常用的illumina芯片測序（Microarray-based methods），還有RRBS一類的我不知道怎么翻譯成中文……有興趣的人自己查吧挺物。
我最常用的芯片測序大概覆蓋了450,000或者850,000個位點，看你選那種芯片（illumina公司設(shè)計的飘弧，全球一個標(biāo)準识藤，誰都改不了……）。也是目前科研屆的最主流工具次伶，有點有很多啦——比如便宜痴昧，其他不重要。測序完的結(jié)果是一種叫做IDAT文件的東西冠王，這就是我程序分析的起點赶撰。這樣一個文件包含的內(nèi)容，就是我剛才說的：* 一群人中每一個的體內(nèi)的450,000或者850,000CpG的甲基化程度。*豪娜，科學(xué)家要做的餐胀，一般來說就是：
1：把它們分開，誰是癌癥瘤载，誰健康否灾？
2：整理數(shù)據(jù)，Normalization一類的跟上鸣奔。
3：開始對兩批人的幾百萬CpG位點進行比較墨技。
4：看比較出來的CpG位點哪些有差異？他們對應(yīng)到什么疾病或者通路中挎狸。(本體論oncology)
OK扣汪，簡單來說就是這么個流程。實際上很復(fù)雜的锨匆，如果不復(fù)雜崭别，這領(lǐng)域也就沒有存在的可能性。
上述是ChAMP包的主要函數(shù)和其組成的研究流程统刮。
其中藍色部分代表的是有關(guān)于數(shù)據(jù)的準備部分紊遵，就是說你要先確認數(shù)據(jù)沒問題，如果有問題侥蒙，你必須停下來暗膜，回頭找什么問題，如果處理錯了重新處理鞭衩，程序錯了該程序学搜，如果實驗做錯了……你這批幾萬美元的數(shù)據(jù)就直接可以報廢了。
其中紅色部分代表的是數(shù)據(jù)分析部分论衍，就是說確認數(shù)據(jù)沒問題之后瑞佩，你可以用它產(chǎn)出結(jié)果了，這些函數(shù)里都封裝著嚴謹?shù)慕y(tǒng)計算法坯台，涉及多重線性回歸（multiple linear regression）炬丸、網(wǎng)絡(luò)分析（network analysis）、細胞中各類分解（Cell Type Detection）等等蜒蕾。這一些函數(shù)大概集中了差不多我領(lǐng)域的所有分析角度了卸留。這就是絕大多數(shù)R包做的事情路翻。
其中黃色部門代表數(shù)據(jù)可視化部分舔株，就是說數(shù)據(jù)分析已經(jīng)完了房铭，你可以用它查看分析結(jié)果，自動生成可交互撤摸，可下載的圖片毅桃。
箭頭指示了分析流程褒纲，首先是load數(shù)據(jù)，然后是QC（quality control）钥飞，然后是normalization莺掠，然后是SVD分析（看有沒有batch effect）。準備完畢之后代承，黑點代表了一批質(zhì)量有保證汁蝶，經(jīng)過處理，可以直接上馬進行數(shù)據(jù)分析的數(shù)據(jù)了论悴。之后的分析掖棉，DMP代表找出Differential Methylation Probe（差異化CpG位點），DMR代表找出Differential Methylation Region（差異化CpG區(qū)域）膀估，Block代表Differential Methylation Block（更大范圍的差異化region區(qū)域）幔亥，RefFree代表細胞差異被修正過后再找的DMP，EpiMod是基于基因作用網(wǎng)絡(luò)的差異化分析察纯。（看到木有帕棉，整個生信研究，就是各種找兩群人間的差異……但是真心很不好做饼记。）
首先通過Bioconductor安裝程序香伴，直接在R中輸入：
source("https://bioconductor.org/biocLite.R")biocLite("ChAMP")
安裝如果報錯很常見，因為我引用了很多其他包具则，而那些包也用到了很多其他包即纲，偶爾會缺失。只需要認真看一下錯誤最后幾行博肋，會顯示出哪個包缺失了低斋，然后重新安裝那個就行了。
程序提供了兩批測試數(shù)據(jù)作為演示匪凡，一批450K的數(shù)據(jù)膊畴，和一批我自己模擬的850K數(shù)據(jù)（沒錯，數(shù)據(jù)分析程序的開發(fā)人自己沒有數(shù)據(jù)用于開發(fā)……）病游。數(shù)據(jù)在我附帶的ChAMPdata包中唇跨。
下面開始一套完成的DNA甲基化數(shù)據(jù)分析流程：
程序過程中會跑出很多中間結(jié)果，有些花時間很長衬衬，比如幾十分鐘买猖，測試數(shù)據(jù)的所有中間結(jié)果，都已經(jīng)存儲在了包中佣耐，可以用下面代碼讀取：
data(testDataSet)
比如只對于Shiny或者Plotly感興趣的同學(xué)唧龄，可以直接用上面代碼讀取數(shù)據(jù)兼砖，調(diào)到文章后部的GUI程序部分奸远，開始把玩看效果。
3.1 首先是數(shù)據(jù)Load
下面的testDir指向的是我提供的450K測試數(shù)據(jù)的地址讽挟，champ.load()函數(shù)可以直接從目錄讀取數(shù)據(jù)懒叛，你只需要把illumina機器跑出來的IDAT文件全部解壓放在目錄下就行了，但是目錄下還必須放一個有關(guān)于樣本信息的csv文件耽梅，也是illumina機器產(chǎn)生的薛窥，這很重要。程序再智能眼姐，一開始的這點小要求還是必須要滿足的诅迷。
testDir=system.file("extdata",package="ChAMPdata")
myLoad <- champ.load(testDir,arraytype="450K")
如果是自己的文件，那么需要把所有的idat和csv文件众旗，放在一個文件夾里面罢杉，然后R的工作目錄定義到該文件夾(自定義為AA）
myLoad <- champ.load(directory="AA", arraytype="450K")
這樣就把數(shù)據(jù)全部讀入了，其中會做一部分檢測贡歧，比如低質(zhì)量的CpG位點滩租，也會做一些過濾，比如位于性染色體上的CpG位點會被刪掉利朵，比如SNP位點（就是基因突變位點）附近的CpG位點會不穩(wěn)定等等……詳細原因算法我就不多說了律想，說起來沒完沒了，一個原因就是一兩篇論文绍弟。
然后我提供了一個小程序CpG.GUI來看一下你這批數(shù)據(jù)的所有CpG的分布技即。這個小程序是我一開始學(xué)習(xí)Shiny和Plotly的時候測試代碼，一直留下來了晌柬。這也是所有GUI里最簡單的一個姥份，最適合用于分析代碼結(jié)構(gòu)，效果如下：
顯示的東西都是年碘，你數(shù)據(jù)的CpG位點是如何覆蓋的澈歉，比如第一張圖（左上）是所有CpG在染色體上的分布，等等……代碼太長屿衅，而且不易讀（不是我的錯埃难，Shiny這個框架寫東西縮進空格太多，結(jié)構(gòu)類似于HTML那種涤久。） 大家可以看Github涡尘。
3.2 下一步做的是Quality Control
代碼如下：
QC.GUI(beta=myLoad$beta,arraytype="450K")
然后會出現(xiàn)一系列的圖像：
會有一個網(wǎng)頁自動打開，其中有幾個tab响迂，每一個tab代表了一些分析圖像考抄，用戶可以點擊那些tab，每一次轉(zhuǎn)到一個圖像上蔗彤，需要等待幾秒鐘川梅，后臺在計算疯兼，計算完畢后，圖像會自動生成贫途。圖像用的是plotly吧彪，可交互，鼠標(biāo)放上去會有顯示丢早。
這一個函數(shù)的功能是告訴用戶姨裸，
tab1: 樣本之間的Clustering情況是怎樣的？
tab2: 測序芯片上的兩種測序技術(shù)（Type-I和Type-II）之間差異如何怨酝？
tab3: 所有樣本的分布式怎樣的傀缩？
tab4: 根據(jù)當(dāng)前數(shù)據(jù)進行了層聚類樹狀圖是怎樣的。
tab5: 這批數(shù)據(jù)中差異性最大的一批CpG的熱力圖凫碌。
3.3 然后做的是Normalization
Normalization就是為了把上面的tab2里的兩種測序方法得到的數(shù)據(jù)不一致現(xiàn)象消除扑毡，否則如果你通過分析函數(shù)知道CpG有差異，你如何知道這個差異是疾病帶來的盛险？還是測序方法差異帶來的瞄摊？這一部至關(guān)重要，現(xiàn)在已經(jīng)有很多方法用于完成這個問題苦掘，還有科學(xué)家在樂此不疲地開發(fā)更多算法完成這一步换帜。
myNorm <- champ.norm(beta=myLoad$beta,arraytype="450K",cores=5)
規(guī)則化之后的結(jié)果是顯而易見的，type2-BMIQ是規(guī)則化之后的分布曲線鹤啡，和紅色的已經(jīng)一致了惯驼，所以測序方法上的差異，不會導(dǎo)致結(jié)果出現(xiàn)問題递瑰。
3.4 然后做的是SVD
這個分析講真很簡單祟牲，但是確實很有用，算法一句話形容就是：把數(shù)據(jù)進行SVD分解抖部，查看Component有沒有與某些重要的Factor相關(guān)说贝。SVD大家都知道，主成分分析第一步就是SVD慎颗，其性質(zhì)類似于把重要的Component給找出來了乡恕。而與Factor的相關(guān)，可以告訴使用者俯萎，你的數(shù)據(jù)中的大部分差異傲宜，是又你想研究的因素（比如癌癥、疾病夫啊、年齡等等）導(dǎo)致的函卒，還是由你不想研究的東西，比如樣本的種族撇眯、數(shù)據(jù)產(chǎn)出的研究所等等一些無關(guān)痛癢的東西決定的报嵌。如果是后者躁愿，那很不幸你的數(shù)據(jù)質(zhì)量也不高……
champ.SVD(beta=myNorm,pd=myLoad$pd)
比如上述我這批測試數(shù)據(jù)的結(jié)果就非常好：因為數(shù)據(jù)中的差異，顯著與樣本的PhenoType互相關(guān)聯(lián)沪蓬。換言之，你如果做針對Sample_Group（就是Cancer和Normal）的分析来候，得到的結(jié)果就是由于疾病導(dǎo)致的跷叉。
額外提一句，如果數(shù)據(jù)與其他無關(guān)的東西相關(guān)营搅，而與你想要研究的東西不相關(guān)云挟，那么你就需要想很多辦法了，比如转质，如果上述的圖片中园欣，紅色那個板塊（顯著相關(guān)的意思，越紅越相關(guān)）是在Array那一行休蟹，而不是在Sample_Group那一行沸枯，那就意味著你的數(shù)據(jù)中的差異更多是Array導(dǎo)致的，而Array是芯片中的樣本排列方式赂弓，換言之就是……你的實驗不夠好绑榴。但有些問題很難避免，比如你的100個病人年齡不同很正常盈魁，老人身體上自然會有一些弱化翔怎，年輕人可能更好，樣本的年齡很自然地會引入一些誤差杨耙。
上述這種情況一旦遇到是非常痛苦的赤套，至今也有很多我領(lǐng)域的人在研究如何把這樣的錯誤修正，如何矯正數(shù)據(jù)等等珊膜。我程序集成了一個別人寫好的用于修正這一類錯誤的函數(shù)Combat()容握。我個人覺得效果一般，而且時間非常非常非常長（該測試數(shù)據(jù)只有8個樣本辅搬，要跑一個半小時Ｎň凇）。但是似乎也沒有更好的辦法……
如果想要用我程序集成的函數(shù)堪遂，用如下代碼修正：
myCombat <- champ.runCombat(beta=myNorm,pd=myLoad$pd,batchname=c("Array"))
1
3.5 尋找差異化CpG
如果上述步驟都能完成介蛉，說明你的數(shù)據(jù)質(zhì)量是挺好的，可以分析了溶褪！最簡單币旧，但同時最重要的一種分析就是差異化CpG，正如我之前說的猿妈，這個步驟的目的就是吹菱，找出幾百萬CpG中的哪些在疾病中發(fā)生了變化巍虫，而這些變化又是如何導(dǎo)致了基因發(fā)生了變化，最終導(dǎo)致了人體生病鳍刷。而做的方法直觀上簡單的可怕：
你有100個癌癥病人占遥，100個正常人，每個人身體中都有450K個CpG的位點在測序出來输瓜，針對其中的每一個CpG瓦胎，你都有200個數(shù)據(jù)對不對？如果這一個CpG在100個正常人中和100個癌癥病人中的甲基化水平都差不多尤揣，你還會繼續(xù)懷疑它嗎搔啊，當(dāng)然不會！
但如果北戏，你的100個癌癥病人普遍在這一個CpG上的甲基化水平高（不太嚴謹?shù)呛苄蜗蟮卣f负芋，就是DNA那個CpG的序列外層越來越多的部分被甲基附集上去了），而那100個普通人的甲基化水平不高嗜愈，那這個位點就很有嫌疑了對吧旧蛾。
為什么需要一定量的樣本，比如100個蠕嫁？因為如果你只找兩個人蚜点，一個德國人一個中國人，那個德國人高拌阴，那個中國人矮绍绘，你能因此就說德國人在人種上比中國人高嗎？當(dāng)然不能……但如果你找到100個德國人迟赃，在找到100個中國人陪拘，比較以后，就比較可信了纤壁，這涉及t.test()里的power的問題左刽。
這就是做研究的艱難：你得找到100個病人，讓他們同意治療酌媒，然后耗資幾萬幾十萬完成測序欠痴，有了數(shù)據(jù)還要祈禱實驗員沒有點錯試管，數(shù)據(jù)下來了如果由于年齡秒咨、人種等原因喇辽，數(shù)據(jù)差異已經(jīng)找不到了，你還得想辦法修正這些問題雨席，然后你可以開始比較菩咨，從幾百萬位點（基因、SNP、CpG都一樣）中找出那些可能有關(guān)系的……等你做完這一切抽米，可能一兩年已經(jīng)過去了特占，而你本科畢業(yè)就進入IT界的同學(xué)可能已經(jīng)工資兩萬多了，這還絕對算是科研中很快節(jié)奏的項目了云茸，所以我覺得社會真的應(yīng)該給科研工作者更多尊重是目，做最難的活，拿最低的工資标捺。
代碼如下：
myDMP <- champ.DMP(beta = myNorm,pheno=myLoad$pd$Sample_Group)
我的程序函數(shù)之間的接口是設(shè)計過的胖笛，就是說，你從上一步跑出來的程序宜岛，可以自動被下一步識別到（當(dāng)然你要是改名字什么的那就找不到了，只能自己在參數(shù)中設(shè)定）功舀。完成這一步以后萍倡，就找出了幾百萬位點中可能有問題的。做這一步我最討厭的就是：有時候一把跑出來辟汰，幾百萬位點中列敲，幾萬甚至十幾萬都是顯著的（做年齡的時候就這樣，因為年齡對于人身體的影響太大帖汞，可謂是全方位無死角的）戴而！
然后是整個包中的最復(fù)雜的GUI函數(shù)：
DMP.GUI(DMP=myDMP,beta=myNorm,pheno=myLoad$pd$Sample_Group)
DMP.GUI(DMP=myDMP[[1]],beta=myNorm,pheno=myLoad$pd$Sample_Group)
會出現(xiàn)一個稍微復(fù)雜一點的網(wǎng)頁，如下：
這個頁面左邊翩蘸，是一些可以調(diào)控的參數(shù)控件所意。比如你可以選擇你想要的位點的p value等等，然后你點擊Sumbit催首，右邊就會出現(xiàn)通過了參數(shù)篩選的所有CpG的信息扶踊。這比起以往的方法方便太多了。
然后如果你選擇第二個Tab郎任，會出現(xiàn)另一個界面(自動的)秧耗。
第二個tab是前一個頁面的所有顯著CpG的熱烈圖。每一行是一個位點舶治，每一列是一個樣本分井，可以很明顯地看出，通過統(tǒng)計分析找到的人體中幾百萬CpG位點中最顯著差異的那些位點霉猛，在癌癥和正常人的甲基化水平完全不同尺锚。
第三個tab如下：
這個中的三個柱狀圖是所找到的顯著CpG的分布，和之前的CpG.GUI類似惜浅。這張圖充分證明了Plotly的強大缩麸，每一個柱子上有三個bar條對吧，你可以點擊右邊的legend把它關(guān)掉。
下一張是整個R包實現(xiàn)起來最復(fù)雜的一個頁面：
用戶可以直接輸入基因的名字來選擇你們想要查詢的基因杭朱，看這些基因里是不是收到了”癌變CpG“的影響阅仔。比如，在圖中的左側(cè)輸入NFIX弧械，點擊Submit八酒，右側(cè)需要計算十秒鐘左右，然后上邊會顯示出這個基因中顯著差異CpG刃唐。下面會顯示出基因的不同區(qū)域羞迷，比如是基因開始，還是在基因后端一類的画饥。我例子中這個基因衔瓮，有大約十多個CpG位點，而且這些位點無一例外有顯著差異抖甘，所以基本可以確定這個基因與癌癥關(guān)系有關(guān)系热鞍。（但是具體關(guān)系是什么很難確定，誰導(dǎo)致了誰（因果）都不一定衔彻，這就必須要去實驗室點試管做實驗才能知道了薇宠。），在下方艰额，我用R語言爬網(wǎng)站的技術(shù)澄港，抓去了這個基因在一個比較權(quán)威的數(shù)據(jù)網(wǎng)站wegene的信息，這樣的話柄沮，用戶一下點擊回梧，就可以同時看出分析圖和基因說明，而不用再去查網(wǎng)頁看這個基因是干嘛的祖搓。
抓取網(wǎng)頁的代碼如下：
webpage <- tags$iframe(src=paste("https://www.wikigenes.org/e/gene/e/",MatchGeneName[which(MatchGeneName[,1]==Gene_repalce()$genename),2],".html",sep=""),width="60%",height="100%")
不過還需要在plotly中調(diào)用htmlOutput()函數(shù)創(chuàng)建一塊用于網(wǎng)頁展示的空頁漂辐。詳細的還是看代碼吧……
基本上這樣的圖，已經(jīng)足夠發(fā)表文獻了棕硫，直接點擊右上角下載髓涯，就可以作為論文圖片的，所以軟件發(fā)布以后哈扮，很多人都說很方便纬纪，這也讓我比較開心，我沒水平做出解決癌癥的方案滑肉，療法包各，至少可以讓全球攻堅這方面的人做得更快更舒心一些吧。
之后的一個tab是：
圖片上半部是根據(jù)tab1選擇出的顯著CpG針對基因的富集情況做的排序圖靶庙。（否則人體中有兩萬多基因问畅，用戶怎么知道應(yīng)該在上一個tab搜索哪一個呢？）下圖就是一個最簡單不過的boxplot，但是它依然是目前科研屆最常見最有用的圖护姆。用戶可以選擇左邊的cpg來選擇想要查看的CpG矾端。
有些需要做羥甲基化的，可以使用下面代碼：
myDMP<-champ.DMP(beta=myNorm,pheno=myLoad$pd$Sample_Group,compare.group=c("oxBS","BS"))
hmc <- myDMP[[1]][myDMP[[1]]$deltaBeta>0,]
# Then you can use above code to extract hydroxymethylation CpGs.
3.6 尋找差異化區(qū)域
Differential Methylation Block(DMB) 和 Differential Methylation Region(DMR)這兩個概念是之前的科學(xué)家提出來的卵皂，就是因為找單一的CpG已經(jīng)不夠顯著了秩铆！比如一個基因跨距幾千位點，通過上一步分析灯变，就一個CpG有問題殴玛，你算他有問題嗎？這就是一個比較尷尬的問題……所以有的科學(xué)家就覺得添祸，單一CpG影響可能不可靠（感覺上有點像木村資生的中性理論滚粟，覺得有些位點為人太隨意了，而有些真正受重大影響的基因刃泌，其實他們的基因序列上會受到一系列的暴擊傷害—— 一連串的CpG都會出現(xiàn)很明顯的差異凡壤，既然這樣，我們就專注于研究那些比較要緊的基因吧蔬咬，有些暫時看不出來的就不管了，比如那些零零散散的差異CpG就當(dāng)它不存在吧沐寺。
所以林艘，比如先用代碼：
myDMR <- champ.DMR(beta=myNorm,pheno=myLoad$pd$Sample_Group,method="Bumphunter")
上面是我寫的DMR函數(shù)，其中集成了三個方法混坞，DMRcate,Bumphunter狐援，和我自己寫的ProbeLasso。至于誰好……互相都說自己是最好的究孕，經(jīng)過我之前用模擬數(shù)據(jù)的評測啥酱，Bumphunter比較可靠，精確度可以有90%以上厨诸，ProbeLasso才有75%左右吧镶殷，DMRcate是我后來集成進去的，沒有評測過微酬。
然后我們可以用DMR的GUI函數(shù)查看一下結(jié)果（有三個tab绘趋，我只放最重要的那個）：
輸入代碼以后，估計需要等待大約20-30秒颗管，那一段時間陷遮，程序正在后臺完成Mapping，把找到的CpG匹配到全基因組上垦江。然后就可以像之前的那個DMP.GUI一樣操作了帽馋。圖中可以看出，DMR就是一個連續(xù)不斷都比較長的差異片段，科學(xué)家們覺得绽族，這樣的連續(xù)差異片段姨涡，對于基因的影響會更加明顯，只找這樣的片段项秉，可以使得計算生物學(xué)的打擊精度更為準確绣溜，也可以讓最終找出來的結(jié)論數(shù)據(jù)更少，便于實驗人員篩選娄蔼。
而其實找出來以后怖喻，差別就不大了，你依然需要把這一些區(qū)域Mapping到基因組上岁诉，找到對應(yīng)的基因锚沸，然后該做什么繼續(xù)做，又回到普通的比較模式了涕癣。這個網(wǎng)頁的其他兩個tab自動完成了匹配哗蜈，這就是一些很簡單的匹配工作而已，不過畢竟省了一些時間坠韩。
另外一個類似的東西就是DMB距潘，那個東西出現(xiàn)的原因是，有的科學(xué)家覺得只搁，DMR這樣的區(qū)域還不夠顯著音比，DNA上的甲基化出現(xiàn)變化，可能是綿延幾千位點的氢惋！而且只會在基因以外的區(qū)域洞翩，但是這些基因以外的區(qū)域發(fā)生變化，卻會導(dǎo)致基因的表達發(fā)生變化焰望。你可以想象成骚亿，北京周邊的河北在大煉鋼鐵，然后北京也跟著霧霾了熊赖，大概就是這意思来屠。
我個人對這種假設(shè)還是有點懷疑的……根據(jù)這個提出這個理論的作者的文章，他通過找到這種Block之后震鹉，將它與一些癌癥形狀進行了correlation的妖，結(jié)果是顯著的，然后這種理論大概就有了足陨，相應(yīng)的算法也出現(xiàn)了嫂粟，雖然并沒有見到這樣的算法確實找到了某些重要的生物突破，但是身為這方面的人墨缘，我也不得不實現(xiàn)了他這個算法星虹，放在這里供用戶使用零抬。
我個人覺得correlation出來的結(jié)果，距離因果還是有很大距離的宽涌。建立在Correlation上的假設(shè)理論平夜，甚至算法，實在是有點懷疑它的統(tǒng)計power卸亮。
我感覺得做計算生物很大一個問題就是不要統(tǒng)計過度忽妒，比如我根據(jù)相關(guān)性得出一些結(jié)論，我在用這批結(jié)論去做另一個聚類或者相關(guān)兼贸，在得到一個統(tǒng)計值結(jié)論段直。這樣的power減弱的很快啊。你做一次統(tǒng)計溶诞，錯誤率可能只是20%鸯檬，那么根據(jù)其結(jié)果再錯一次統(tǒng)計，錯誤率就會上升到 1-80%*80%=36%螺垢！除非在一次統(tǒng)計過后喧务，用相對嚴格的標(biāo)準去篩選統(tǒng)計值，也許有助于這種情況好轉(zhuǎn)枉圃。
Differential Methylation Block的代碼和圖都不展示了功茴，程序的Vignette里邊有。
3.7 尋找作用通路網(wǎng)絡(luò)中的疾病關(guān)聯(lián)小網(wǎng)絡(luò)
這貨根本就沒有中文譯名孽亲，上邊是我自己翻譯的坎穿。簡而言之就是，人體內(nèi)的作用網(wǎng)絡(luò)實在是太多了墨林，幾百幾千吧赁酝，每一個都涉及了一系列基因犯祠。這就是過往的科學(xué)家研究出來的旭等，比如某個通路會導(dǎo)致頭疼，然后有幾百個蛋白（及其背后的基因）都是通過共同作用導(dǎo)致了這一場頭疼的衡载。但是如果你頭疼搔耕，不見得是所有這些基因都出了問題，而是可能其中的某一部分痰娱，甚至于只有一兩個出了問題弃榨。所以你就可以基于已經(jīng)存在的那個網(wǎng)絡(luò)，再集合數(shù)據(jù)梨睁，找出哪些網(wǎng)絡(luò)可能出問題了鲸睛？或者是這些大網(wǎng)絡(luò)中的哪些基因具體除了問題。
這個功能很重要坡贺，否則做完上邊幾步官辈，用戶只會知道那些基因有問題箱舞，至于他們之間有沒有關(guān)系，是不是會同時作用于某些網(wǎng)絡(luò)拳亿，就沒法知道的了晴股。但其實也不復(fù)雜，只需要自己寫個程序肺魁，匹配一樣基因和網(wǎng)絡(luò)就完了电湘，只不過數(shù)據(jù)的準備啊，洗啊鹅经，匹配啊寂呛，也是夠煩的，所以這個函數(shù)就提供了全套的分析瞬雹。
myEpiMod <- champ.EpiMod(beta=myNorm,pheno=myLoad$pd$Sample_Group)
上述代碼之后可以輸出一些圖片到當(dāng)前目錄的文件夾下：
其中的顏色代表了基因的甲基化上下調(diào)的趨勢昧谊。（人們得了某個病，不一定是導(dǎo)致那個病的基因高表達了酗捌，也有可能是抑制那種病的某個基因低表達了呢诬，甚至于可能是抑制那個抑制生病的基因的基因突然高表達了，這可以無窮下去……）胖缤。所以一個網(wǎng)絡(luò)了尚镰，不見得都是同樣的表達趨勢，藍點代表甲基化高表達哪廓，而黃點代表低表達狗唉。一言以蔽之，可能是這個網(wǎng)絡(luò)中的這一些基因”同謀“導(dǎo)致了你在研究的這個疾病發(fā)生了涡真。
3.8 拷貝數(shù)變異
這個的意思就是可能你的測序數(shù)據(jù)里分俯，表現(xiàn)出某些疾病導(dǎo)致的原因，是有些基因片段被復(fù)制的此處過多或者過少哆料。檢測方法很聰明缸剪，基因人們是看不見的，但是測序的原理，是用一些片段去把與CpG有關(guān)的片段”粘“下來。然后對這些片段進行測序兵睛。一般來說，我是不關(guān)心粘下來多少的奋渔，因為我研究的基因的甲基化水平，其實就是某一個為點的被甲基化和未被甲基化水平的比例壮啊。例如嫉鲸，假設(shè)正常的情況下，某個CpG所在片段會被”粘“下來一萬次歹啼，其中8000片上的CpG都沒有被甲基化玄渗，而2000片被甲基化了减江。那么最后到我手里的數(shù)據(jù)，就是這個CpG的甲基化值為2000/8000=0.25捻爷。如果被粘下來兩萬次辈灼，那么必然是16000片上沒有被甲基化，而4000片上被甲基化了也榄，我得到的比例還是0.25巡莹。完全感覺不到區(qū)別啊甜紫！
所以檢測辦法降宅，就是去看，到底多少片被拉下來了囚霸，然后與正常人作比較腰根。這個函數(shù)不是我寫的，是以前一個教授拓型，但是我覺得有點粗糙额嘿，因為它只能統(tǒng)計出每一條染色體上的拷貝數(shù)變化差異，我覺得這不夠的劣挫。如果有機會册养，我希望更新一下程序，讓它打擊更為定點压固，就是找出來是哪些CpG發(fā)生了突變導(dǎo)致了拷貝數(shù)變異球拦，然后與其他的指標(biāo)，比如表達值一類的作比較帐我，這樣應(yīng)該會有新發(fā)現(xiàn)坎炼。
myCNA <- champ.CNA(intensity=myLoad$intensity,pheno=myLoad$pd$Sample_Group)
代碼過后，只會生成一些這樣的圖片：
圖的大概意思就是拦键，C這個疾病情況下谣光，相比于正常人，在全基因上那些位置發(fā)生拷貝數(shù)變異矿咕。我從來沒有用過這個功能抢肛，感覺精度根本不夠狼钮，我就努力確保了它程序沒問題碳柱，運轉(zhuǎn)順利。
3.9 細胞種類分離問題
這個堪稱是我這個研究方向上目前全球在最努力攻關(guān)的問題……我做它快三年了熬芜，也完全沒有實質(zhì)性突破和進展莲镣，好在全球范圍內(nèi)好像也沒有人做出來。簡而言之是問題是這樣的涎拉，人體內(nèi)充斥著各種細胞瑞侮，你給病人采血取樣的圆，其實取到的也是細胞的組合。然后問題就來了半火，你確定某個疾病是所有的細胞都出問題了嗎越妈？還是某幾種細胞出了問題？但是你的數(shù)據(jù)是細胞們混合在一起的钮糖，你怎么分開梅掠？比如下面的公式：
X=C.F
這里的X
是你的身體中采集到了CpG甲基化值，或者基因表達值店归，或者各種亂七八糟的數(shù)值（只要不是基因序列本身阎抒，應(yīng)該都存在這個問題吧，基因序列確實穩(wěn)定消痛，不會發(fā)生變化且叁，這也就是為什么SNP測序可以發(fā)展成為產(chǎn)業(yè)，而其他的不可能秩伞。）逞带。這個X其實是又兩個部分組成的：C和F。C是每一個位點（或者基因）在單一細胞中的表達值纱新，比如你的血液里有紅細胞掰担，白細胞和血小板，C就是每一個位點在這些細胞中的表達值矩陣怒炸，F(xiàn)
是你采集的樣本中带饱，各種細胞的比例。比如你的血液中紅細胞阅羹，白細胞和血小板的比例勺疼。
這帶來的問題就是。每個人體內(nèi)的細胞比例是不同的捏鱼。甚至于执庐，有些疾病本身就是因為細胞比例發(fā)生變化導(dǎo)致的，血小板太少會導(dǎo)致流血不止對吧导梆？進而轨淌，導(dǎo)致疾病發(fā)生的基因也是在不同細胞中的，如果他們混在一起看尼，你怎么知道递鹉，是誰具體導(dǎo)致了癌癥？其性質(zhì)就像是之前在SVD的時候我提過的干擾因素藏斩，你找了兩群人想找癌癥和正常人的區(qū)別躏结，結(jié)果一群人是老人，一群人是小孩狰域，這怎么可能找的對媳拴？年齡還可以修正黄橘，只要你知道大家的年齡就可以修正。如果你找的100個樣本血細胞混合是不同的屈溉，你怎么繼續(xù)做研究塞关，而且你也不知道大家的細胞比例是什么，連修正的辦法都沒有子巾。
所以就是X
是你的數(shù)據(jù)描孟，C是純化過的單一細胞內(nèi)的甲基化值，F(xiàn)是樣本的細胞混合比例砰左。X是C和F的矩陣乘匿醒。目的就是為了識別出F
。
這個問題還要分為兩個方向缠导，是關(guān)于C
的廉羔，很明顯，在上邊的算式中僻造，C和F兩個矩陣決定了最終的數(shù)據(jù)矩陣X憋他，那么只要有一個，就能大概解出另一個髓削。換言之竹挡，如果你知道一群人的數(shù)據(jù)里的細胞比例，那么你就可以推測出純化過的細胞甲基化值立膛【竞保或者如果你知道純化過的細胞甲基化值，我也可以推測出F
——人群的細胞比例宝泵。
上述這種情況已經(jīng)有解好啰，而且解很不錯，Houseman教授的Refbase算法很好用儿奶，他的辦法就是框往，在你知道C
的情況下，解開F
闯捎，然后該修正修正椰弊，該分析分析。
代碼：
myRefBase<- champ.refbase(beta=myNorm,arraytype="450K")
1
該函數(shù)可以在已知C的情況下瓤鼻，解出F秉版。效果很好，但是我也沒有什么可以展示的娱仔，因為這畢竟不是一個什么可以可視化的過程沐飘。
真正困難的是根據(jù)X進行C和F的同時估測游桩，這個到目前為止似乎沒有什么太好用的辦法牲迫。Houseman教授提出了替代變量發(fā)去將已知的一些變量替代為PCA估測出來的Component耐朴，并且默認這些Component中有細胞比例。這個感覺有道理盹憎，但是應(yīng)該對數(shù)據(jù)質(zhì)量要求很高陪每，如果你還有其他可能的沒有消除干凈的變量，就會導(dǎo)致你的Component其實根本不是細胞比例……然后就錯掉了。
他的函數(shù)我也集成進來了灌灾，但是不見得很好用：
myRefFree <- champ.reffree(beta=myNorm,pheno=myLoad$pd$Sample_Group)
1
大概就是覺得自己可以嘗試開一個技術(shù)博客，記錄一下平時用到的技術(shù)轴总。這個項目是去年下半年做的一個小東西往堡，不是什么正式項目痹兜，科研屆不會在乎什么用戶可視化的無聊軟件的，花時間也不長。發(fā)在CSDA上主要是覺得Shiny和Plotly挺好用的，也許可以給對R語言有興趣的同學(xué)一些啟發(fā)和靈感。我學(xué)的也不好肮之，尤其有些算法在深度和概念上理解可能有錯，歡迎同學(xué)們指正柑土。
最后總結(jié)幾句：計算生物學(xué)真的不好做狐榔，對統(tǒng)計編程要求挺高的庵楷，而且很多技術(shù)可以與很多領(lǐng)域?qū)印１热缱钚路脚d未艾的那個”Data Science“看疙。
搞這行的人，至少接觸過幾個G以上的數(shù)據(jù)脚线，這些數(shù)據(jù)說太大不至于搁胆，但是能把它們處理清楚，意味著他處理百萬行以下的超小數(shù)據(jù)集是沒問題的邮绿。
做生信統(tǒng)計嚴謹要求很高渠旁，每一個同學(xué)都必須搞清楚什么時候用Wilcox，什么時候用t.test()船逮，都能知道correlation不等于因果顾腊，所以不會犯里約奧運會時候那種幼稚的錯誤：”里約奧運會是社交網(wǎng)絡(luò)上被談?wù)摰淖疃嗟囊粚茫 埃峙虏皇菉W運會進步了挖胃，而是社交網(wǎng)絡(luò)發(fā)展了）杂靶。
有些IT技術(shù)挺好的，但是一直進不到科研屆來酱鸭，我曾經(jīng)需要跑一個數(shù)據(jù)吗垮，用R不可能跑完，用Python估計兩三天凹髓，最后用了Spark烁登，3分鐘出結(jié)果……那個時候我就覺得，領(lǐng)域之間應(yīng)該互相多學(xué)習(xí)學(xué)習(xí)蔚舀。金融IT就業(yè)單位防泵，不要看到掛著生物兩個字的人，就推到門外蝗敢，這其實可能是世界上最接近”大數(shù)據(jù)“和”數(shù)據(jù)科學(xué)“的專業(yè)了捷泞。
生物數(shù)據(jù)和金融數(shù)據(jù)差別能有多大呢？如果細胞識別算法能出來寿谴，是不是一只股票的投資維度分析也會發(fā)生變化锁右？通信問題中的音頻混合問題是不是可以被解決？Shiny和Plotly這個技術(shù)用于數(shù)據(jù)展示不也很方便嗎？SVD分解針對各種因素的Correlation與社會科學(xué)家分析城市省份經(jīng)濟發(fā)展的方式應(yīng)該是一致的吧咏瑟？生物數(shù)據(jù)有干擾因素的性質(zhì)拂到，與海洋大氣采集的數(shù)據(jù)干擾有沒有共通性？問題往往是相通的码泞，技術(shù)更經(jīng)常是一樣兄旬。
