原貼地址：https://blog.csdn.net/joshua_hit/article/details/54982018
?????????????? https://bioconductor.org/packages/release/bioc/vignettes/ChAMP/inst/doc/ChAMP.html
我目前的研究領域是DNA甲基化，這個領域并沒有很好的分析軟件，頂多一些R包吧，而且很多學生物的人，想要分析DNA甲基化數(shù)據(jù)并不容易，雖然這樣一些包已經(jīng)存在的，但是首先，他們不一定易用，例如minfi（我領域第一號分析包），開發(fā)人在寫的時候，就已經(jīng)默認使用者有相當?shù)腞語言知識和統(tǒng)計知識。其次是，包的功能分散，并 不集中，經(jīng)常需要從一個軟件跑出來的東西，弄來弄去導入到其他軟件中去進行下一步分析。
我在英國留學的時候，遇到了ChAMP包的原作者，她已經(jīng)不想再管這個軟件了。所以我接手了，在用了一段時間之后，我實在感覺這個包沒法用……所以將其升級了一下。用到了Shiny和Plotly兩個我超級喜歡的R工具。
目前最新的軟件已經(jīng)發(fā)布在了Bioconductor上。我也負責回答有關于這個軟件的所有問題，處理bug，幫助全球科學家做分析等等。不過個人覺得，寫這個軟件也算是一個小工程問題。
先上一張圖：
這就是最后的R包運行結(jié)果。其中的GUI展示窗口。Shiny并不是什么很新的技術了，我很喜歡。只需要幾行代碼，就可以形成一個上述的網(wǎng)站，非常有用，其中還可以有動態(tài)控件讓用戶調(diào)控參數(shù)。后臺的接入其實就是普通的R程序，你要是有txt文件，可以直接讀取進來；你要是有數(shù)據(jù)庫，可以用RSQLite等包從數(shù)據(jù)庫里直接調(diào)用……完全不受限制。另一個很好用的工具是Plotly。這個工具有很多語言的接口，我純粹使用了它的R語言接口而已，效果是非常好的，直接生成了可交互的圖片。對于R語言熟悉的同學，用會這兩個包不會太久。
我之所有用這個辦法寫了這個包有幾個原因：
首先這個包畢竟還是用來幫助科學家研究DNA甲基化數(shù)據(jù)的，并不是什么用于Visualization的花殼子，后邊的統(tǒng)計分析嚴謹認真，也是我研究領域一個小有名氣的軟件了。
其次，科學家對于可視化絕對是有要求的，且不說圖片能不能達到文章要求，但是可交互的圖片，可調(diào)整的參數(shù)，可以節(jié)省人們大量的時間。
第三，有人問我為什么不像很多人一樣，做網(wǎng)頁數(shù)據(jù)庫，讓用戶上傳數(shù)據(jù)，然后再后臺運算分析（這是科研屆太常見的一種模式）？我的回答是，生物數(shù)據(jù)太大，一批數(shù)據(jù)如果有500個以上樣本，上傳都要好久，我也不知道怎么搞到那么大內(nèi)存的服務器？也不知道如果很多人同時上傳會不會有問題？而且程序參數(shù)眾多（整個包參數(shù)加起來幾百個），有人要調(diào)怎么辦？所以我最后想到了這個辦法——在每個人自己的服務器上實現(xiàn)可視化，本地就可以基于數(shù)據(jù)生成可交互的網(wǎng)頁。這樣我不用出錢租服務器，不用考慮數(shù)據(jù)大小問題，又成功把用戶體驗做出來了。
鑒于這是CSDN博客，我之前寫了更多R語言技術，現(xiàn)在還是介紹一些R包背景，也許對于國內(nèi)做計算生物分析DNA甲基化的同學有一些幫助。、
DNA甲基化指的是包裹在DNA上CG堿基（又叫CpG）的外周的一些蛋白發(fā)生變化，進而導致基因或者一些調(diào)控因子的表達出現(xiàn)變化，進而導致那些基因或者調(diào)控因子控制的表型出現(xiàn)變化，進而導致疾病或者差異的發(fā)生。（沒錯，生物研究就這樣，一層關聯(lián)一層，一層決定一層，非常難研究，比起寫一個網(wǎng)站直接用鼠標點點看好不好用難太多了……）。我研究的部分就是上述鏈條的第一步：DNA最初的甲基化是如何導致基因或者調(diào)控因子的表達出現(xiàn)變化的。
正如大家所聽聞，測序是目前生物研究的“標配”了。同樣科學家也可以把DNA甲基化測出來，然后比較兩群人的差異，簡而言之：100個癌癥病人，100個正常人，把它們每一個人的DNA甲基化測出來，然后比較一下，看看人體那幾百萬個CpG位點差異如何，哪些有差異。這就是大部分生信人的工作。但是這并不容易，因為統(tǒng)計方法需要學習，數(shù)據(jù)處理需要學習，如果有雜音需要消除，如果有結(jié)果要看相關，如果有相關要看因果……所以，個人覺得，計算生物學（或者生信）做得好的人，絕對在IT和統(tǒng)計界也算高手了。
目前主流的測序方法大概有全基因組測序（whole genome bisulfite sequencing）、我最常用的illumina芯片測序（Microarray-based methods），還有RRBS一類的我不知道怎么翻譯成中文……有興趣的人自己查吧。
我最常用的芯片測序大概覆蓋了450,000或者850,000個位點，看你選那種芯片（illumina公司設計的，全球一個標準，誰都改不了……）。也是目前科研屆的最主流工具，有點有很多啦——比如便宜，其他不重要。測序完的結(jié)果是一種叫做IDAT文件的東西，這就是我程序分析的起點。這樣一個文件包含的內(nèi)容，就是我剛才說的：* 一群人中每一個的體內(nèi)的450,000或者850,000CpG的甲基化程度。*，科學家要做的，一般來說就是：
1：把它們分開，誰是癌癥，誰健康？
2：整理數(shù)據(jù)，Normalization一類的跟上。
3：開始對兩批人的幾百萬CpG位點進行比較。
4：看比較出來的CpG位點哪些有差異？他們對應到什么疾病或者通路中。(本體論oncology)
OK，簡單來說就是這么個流程。實際上很復雜的，如果不復雜，這領域也就沒有存在的可能性。
上述是ChAMP包的主要函數(shù)和其組成的研究流程。
其中藍色部分代表的是有關于數(shù)據(jù)的準備部分，就是說你要先確認數(shù)據(jù)沒問題，如果有問題，你必須停下來，回頭找什么問題，如果處理錯了重新處理，程序錯了該程序，如果實驗做錯了……你這批幾萬美元的數(shù)據(jù)就直接可以報廢了。
其中紅色部分代表的是數(shù)據(jù)分析部分，就是說確認數(shù)據(jù)沒問題之后，你可以用它產(chǎn)出結(jié)果了，這些函數(shù)里都封裝著嚴謹?shù)慕y(tǒng)計算法，涉及多重線性回歸（multiple linear regression）、網(wǎng)絡分析（network analysis）、細胞中各類分解（Cell Type Detection）等等。這一些函數(shù)大概集中了差不多我領域的所有分析角度了。這就是絕大多數(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)絡的差異化分析。（看到木有，整個生信研究，就是各種找兩群人間的差異……但是真心很不好做。）
首先通過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)存儲在了包中，可以用下面代碼讀?。?/p>
data(testDataSet)
比如只對于Shiny或者Plotly感興趣的同學，可以直接用上面代碼讀取數(shù)據(jù)，調(diào)到文章后部的GUI程序部分，開始把玩看效果。
3.1 首先是數(shù)據(jù)Load
下面的testDir指向的是我提供的450K測試數(shù)據(jù)的地址，champ.load()函數(shù)可以直接從目錄讀取數(shù)據(jù)，你只需要把illumina機器跑出來的IDAT文件全部解壓放在目錄下就行了，但是目錄下還必須放一個有關于樣本信息的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的分布。這個小程序是我一開始學習Shiny和Plotly的時候測試代碼，一直留下來了。這也是所有GUI里最簡單的一個，最適合用于分析代碼結(jié)構，效果如下：
顯示的東西都是，你數(shù)據(jù)的CpG位點是如何覆蓋的，比如第一張圖（左上）是所有CpG在染色體上的分布，等等……代碼太長，而且不易讀（不是我的錯，Shiny這個框架寫東西縮進空格太多，結(jié)構類似于HTML那種。） 大家可以看Github。
3.2 下一步做的是Quality Control
代碼如下：
QC.GUI(beta=myLoad$beta,arraytype="450K")
然后會出現(xiàn)一系列的圖像：
會有一個網(wǎng)頁自動打開，其中有幾個tab，每一個tab代表了一些分析圖像，用戶可以點擊那些tab，每一次轉(zhuǎn)到一個圖像上，需要等待幾秒鐘，后臺在計算，計算完畢后，圖像會自動生成。圖像用的是plotly，可交互，鼠標放上去會有顯示。
這一個函數(shù)的功能是告訴用戶，
tab1: 樣本之間的Clustering情況是怎樣的？
tab2: 測序芯片上的兩種測序技術（Type-I和Type-II）之間差異如何？
tab3: 所有樣本的分布式怎樣的？
tab4: 根據(jù)當前數(shù)據(jù)進行了層聚類樹狀圖是怎樣的。
tab5: 這批數(shù)據(jù)中差異性最大的一批CpG的熱力圖。
3.3 然后做的是Normalization
Normalization就是為了把上面的tab2里的兩種測序方法得到的數(shù)據(jù)不一致現(xiàn)象消除，否則如果你通過分析函數(shù)知道CpG有差異，你如何知道這個差異是疾病帶來的？還是測序方法差異帶來的？這一部至關重要，現(xiàn)在已經(jīng)有很多方法用于完成這個問題，還有科學家在樂此不疲地開發(fā)更多算法完成這一步。
myNorm <- champ.norm(beta=myLoad$beta,arraytype="450K",cores=5)
規(guī)則化之后的結(jié)果是顯而易見的，type2-BMIQ是規(guī)則化之后的分布曲線，和紅色的已經(jīng)一致了，所以測序方法上的差異，不會導致結(jié)果出現(xiàn)問題。
3.4 然后做的是SVD
這個分析講真很簡單，但是確實很有用，算法一句話形容就是：把數(shù)據(jù)進行SVD分解，查看Component有沒有與某些重要的Factor相關。SVD大家都知道，主成分分析第一步就是SVD，其性質(zhì)類似于把重要的Component給找出來了。而與Factor的相關，可以告訴使用者，你的數(shù)據(jù)中的大部分差異，是又你想研究的因素（比如癌癥、疾病、年齡等等）導致的，還是由你不想研究的東西，比如樣本的種族、數(shù)據(jù)產(chǎn)出的研究所等等一些無關痛癢的東西決定的。如果是后者，那很不幸你的數(shù)據(jù)質(zhì)量也不高……
champ.SVD(beta=myNorm,pd=myLoad$pd)
比如上述我這批測試數(shù)據(jù)的結(jié)果就非常好：因為數(shù)據(jù)中的差異，顯著與樣本的PhenoType互相關聯(lián)。換言之，你如果做針對Sample_Group（就是Cancer和Normal）的分析，得到的結(jié)果就是由于疾病導致的。
額外提一句，如果數(shù)據(jù)與其他無關的東西相關，而與你想要研究的東西不相關，那么你就需要想很多辦法了，比如，如果上述的圖片中，紅色那個板塊（顯著相關的意思，越紅越相關）是在Array那一行，而不是在Sample_Group那一行，那就意味著你的數(shù)據(jù)中的差異更多是Array導致的，而Array是芯片中的樣本排列方式，換言之就是……你的實驗不夠好。但有些問題很難避免，比如你的100個病人年齡不同很正常，老人身體上自然會有一些弱化，年輕人可能更好，樣本的年齡很自然地會引入一些誤差。
上述這種情況一旦遇到是非常痛苦的，至今也有很多我領域的人在研究如何把這樣的錯誤修正，如何矯正數(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ā)生了變化，而這些變化又是如何導致了基因發(fā)生了變化，最終導致了人體生病。而做的方法直觀上簡單的可怕：
你有100個癌癥病人，100個正常人，每個人身體中都有450K個CpG的位點在測序出來，針對其中的每一個CpG，你都有200個數(shù)據(jù)對不對？如果這一個CpG在100個正常人中和100個癌癥病人中的甲基化水平都差不多，你還會繼續(xù)懷疑它嗎，當然不會！
但如果，你的100個癌癥病人普遍在這一個CpG上的甲基化水平高（不太嚴謹?shù)呛苄蜗蟮卣f，就是DNA那個CpG的序列外層越來越多的部分被甲基附集上去了），而那100個普通人的甲基化水平不高，那這個位點就很有嫌疑了對吧。
為什么需要一定量的樣本，比如100個？因為如果你只找兩個人，一個德國人一個中國人，那個德國人高，那個中國人矮，你能因此就說德國人在人種上比中國人高嗎？當然不能……但如果你找到100個德國人，在找到100個中國人，比較以后，就比較可信了，這涉及t.test()里的power的問題。
這就是做研究的艱難：你得找到100個病人，讓他們同意治療，然后耗資幾萬幾十萬完成測序，有了數(shù)據(jù)還要祈禱實驗員沒有點錯試管，數(shù)據(jù)下來了如果由于年齡、人種等原因，數(shù)據(jù)差異已經(jīng)找不到了，你還得想辦法修正這些問題，然后你可以開始比較，從幾百萬位點（基因、SNP、CpG都一樣）中找出那些可能有關系的……等你做完這一切，可能一兩年已經(jīng)過去了，而你本科畢業(yè)就進入IT界的同學可能已經(jīng)工資兩萬多了，這還絕對算是科研中很快節(jié)奏的項目了，所以我覺得社會真的應該給科研工作者更多尊重，做最難的活，拿最低的工資。
代碼如下：
myDMP <- champ.DMP(beta = myNorm,pheno=myLoad$pd$Sample_Group)
我的程序函數(shù)之間的接口是設計過的，就是說，你從上一步跑出來的程序，可以自動被下一步識別到（當然你要是改名字什么的那就找不到了，只能自己在參數(shù)中設定）。完成這一步以后，就找出了幾百萬位點中可能有問題的。做這一步我最討厭的就是：有時候一把跑出來，幾百萬位點中，幾萬甚至十幾萬都是顯著的（做年齡的時候就這樣，因為年齡對于人身體的影響太大，可謂是全方位無死角的）！
然后是整個包中的最復雜的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)一個稍微復雜一點的網(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把它關掉。
下一張是整個R包實現(xiàn)起來最復雜的一個頁面：
用戶可以直接輸入基因的名字來選擇你們想要查詢的基因，看這些基因里是不是收到了”癌變CpG“的影響。比如，在圖中的左側(cè)輸入NFIX，點擊Submit，右側(cè)需要計算十秒鐘左右，然后上邊會顯示出這個基因中顯著差異CpG。下面會顯示出基因的不同區(qū)域，比如是基因開始，還是在基因后端一類的。我例子中這個基因，有大約十多個CpG位點，而且這些位點無一例外有顯著差異，所以基本可以確定這個基因與癌癥關系有關系。（但是具體關系是什么很難確定，誰導致了誰（因果）都不一定，這就必須要去實驗室點試管做實驗才能知道了。），在下方，我用R語言爬網(wǎng)站的技術，抓去了這個基因在一個比較權威的數(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針對基因的富集情況做的排序圖。（否則人體中有兩萬多基因，用戶怎么知道應該在上一個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)這兩個概念是之前的科學家提出來的，就是因為找單一的CpG已經(jīng)不夠顯著了！比如一個基因跨距幾千位點，通過上一步分析，就一個CpG有問題，你算他有問題嗎？這就是一個比較尷尬的問題……所以有的科學家就覺得，單一CpG影響可能不可靠（感覺上有點像木村資生的中性理論，覺得有些位點為人太隨意了，而有些真正受重大影響的基因，其實他們的基因序列上會受到一系列的暴擊傷害—— 一連串的CpG都會出現(xiàn)很明顯的差異，既然這樣，我們就專注于研究那些比較要緊的基因吧，有些暫時看不出來的就不管了，比如那些零零散散的差異CpG就當它不存在吧。
所以，比如先用代碼：
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ù)不斷都比較長的差異片段，科學家們覺得，這樣的連續(xù)差異片段，對于基因的影響會更加明顯，只找這樣的片段，可以使得計算生物學的打擊精度更為準確，也可以讓最終找出來的結(jié)論數(shù)據(jù)更少，便于實驗人員篩選。
而其實找出來以后，差別就不大了，你依然需要把這一些區(qū)域Mapping到基因組上，找到對應的基因，然后該做什么繼續(xù)做，又回到普通的比較模式了。這個網(wǎng)頁的其他兩個tab自動完成了匹配，這就是一些很簡單的匹配工作而已，不過畢竟省了一些時間。
另外一個類似的東西就是DMB，那個東西出現(xiàn)的原因是，有的科學家覺得，DMR這樣的區(qū)域還不夠顯著，DNA上的甲基化出現(xiàn)變化，可能是綿延幾千位點的！而且只會在基因以外的區(qū)域，但是這些基因以外的區(qū)域發(fā)生變化，卻會導致基因的表達發(fā)生變化。你可以想象成，北京周邊的河北在大煉鋼鐵，然后北京也跟著霧霾了，大概就是這意思。
我個人對這種假設還是有點懷疑的……根據(jù)這個提出這個理論的作者的文章，他通過找到這種Block之后，將它與一些癌癥形狀進行了correlation，結(jié)果是顯著的，然后這種理論大概就有了，相應的算法也出現(xiàn)了，雖然并沒有見到這樣的算法確實找到了某些重要的生物突破，但是身為這方面的人，我也不得不實現(xiàn)了他這個算法，放在這里供用戶使用。
我個人覺得correlation出來的結(jié)果，距離因果還是有很大距離的。建立在Correlation上的假設理論，甚至算法，實在是有點懷疑它的統(tǒng)計power。
我感覺得做計算生物很大一個問題就是不要統(tǒng)計過度，比如我根據(jù)相關性得出一些結(jié)論，我在用這批結(jié)論去做另一個聚類或者相關，在得到一個統(tǒng)計值結(jié)論。這樣的power減弱的很快啊。你做一次統(tǒng)計，錯誤率可能只是20%，那么根據(jù)其結(jié)果再錯一次統(tǒng)計，錯誤率就會上升到 1-80%*80%=36%！除非在一次統(tǒng)計過后，用相對嚴格的標準去篩選統(tǒng)計值，也許有助于這種情況好轉(zhuǎn)。
Differential Methylation Block的代碼和圖都不展示了，程序的Vignette里邊有。
3.7 尋找作用通路網(wǎng)絡中的疾病關聯(lián)小網(wǎng)絡
這貨根本就沒有中文譯名，上邊是我自己翻譯的。簡而言之就是，人體內(nèi)的作用網(wǎng)絡實在是太多了，幾百幾千吧，每一個都涉及了一系列基因。這就是過往的科學家研究出來的，比如某個通路會導致頭疼，然后有幾百個蛋白（及其背后的基因）都是通過共同作用導致了這一場頭疼的。但是如果你頭疼，不見得是所有這些基因都出了問題，而是可能其中的某一部分，甚至于只有一兩個出了問題。所以你就可以基于已經(jīng)存在的那個網(wǎng)絡，再集合數(shù)據(jù)，找出哪些網(wǎng)絡可能出問題了？或者是這些大網(wǎng)絡中的哪些基因具體除了問題。
這個功能很重要，否則做完上邊幾步，用戶只會知道那些基因有問題，至于他們之間有沒有關系，是不是會同時作用于某些網(wǎng)絡，就沒法知道的了。但其實也不復雜，只需要自己寫個程序，匹配一樣基因和網(wǎng)絡就完了，只不過數(shù)據(jù)的準備啊，洗啊，匹配啊，也是夠煩的，所以這個函數(shù)就提供了全套的分析。
myEpiMod <- champ.EpiMod(beta=myNorm,pheno=myLoad$pd$Sample_Group)
上述代碼之后可以輸出一些圖片到當前目錄的文件夾下：
其中的顏色代表了基因的甲基化上下調(diào)的趨勢。（人們得了某個病，不一定是導致那個病的基因高表達了，也有可能是抑制那種病的某個基因低表達了，甚至于可能是抑制那個抑制生病的基因的基因突然高表達了，這可以無窮下去……）。所以一個網(wǎng)絡了，不見得都是同樣的表達趨勢，藍點代表甲基化高表達，而黃點代表低表達。一言以蔽之，可能是這個網(wǎng)絡中的這一些基因”同謀“導致了你在研究的這個疾病發(fā)生了。
3.8 拷貝數(shù)變異
這個的意思就是可能你的測序數(shù)據(jù)里，表現(xiàn)出某些疾病導致的原因，是有些基因片段被復制的此處過多或者過少。檢測方法很聰明，基因人們是看不見的，但是測序的原理，是用一些片段去把與CpG有關的片段”粘“下來。然后對這些片段進行測序。一般來說，我是不關心粘下來多少的，因為我研究的基因的甲基化水平，其實就是某一個為點的被甲基化和未被甲基化水平的比例。例如，假設正常的情況下，某個CpG所在片段會被”粘“下來一萬次，其中8000片上的CpG都沒有被甲基化，而2000片被甲基化了。那么最后到我手里的數(shù)據(jù)，就是這個CpG的甲基化值為2000/8000=0.25。如果被粘下來兩萬次，那么必然是16000片上沒有被甲基化，而4000片上被甲基化了，我得到的比例還是0.25。完全感覺不到區(qū)別??！
所以檢測辦法，就是去看，到底多少片被拉下來了，然后與正常人作比較。這個函數(shù)不是我寫的，是以前一個教授，但是我覺得有點粗糙，因為它只能統(tǒng)計出每一條染色體上的拷貝數(shù)變化差異，我覺得這不夠的。如果有機會，我希望更新一下程序，讓它打擊更為定點，就是找出來是哪些CpG發(fā)生了突變導致了拷貝數(shù)變異，然后與其他的指標，比如表達值一類的作比較，這樣應該會有新發(fā)現(xiàn)。
myCNA <- champ.CNA(intensity=myLoad$intensity,pheno=myLoad$pd$Sample_Group)
代碼過后，只會生成一些這樣的圖片：
圖的大概意思就是，C這個疾病情況下，相比于正常人，在全基因上那些位置發(fā)生拷貝數(shù)變異。我從來沒有用過這個功能，感覺精度根本不夠，我就努力確保了它程序沒問題，運轉(zhuǎn)順利。
3.9 細胞種類分離問題
這個堪稱是我這個研究方向上目前全球在最努力攻關的問題……我做它快三年了，也完全沒有實質(zhì)性突破和進展，好在全球范圍內(nèi)好像也沒有人做出來。簡而言之是問題是這樣的，人體內(nèi)充斥著各種細胞，你給病人采血取樣，其實取到的也是細胞的組合。然后問題就來了，你確定某個疾病是所有的細胞都出問題了嗎？還是某幾種細胞出了問題？但是你的數(shù)據(jù)是細胞們混合在一起的，你怎么分開？比如下面的公式：
X=C.F
這里的X
是你的身體中采集到了CpG甲基化值，或者基因表達值，或者各種亂七八糟的數(shù)值（只要不是基因序列本身，應該都存在這個問題吧，基因序列確實穩(wěn)定，不會發(fā)生變化，這也就是為什么SNP測序可以發(fā)展成為產(chǎn)業(yè)，而其他的不可能。）。這個X其實是又兩個部分組成的：C和F。C是每一個位點（或者基因）在單一細胞中的表達值，比如你的血液里有紅細胞，白細胞和血小板，C就是每一個位點在這些細胞中的表達值矩陣，F(xiàn)
是你采集的樣本中，各種細胞的比例。比如你的血液中紅細胞，白細胞和血小板的比例。
這帶來的問題就是。每個人體內(nèi)的細胞比例是不同的。甚至于，有些疾病本身就是因為細胞比例發(fā)生變化導致的，血小板太少會導致流血不止對吧？進而，導致疾病發(fā)生的基因也是在不同細胞中的，如果他們混在一起，你怎么知道，是誰具體導致了癌癥？其性質(zhì)就像是之前在SVD的時候我提過的干擾因素，你找了兩群人想找癌癥和正常人的區(qū)別，結(jié)果一群人是老人，一群人是小孩，這怎么可能找的對？年齡還可以修正，只要你知道大家的年齡就可以修正。如果你找的100個樣本血細胞混合是不同的，你怎么繼續(xù)做研究，而且你也不知道大家的細胞比例是什么，連修正的辦法都沒有。
所以就是X
是你的數(shù)據(jù)，C是純化過的單一細胞內(nèi)的甲基化值，F(xiàn)是樣本的細胞混合比例。X是C和F的矩陣乘。目的就是為了識別出F
。
這個問題還要分為兩個方向，是關于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中有細胞比例。這個感覺有道理，但是應該對數(shù)據(jù)質(zhì)量要求很高，如果你還有其他可能的沒有消除干凈的變量，就會導致你的Component其實根本不是細胞比例……然后就錯掉了。
他的函數(shù)我也集成進來了，但是不見得很好用：
myRefFree <- champ.reffree(beta=myNorm,pheno=myLoad$pd$Sample_Group)
1
大概就是覺得自己可以嘗試開一個技術博客，記錄一下平時用到的技術。這個項目是去年下半年做的一個小東西，不是什么正式項目，科研屆不會在乎什么用戶可視化的無聊軟件的，花時間也不長。發(fā)在CSDA上主要是覺得Shiny和Plotly挺好用的，也許可以給對R語言有興趣的同學一些啟發(fā)和靈感。我學的也不好，尤其有些算法在深度和概念上理解可能有錯，歡迎同學們指正。
最后總結(jié)幾句：計算生物學真的不好做，對統(tǒng)計編程要求挺高的，而且很多技術可以與很多領域?qū)?。比如最新方興未艾的那個”Data Science“。
搞這行的人，至少接觸過幾個G以上的數(shù)據(jù)，這些數(shù)據(jù)說太大不至于，但是能把它們處理清楚，意味著他處理百萬行以下的超小數(shù)據(jù)集是沒問題的。
做生信統(tǒng)計嚴謹要求很高，每一個同學都必須搞清楚什么時候用Wilcox，什么時候用t.test()，都能知道correlation不等于因果，所以不會犯里約奧運會時候那種幼稚的錯誤：”里約奧運會是社交網(wǎng)絡上被談論的最多的一屆！“（恐怕不是奧運會進步了，而是社交網(wǎng)絡發(fā)展了）。
有些IT技術挺好的，但是一直進不到科研屆來，我曾經(jīng)需要跑一個數(shù)據(jù)，用R不可能跑完，用Python估計兩三天，最后用了Spark，3分鐘出結(jié)果……那個時候我就覺得，領域之間應該互相多學習學習。金融IT就業(yè)單位，不要看到掛著生物兩個字的人，就推到門外，這其實可能是世界上最接近”大數(shù)據(jù)“和”數(shù)據(jù)科學“的專業(yè)了。
生物數(shù)據(jù)和金融數(shù)據(jù)差別能有多大呢？如果細胞識別算法能出來，是不是一只股票的投資維度分析也會發(fā)生變化？通信問題中的音頻混合問題是不是可以被解決？Shiny和Plotly這個技術用于數(shù)據(jù)展示不也很方便嗎？SVD分解針對各種因素的Correlation與社會科學家分析城市省份經(jīng)濟發(fā)展的方式應該是一致的吧？生物數(shù)據(jù)有干擾因素的性質(zhì)，與海洋大氣采集的數(shù)據(jù)干擾有沒有共通性？問題往往是相通的，技術更經(jīng)常是一樣。
