GEO數據挖掘技術可以應用到表達芯片也可以是轉錄組測序

2019 年 10 月 31 日
筆記

GEO數據挖掘技巧，基本上該分享的都在B站和GitHub了，目錄如下：

第一講：GEO，表達芯片與R
第二講：從GEO下載數據得到表達量矩陣
第三講：對表達量矩陣用GSEA軟件做分析
第四講：根據分組信息做差異分析
第五講：對差異基因結果做GO/KEGG超幾何分佈檢驗富集分析
第六講：指定基因分組boxplot指定基因list畫熱圖
第七講：根據差異基因list獲取string數據庫的PPI網絡數據
第八講：PPI網絡數據用R或者cytoscape畫網絡圖
第九講：網絡圖的子網絡獲取
第十講：hug genes如何找

雖然一直演示的表達芯片數據分析，這些芯片分析難點主要是在ID轉換，因為不同公司設計的探針命名都不一樣，在我4年前博客整理的芯片平台對應R包找：(16)芯片探針與基因的對應關係-生信菜鳥團博客2周年精選文章集

基本上你使用我的標準數據分析代碼，下載到表達矩陣，走標準分析流程，火山圖，熱圖，GO/KEGG數據庫注釋等等，肯定可以出對應的圖表。最重要的3張圖見：你確定你的差異基因找對了嗎？

那如果是RNA-seq測序數據呢？

通常呢，RNA-seq測序數據並不會把其表達矩陣存儲在Series Matrix File(s) 裏面，所以你使用我的標準代碼：

rm(list = ls())  ## 魔幻操作，一鍵清空~  options(stringsAsFactors = F)#在調用as.data.frame的時，將stringsAsFactors設置為FALSE可以避免character類型自動轉化為factor類型  # 注意查看下載文件的大小，檢查數據  f='GSE103611_eSet.Rdata'  # https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE103611  library(GEOquery)  # 這個包需要注意兩個配置，一般來說自動化的配置是足夠的。  #Setting options('download.file.method.GEOquery'='auto')  #Setting options('GEOquery.inmemory.gpl'=FALSE)  if(!file.exists(f)){    gset <- getGEO('GSE103611', destdir=".",                   AnnotGPL = F,     ## 注釋文件                   getGPL = F)       ## 平台文件    save(gset,file=f)   ## 保存到本地  }  load('GSE103611_eSet.Rdata')  ## 載入數據  class(gset)  #查看數據類型  length(gset)  #  class(gset[[1]])  gset  # assayData: 352859 features, 48 samples

只需要把上面的GSE號替換即可，當然如果你不懂GSE號，就需要再細讀表達芯片的公共數據庫挖掘系列推文感興趣的也可以去看看；

比如對 https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE106292 上面的代碼就拿不到表達矩陣

因為，這個是RNA-seq數據，作者會把自己的表達矩陣變成Excel表格，方便大家探索！

記住，我這裡強調了是作者自己的表達矩陣，因為RNA-seq數據分析流程還不一樣！參數不一樣，軟件不一樣，數據庫不一樣，而且最後的表達矩陣的表現形式又不一樣！是原始的counts還是RPKM，TPM都不一樣！如果作者確實不上傳其表達矩陣，你也沒辦法，只能是自己走RNA-seq數據分析流程：

這裏面的知識細節太複雜了，我就不一一展開！建議大家看我們閱讀量過10萬的RNA-seq系列推文，比如：表達矩陣的歸一化和標準化，去除極端值，異常值

現在給大家一個作業

我前面舉例的 https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE106292 數據集，其實是一個WGCNA文獻，你可以看看下載到的 Excel表格如何讀入R裏面，做出作者文章的那樣的圖，可以參考關鍵問題答疑：WGCNA的輸入矩陣到底是什麼格式，詳細教程見：一文看懂WGCNA 分析(2019更新版)

這兩個圖難度非常大，基本上相當於半年作業的生信工程師經驗了，如果你能做出來，發郵件給我你的全部思考分析過程，你可以獲得我認可，畢竟相當於有了我7.6%的功力，已經是非常的了不起了！