轉(zhuǎn)錄組學(xué)習(xí)八（功能富集分析）

轉(zhuǎn)錄組學(xué)習(xí)一（軟件安裝）
轉(zhuǎn)錄組學(xué)習(xí)二（數(shù)據(jù)下載）
轉(zhuǎn)錄組學(xué)習(xí)三（數(shù)據(jù)質(zhì)控）
轉(zhuǎn)錄組學(xué)習(xí)四（參考基因組及gtf注釋探究）
轉(zhuǎn)錄組學(xué)習(xí)五（reads的比對(duì)與samtools排序）
轉(zhuǎn)錄組學(xué)習(xí)六（reads計(jì)數(shù)與標(biāo)準(zhǔn)化）
轉(zhuǎn)錄組學(xué)習(xí)七（差異基因分析）
轉(zhuǎn)錄組學(xué)習(xí)八（功能富集分析）

任務(wù)

  • 選擇p<0.05而且abs(log2FC)大于1的基因?yàn)轱@著差異表達(dá)基因集，對(duì)這個(gè)基因集用R包做KEGG/GO超幾何分布檢驗(yàn)分析。
  • 把表達(dá)矩陣和分組信息分別作出cls和gct文件，導(dǎo)入到GSEA軟件分析。

<font color=orange>GO富集分析</font>

一、bioconductor注釋數(shù)據(jù)庫的探究

JIMMY_數(shù)據(jù)包library(org.Hs.eg.db)簡(jiǎn)介

  • 官網(wǎng)上已有19個(gè)注釋數(shù)據(jù)庫，網(wǎng)址bioconductor.org/packages/release/ 看了一下植物只有一個(gè)擬南芥的。
  • 如果不在注釋數(shù)據(jù)庫里，使用AnnotationHub來構(gòu)建自己的Org.db。
library(AnnotationHub)
hub <- AnnotationHub()
# 可以用query（）函數(shù)來查找你要的物種注釋信息
# 選擇的格式是OrgDb.
query（hub, "Solanum lycopersicum"）
sl <- hub[["AH55774"]]
  • 對(duì)ORG.Mm.eg.db進(jìn)行簡(jiǎn)單的探究
library(org.Mm.eg.db)
keytypes(org.Mm.eg.db) ##查看有哪些數(shù)據(jù)類型，包含著各大主流數(shù)據(jù)庫的數(shù)據(jù)。
###用select函數(shù)，就可以把任意公共數(shù)據(jù)庫的數(shù)據(jù)進(jìn)行一一對(duì)應(yīng)。
### keys是原始的ID，columns是轉(zhuǎn)換之后的ID，keytype是要指定的原始ID類型
select(org.Mm.eg.db,keys = "ENSMUSG00000031762",columns = c("SYMBOL","GENENAME","UNIGENE","REFSEQ"),keytype = "ENSEMBL")
  • 對(duì)ENSEMBL-ID進(jìn)行轉(zhuǎn)換
diff_gene_DESeq_raw <- subset(res_deseq, padj<0.05 & (log2FoldChange > 1 | log2FoldChange < -1))
diff_gene_DESeq_name <- row.names(diff_gene_DESeq_raw)
diff_gene_DESeq_transID<- select(org.Mm.eg.db, keys= diff_gene_DESeq_name, columns= c("SYMBOL", "GENENAME", "UNIGENE", "REFSEQ"), keystype="ENSEMBL")
image
  • 結(jié)果可知幾大數(shù)據(jù)庫的基因?qū)?yīng)關(guān)系，其中有部分的數(shù)據(jù)是未知的？還不知道是什么原因，查找了下原始的小鼠gtf文件，發(fā)現(xiàn)注釋的基因是最近的？所以說可能是數(shù)據(jù)還不全導(dǎo)致的吧。

二、基因GO分析：

利用clusterProfiler的R包進(jìn)行GO分析

enrichGO(gene, OrgDb, keytype = "ENTREZID", ont = "MF",
  pvalueCutoff = 0.05, pAdjustMethod = "BH", universe, qvalueCutoff = 0.2,
  minGSSize = 10, maxGSSize = 500, readable = FALSE, pool = FALSE)
  • 基本參數(shù)：gene:即基因ID；OrgDb:物種注釋數(shù)據(jù)；keytype:ID的類型 ont 三個(gè)層面來闡述基因功能，生物學(xué)過程（BP），細(xì)胞組分（CC），分子功能（MF）；pAdjustMethod:P值校正方法
  • 進(jìn)行GO分析的基本代碼
