有一個term注釋了100個差異表達基因參與了哪個過程，注釋完之后（模式生物都有現(xiàn)成的注釋包，不用我們自己注釋），計算相對于背景它是否顯著集中在某條通路、某一個細胞學(xué)定位、某一種生物學(xué)功能。
黃晶_id122019.05.09 09:12:57字數(shù) 2,030閱讀 30,597
這是我聽B站
鯪魚不會飛視頻（GO，KEGG，DO富集分析）里的筆記哦~
當(dāng)然有的地方加上的是我自己的理解，如果哪里和視頻上講的不一樣那是我自己發(fā)揮的，自行忽略....
市面上公司做RNA-Seq的一般流程是：
tophat2 ---> Cufflinks ---> Cuffdiff ---> R
tophat2是把reads回帖到基因組上；
Cufflinks在計算基因表達量；
Cuffdiff比較control和treatment找差異基因（生成一個數(shù)據(jù)框）
后面的富集分析，一般只做GO分析，KEGG pathway 分析，最多再做一個DO分析，公司一般用的是已經(jīng)成熟的database，這就導(dǎo)致數(shù)據(jù)分析不完全，而且公司用的數(shù)據(jù)庫很多時候都已經(jīng)過時了，所以我們需要自己學(xué)會做下游的富集分析。
GO分析的理論知識
what is Gene Ontology（GO）? 基因"本體論"
基因本體論是對基因在不同維度和不同層次上的描述。
對基因的描述一般從三個層面進行：
Cellular component，CC 細胞成分
Biological process， BP 生物學(xué)過程
Molecular function，MF 分子功能
這三個層面具體是指：
Cellular component解釋的是基因存在在哪里，在細胞質(zhì)還是在細胞核？如果存在細胞質(zhì)那在哪個細胞器上？如果是在線粒體中那是存在線粒體膜上還是在線粒體的基質(zhì)當(dāng)中？這些信息都叫Cellular component。
Biological process是在說明該基因參與了哪些生物學(xué)過程，比如，它參與了rRNA的加工或參與了DNA的復(fù)制，這些信息都叫Biological process
Molecular function在講該基因在分子層面的功能是什么？它是催化什么反應(yīng)的？
So, we will have a gene annotation infarmation.
立足于這三個方面，我們將得到基因的注釋信息。
得到GO注釋
model organism ---> annotated database
non-model organism ---> search database or blast
模式生物 ---> 有標準的注釋數(shù)據(jù)庫；
非模式生物 ---> 自己搜注釋數(shù)據(jù)庫(怎們搜后面具體介紹)，搜不到就用blast的辦法解決。
做GO分析的思路：
control VS treatment ---> DEG ---> GO enrichment analysis
也就是RNA-Seq先測出各組的基因表達分布：
control gene distribution
treatment gene distribution
control VS treatment ---> DEG : differential genes
通過比較 control 和 treatment 得到差異表達基因
再去做GO富集分析：
DEG ---> GO enrichment analysis
用找到的差異基因去做GO富集分析，希望能從這三方面找到和我們背景不一樣的地方。
比如，在疾病研究的時候，進行藥物治療之后某些基因的表達量明顯的發(fā)生了變化，拿這些基因去做GO分析發(fā)現(xiàn)在Biological process過程當(dāng)中集中在RNA修飾上，然后在此基礎(chǔ)上繼續(xù)進行挖掘。這個例子就是想啟示大家拿到差異表達基因DEG只是一個開始，接下來就應(yīng)該去做GO注釋，之后需要進行一個分析看這些注釋主要集中在哪個地方。假如我們有100個差異表達基因其中有99個都集中在細胞核里，那我們通過GO分析就得到了一個顯著的分布。
GO富集分析原理：有一個term注釋了100個差異表達基因參與了哪個過程，注釋完之后（模式生物都有現(xiàn)成的注釋包，不用我們自己注釋），計算相對于背景它是否顯著集中在某條通路、某一個細胞學(xué)定位、某一種生物學(xué)功能。
KEGG enrichment analysis？
把生物體中所有的pathway都要進行富集分析
DO enrichment analysis？
看目標基因是否在某個疾病或某一類疾病當(dāng)中富集
代碼部分
RNA-seq分析中第一步是：fastq ---> bam (tophat2 , hisat2 , star....)；
使用 cufflink 輸入文件是bam；
使用 cutffdiff 做差異表達分析，輸入文件是 bam GTF注釋文件
這套流程是上游分析，拿到cutffdiff結(jié)果之后就可以轉(zhuǎn)到R里進行下一步分析：
1. load cutffdiff result
cuffdiff_result = read.table(file="./hela_gene_exp.diff",header = T,sep = "\t")cuffdiff_result$sample_1 = "treat"cuffdiff_result$sample_2 = "ctrl"
2. select DEG
I. FPKM1 or FPKM2 > 1
II. log2(fold change) >1 or < -1
III. p-value < 0.05
select_vector = (cuffdiff_result$value_1 > 1 | cuffdiff_result$value_2 > 1) & (abs(cuffdiff_result$log2.fold_change.) >= 1) & (cuffdiff_result$p_value < 0.05)
得到差異表達基因，賦值給一個新的數(shù)據(jù)框 cuffdiff_result.sign
cuffdiff_result.sign = cuffdiff_result[select_vector,]> dim(cuffdiff_result.sign)[1] 2739 14
這就說明在這種條件下我們篩選出了2739個差異表達基因，這個其實有點多了。我們只是為了走一下富集分析的流程，所以把條件再加緊一下，篩出一千來個基因去做分析正好：
> select_vector = (cuffdiff_result$value_1 > 1 | cuffdiff_result$value_2 > 1) & (abs(cuffdiff_result$log2.fold_change.) >= 1.5) & (cuffdiff_result$p_value < 0.05)> cuffdiff_result.sign = cuffdiff_result[select_vector,]> dim(cuffdiff_result.sign)[1] 1268 14
網(wǎng)頁工具做GO分析
david
打開谷歌 ---> 搜david --->第一個點進去 ---> 就是做GO分析的最常用的網(wǎng)站
做GO分析的最常用的網(wǎng)站-david
如圖，點進去后，把gene list 放進白色的框里
寫個代碼把 gene id 那一列單獨提取出來并保存到本地。output.gene_id = data.frame(gene_id = cuffdiff_result.sign$gene_id)write.table(output.gene_id,file="./sign_gene_id.txt",col.names = F,row.names = F,sep = "\t",quote = F)
這時當(dāng)前文件夾就多了一個名為sign_gene_id.txt里面裝有所有g(shù)ene_id 的txt文件。
所有差異表達基因的數(shù)據(jù)框
Enrichr
還推薦了一個常用的網(wǎng)站 Enrichr
R代碼做GO分析
用R可以做一些網(wǎng)站上不能做的東西。
1.準備工作——安裝R包
# 安裝包source("https://bioconductor.org/biocLite.R")BiocManager::install("clusterProfiler") #用來做富集分析BiocManager::install("topGO") #畫GO圖用的BiocManager::install("Rgraphviz")BiocManager::install("pathview") #看KEGG pathway的BiocManager::install("org.Hs.eg.db") #這個包里存有人的注釋文件# 載入包dianlibrary(clusterProfiler)library(topGO)library(Rgraphviz)library(pathview)library(org.Hs.eg.db)
2.作圖前處理——提取symbol ID --> 轉(zhuǎn)換為ENTREZID
DEG.gene_symbol = as.character(output.gene_id$gene_id) #獲得基因 symbol ID
防止在做GO分析的時候出現(xiàn)報錯，需要將symbolID轉(zhuǎn)換成ENTREZID：用mapIds函數(shù)就可以轉(zhuǎn)換ID。
DEG.entrez_id = mapIds(x = org.Hs.eg.db, keys = DEG.gene_symbol, keytype = "SYMBOL", column = "ENTREZID")
這時就已經(jīng)把symbolID轉(zhuǎn)換成ENTREZID了，但會出現(xiàn)個別的轉(zhuǎn)換不成功的情況，就是圖中NA的地方，我們進行以下操作即可去掉：
DEG.entrez_id = na.omit(DEG.entrez_id)
做好準備工作，我們就開始做富集分析
3.GO分析代碼
BP（Biological process）層面上的富集分析：
erich.go.BP = enrichGO(gene = DEG.entrez_id, OrgDb = org.Hs.eg.db, keyType = "ENTREZID", ont = "BP", pvalueCutoff = 0.5, qvalueCutoff = 0.5)##分析完成后，作圖dotplot(erich.go.BP)
解讀BP層面富集分析圖：
橫坐標是GeneRatio，意思是說輸入進去的基因，它每個term（縱坐標）站整體基因的百分之多少。圓圈的大小代表基因的多少，圖中給出了最大的圓圈代表60個基因，圓圈的顏色代表P-value，也就是說P-value越小gene count圈越大，這事就越可信。
dotplot(erich.go.BP)
CC（Cellular component）層面上的富集分析：
erich.go.CC = enrichGO(gene = DEG.entrez_id, OrgDb = org.Hs.eg.db, keyType = "ENTREZID", ont = "CC", pvalueCutoff = 0.5, qvalueCutoff = 0.5)## 畫圖barplot(erich.go.CC)
barplot(erich.go.CC)
一般GO分析畫這兩個圖就可以了，有時也把GO分析畫成樹形圖，可以更加幫助我們理解。
plotGOgraph(erich.go.BP)
plotGOgraph(erich.go.BP)
樹狀圖很大，所以我們用代碼把它存成pdf，學(xué)習(xí)下如何用代碼
pdf(file="./enrich.go.bp.tree.pdf",width = 10,height = 15)plotGOgraph(erich.go.BP)dev.off()
至此，GO分析就做完了 ----> over
KEGG pathway介紹
KEGG: Kyoto Encyclopedia of Genes and Genomes
KEGG是日本主導(dǎo)的一個項目對gene和genome進行了非常詳細的注釋
KEGG網(wǎng)頁分析里面有非常全的注釋。
KEGG pathway 分析和上面介紹的GO分析是一樣的只是把enrichGO()函數(shù)改成 enrichKEGG()
GO分析：enrichGO() —---—> KEGG pathway分析：enrichKEGG()
所以不細講啦~
DO分析介紹
DO分析用enrichDO()函數(shù)，是做疾病的，這里我們做一下：
enrichDO(gene = DEG.entrez_id,ont = "DO",pvalueCutoff = 0.5,qvalueCutoff = 0.5)
非模式生物如何做富集分析
其實這個問題的核心是非模式生物怎樣找到org.db數(shù)據(jù)庫（標準注釋庫）？因為有了注釋庫后面的分析都一樣一樣的~
search org.db ----> 套路分析
非模式生物但有參考基因組的情況
以番茄為例，番茄有參考基因組但不在標準注釋庫里
先安裝兩個包
source("https://bioconductor.org/biocLite.R")BiocManager::install("AnnotationHub")BiocManager::install("biomaRt")# 載入包library(AnnotationHub)library(biomaRt)
自己制作一個OrgDB
hub <- AnnotationHub::AnnotationHub()
使用query在我們制作的OrgDB --> hub里面找到番茄相關(guān)的database即org.Solanum_lycopersicum.eg.sqlite 注：Solanum_lycopersicum是番茄的拉丁名和它對應(yīng)的編號AH59087
query(hub, "Solanum") # Solanum番茄的拉丁名
找到之后把它下載下來：
Solanum.OrgDb <- hub[["AH59087"]]
此時，番茄的database就會賦值到變量Solanum.OrgDb
解決完標準注釋庫的問題，剩下的和模式生物做富集分析完全一樣了~~