ego <- enrichGO(gene = row.names(diff_gene_deseq2),  OrgDb = org.Mm.eg.db, keytype = "ENSEMBL", ont = "MF")

###氣泡圖
dotplot(ego, font)

### 網(wǎng)絡(luò)圖
enrichMap(ego, vertex.label.cex=1.2, layout=igraph::layout.kamada.kawai)

###GO圖
plotGOgraph(ego)
image

image

image

三、基因KEGG分析：

  • KEGG介紹：KEGG數(shù)據(jù)庫（Kyoto Encyclopedia of Genes and Genomes）京都基因與基因組百科全書：KEGG PATHWAY數(shù)據(jù)庫由一系列經(jīng)人手繪制而成KEGG代謝路徑圖的構(gòu)成，以代表關(guān)于代謝以及其他細(xì)胞與生物機(jī)能的實(shí)驗(yàn)成果。每一張路徑圖俱包括了一個(gè)分子互動(dòng)與反應(yīng)的網(wǎng)絡(luò)圖，為連系基因組內(nèi)之基因及路徑所示過程生成的基因產(chǎn)物（以蛋白質(zhì)為主）而設(shè)。行使相同功能的基因聚在一起

  • 目前KEGG等數(shù)據(jù)庫，收錄的是已有的研究結(jié)果。但這些pathway的信息，遠(yuǎn)沒有到達(dá)完善的水準(zhǔn)。大部分通路只是了解1個(gè)大概的調(diào)控途徑，而中間有什么轉(zhuǎn)錄因子參與、是否還有其他代謝物的生成，都是不知道的。這些通路的完整性，也會(huì)影響pathway富集分析結(jié)果。例如，基因A發(fā)生變化了，看起來下游基因沒有變化。也許是還有其他的調(diào)控在起作用，只是這些調(diào)控作用現(xiàn)在還不知道而已。

  • pathway 和 GO富集分析結(jié)果的解讀，應(yīng)該從生物學(xué)意義的角度出發(fā)，P value 和 Q value只是個(gè)參考而已，那些不顯著的通路也值得解讀（從功能注釋的角度解讀，而不是從富集分析的角度解讀）。只要結(jié)果可以解釋，有意義，不用太迷信P value。

  • 縮寫查看官網(wǎng)：http://www.genome.jp/kegg/catalog/org_list.html

  • KEGG基本過程

diff_gene_deseq2_transID_kegg <- diff_gene_deseq2_transID[,4]
ekegg <- enrichKEGG(diff_gene_deseq2_transID_kegg,keyType = "kegg",organism = "mmu",pvalueCutoff = 0.05,pAdjustMethod = "BH",qvalueCutoff = 0.1)
### 畫氣泡圖：
dotplot(ekegg,font.size=8)

### 顯示通路圖
browseKEGG(ekegg,'mmu01100')
image

image

<font color=orange>基因集富集分析GSEA</font>

參考文章GSEA分析是個(gè)什么鬼？（上）, GSEA分析是個(gè)什么鬼？（下）。文章將GSEA分析做了詳細(xì)的揭示，目前僅對(duì)看懂的部分做記錄，知道有這個(gè)分析方法，以后有需要再做詳細(xì)學(xué)習(xí)吧

一、基因富集分析概念

  • 基因富集分析(Gene set enrichment analysis, GSEA)的基本思想是使用預(yù)定義的基因集（通常是來自于功能注釋或先前實(shí)驗(yàn)的結(jié)果），將基因按照在兩類樣本中的差異表達(dá)程度排序，然后檢驗(yàn)預(yù)先設(shè)定的基因集合是否在這個(gè)排序表的頂端或者底端富集?；蚣系母患治鍪菣z測(cè)基因集合而不是單個(gè)基因的表達(dá)變化，因此可以包含這些細(xì)微的表達(dá)變化，預(yù)期得到更為理想的結(jié)果。
    • 分析的基因集合而不是單個(gè)基因
    • 將基因與預(yù)定義的基因集進(jìn)行比較
    • 富集分析

二、與通常富集方法GO和KEGG的比較：

  • GO和 KEGG 側(cè)重于比較兩組間的基因表達(dá)差異，集中關(guān)注少數(shù)幾個(gè)顯著上調(diào)或下調(diào)的基因，這容易遺漏部分差異表達(dá)不顯著卻有重要生物學(xué)意義的基因，忽略一些基因的生物特性、基因調(diào)控網(wǎng)絡(luò)之間的關(guān)系及基因功能和意義等有價(jià)值的信息。

  • GSEA分析不需要指定明確的差異基因閾值，算法會(huì)根據(jù)實(shí)際數(shù)據(jù)的整體趨勢(shì)， 為研究者們提供了一種合理地解決目前芯片分析瓶頸問題的方法,即使在沒有先驗(yàn)經(jīng)驗(yàn)存在的情況下也能在表達(dá)譜整體層次上對(duì)數(shù)條基因進(jìn)行分析,從而從數(shù)理統(tǒng)計(jì)上把表達(dá)譜芯片數(shù)據(jù)與生物學(xué)意義很好地銜接起來,使得研究者們能夠更輕松、更合理地解讀芯片結(jié)果。

三、通常做差異分析設(shè)定閾值與后續(xù)KEGG與GO分析的問題**：

  • 差異基因列表后，都會(huì)在此之上提供給客戶Pathway 以及GO 富集分析。這里這些篩選出來的差異基因本身就會(huì)存在一些問題。

  • 在實(shí)際應(yīng)用于生物學(xué)高通量數(shù)據(jù)時(shí)，對(duì)于差異基因檢出的閾值，異常的敏感，客戶需要給出差異基因的一個(gè)明確的定義(閾值)，例如abs(FC) ≧2.0 & p ≦ 0.05。這種一刀切的閾值，對(duì)于發(fā)現(xiàn)真正的生物學(xué)效應(yīng)，許多時(shí)候是一種障礙，因?yàn)?strong>實(shí)際通過芯片觀測(cè)到的RNA 表達(dá)變化，往往是層層的負(fù)反饋調(diào)控后的結(jié)果，并且不同組織對(duì)于表達(dá)差異的敏感度是不同的：在神經(jīng)遞質(zhì)系統(tǒng)內(nèi)，一個(gè)1.2 倍的表達(dá)差異即可能產(chǎn)生及其顯著的效應(yīng)。

四、GSEA富集過程的基本步驟

  1. 計(jì)算富集分?jǐn)?shù)（Enrichment Score）
  2. 估計(jì)富集分?jǐn)?shù)的顯著性水平
  3. 矯正多重假設(shè)檢驗(yàn)
    image

    image

五、基本GSEA分析過程

GSEA_genelist <- diff_gene_deseq2_raw$log2FoldChange ### 對(duì)diff_gene的結(jié)果進(jìn)行分析
names(GSEA_genelist)<- rownames(diff_gene_deseq2_raw) ### 設(shè)置名字
GSEA_genelist<- sort(GSEA_genelist,decreasing = TRUE) ### 排序


gsem_gene <- gseGO(geneList = GSEA_genelist,OrgDb = org.Mm.eg.db, keyType = "ENSEMBL", ont = "MF")

gseaplot(gsem_gene,geneSetID = "GO:0000977")

image

看不懂在說啥，以后再慢慢研究這一類的圖吧。

總結(jié)：從2017年10月7日~12月4號(hào)，軟件安裝——數(shù)據(jù)下載——原始數(shù)據(jù)的質(zhì)控——參考基因組與注釋GTF文件的探究——READS比對(duì)，排序——計(jì)數(shù)，標(biāo)準(zhǔn)化——差異基因分析——功能富集分析。Linux基本操作、shell腳本編寫、Perl腳本編寫、軟件參數(shù)的具體含義了解。兩個(gè)月的時(shí)間終于跟著大神們的步伐將轉(zhuǎn)錄組的流程給學(xué)習(xí)了一遍。途中遇到了不少艱難的事，也花了不少深夜與周末的時(shí)間。好在終于把一個(gè)個(gè)知識(shí)點(diǎn)給攻克，一章章的任務(wù)分析給堅(jiān)持了下來。仍然有許多待學(xué)習(xí)的地方，比如一些軟件的進(jìn)階參數(shù)的選擇，結(jié)果的更加準(zhǔn)確解讀，R語言的語法、各種可視化圖片的繪制以及如何解讀，還有越來越覺得重要的一個(gè)大坑：統(tǒng)計(jì)學(xué)背景知識(shí)。這些都將是后續(xù)學(xué)習(xí)的方向。這不是結(jié)束，以后還會(huì)根據(jù)學(xué)習(xí)到的各種新知識(shí)來更好的完善這個(gè)分析流程大框架。加油吧。

參考文章

  1. 【基因富集分析_學(xué)習(xí)筆記】https://mp.weixin.qq.com/s?__biz=MzIwNTEwMTUyOQ==&mid=2649693906&idx=1&sn=341682dad10a9b52f3290239042c30f5&chksm=8f2dbe64b85a3772d3bc439498560ec22638783cb321be71e7ebb383a4de0186b3ea9b384475&scene=21#wechat_redirect
  2. 【PANDA姐的轉(zhuǎn)錄組入門(8)：差異基因結(jié)果注釋】https://mp.weixin.qq.com/s?__biz=MzIwNTEwMTUyOQ==&mid=2649694917&idx=1&sn=a318f8cf98f306d46963986011c73600&chksm=8f2d8273b85a0b65270a986f7bb28e8efa7fd15309505a5a026bbedaafaff7e30e4d4f13dd02&scene=21#wechat_redirect
  3. 【(偽）從零開始學(xué)轉(zhuǎn)錄組（8）：富集分析】https://mp.weixin.qq.com/s?__biz=MzI1MjU5MjMzNA==&mid=2247484528&idx=1&sn=3af297d4163a70b049f861fc28c85bc2&chksm=e9e02dd1de97a4c79ab04a2c012ee6fe4b0aae5e1df59bdf2bf6b19e6ac1cb08ce379a6694e8&scene=21#wechat_redirect
  4. 【差異基因結(jié)果注釋】http://www.itdecent.cn/p/4910d7cec5c8
  5. 【GSEA分析是個(gè)什么鬼？（上）】https://mp.weixin.qq.com/s?__biz=MzAwMzY4MTYxNw==&mid=2655753566&idx=2&sn=5b5b2c93a7618a69da2cbc6638f03da0&chksm=80884960b7ffc076af53ae74caadb5dbb25d240c31660792e8727964d0177d6a17af7ca5fc5c&mpshare=1&scene=1&srcid=0816ADpKId3sPzgbYfubrFCf#rd
  6. 【GSEA是個(gè)什么鬼?(下)】https://mp.weixin.qq.com/s?__biz=MzAwMzY4MTYxNw==&mid=2655754973&idx=1&sn=3b87d5cb8ddd2d5d77e413e9a87342da&chksm=808846e3b7ffcff5a6b41985b707f52170f20eabe15fc43264b3d14a3ccf4100263789eab856&mpshare=1&scene=1&srcid=0816gHxusewlJeILw0fWxgi3#rd
  7. 【RNA-seq結(jié)果圖片如何解讀？（第一彈）】http://mp.weixin.qq.com/s/OFuP7nGGM3V9ghZ6lI1QuA
  8. 【RNA-seq中GO、KEGG結(jié)果圖如何解讀】http://mp.weixin.qq.com/s/UowQnL4bD7QUFHIXQ_JQKQ
最后編輯于 ：
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請(qǐng)聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成，瀏覽時(shí)請(qǐng)結(jié)合常識(shí)與多方信息審慎甄別。
平臺(tái)聲明：文章內(nèi)容（如有圖片或視頻亦包括在內(nèi)）由作者上傳并發(fā)布，文章內(nèi)容僅代表作者本人觀點(diǎn)，簡(jiǎn)書系信息發(fā)布平臺(tái)，僅提供信息存儲(chǔ)服務(wù)。

相關(guān)閱讀更多精彩內(nèi)容

友情鏈接更多精彩內(nèi)容