市面上公司做RNA-Seq的一般流程是：

tophat2 ---> Cufflinks ---> Cuffdiff ---> R

tophat2是把reads回帖到基因組上；
Cufflinks在計算基因表達(dá)量；
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）? 基因"本體論"
基因本體論是對基因在不同維度和不同層次上的描述。

對基因的描述一般從三個層面進(jìn)行：

Cellular component，CC 細(xì)胞成分
Biological process， BP 生物學(xué)過程
Molecular function，MF 分子功能

這三個層面具體是指：

Cellular component解釋的是基因存在在哪里，在細(xì)胞質(zhì)還是在細(xì)胞核？如果存在細(xì)胞質(zhì)那在哪個細(xì)胞器上？如果是在線粒體中那是存在線粒體膜上還是在線粒體的基質(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

模式生物 ---> 有標(biāo)準(zhǔn)的注釋數(shù)據(jù)庫；
非模式生物 ---> 自己搜注釋數(shù)據(jù)庫(怎們搜后面具體介紹)，搜不到就用blast的辦法解決。

做GO分析的思路：

control VS treatment ---> DEG ---> GO enrichment analysis

也就是RNA-Seq先測出各組的基因表達(dá)分布：
control gene expression distribution
treatment gene expression distribution
control VS treatment ---> DEG : differential expression genes
通過比較 control 和 treatment 得到差異表達(dá)基因
再去做GO富集分析：
DEG ---> GO enrichment analysis
用找到的差異基因去做GO富集分析，希望能從這三方面找到和我們背景不一樣的地方。
比如，在疾病研究的時候，進(jìn)行藥物治療之后某些基因的表達(dá)量明顯的發(fā)生了變化，拿這些基因去做GO分析發(fā)現(xiàn)在Biological process過程當(dāng)中集中在RNA修飾上，然后在此基礎(chǔ)上繼續(xù)進(jìn)行挖掘。這個例子就是想啟示大家拿到差異表達(dá)基因DEG只是一個開始，接下來就應(yīng)該去做GO注釋，之后需要進(jìn)行一個分析看這些注釋主要集中在哪個地方。假如我們有100個差異表達(dá)基因其中有99個都集中在細(xì)胞核里，那我們通過GO分析就得到了一個顯著的分布。

GO富集分析原理：有一個term注釋了100個差異表達(dá)基因參與了哪個過程，注釋完之后（模式生物都有現(xiàn)成的注釋包，不用我們自己注釋），計算相對于背景它是否顯著集中在某條通路、某一個細(xì)胞學(xué)定位、某一種生物學(xué)功能。

KEGG enrichment analysis？
把生物體中所有的pathway都要進(jìn)行富集分析
DO enrichment analysis？
看目標(biāo)基因是否在某個疾病或某一類疾病當(dāng)中富集

代碼部分

RNA-seq分析中第一步是：fastq ---> bam (tophat2 , hisat2 , star....)；
使用 cufflink 輸入文件是bam；
使用 cutffdiff 做差異表達(dá)分析，輸入文件是 bam GTF注釋文件

這套流程是上游分析，拿到cutffdiff結(jié)果之后就可以轉(zhuǎn)到R里進(jìn)行下一步分析：

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)

得到差異表達(dá)基因，賦值給一個新的數(shù)據(jù)框 cuffdiff_result.sign

cuffdiff_result.sign = cuffdiff_result[select_vector,]
> dim(cuffdiff_result.sign)
[1] 2739   14

這就說明在這種條件下我們篩選出了2739個差異表達(dá)基因，這個其實有點多了。我們只是為了走一下富集分析的流程，所以把條件再加緊一下，篩出一千來個基因去做分析正好：

> 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 --->第一個點進(jìn)去 ---> 就是做GO分析的最常用的網(wǎng)站

image

如圖，點進(jìn)去后，把gene list 放進(jìn)白色的框里

image

寫個代碼把 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文件。

image

Enrichr

還推薦了一個常用的網(wǎng)站 Enrichr

R代碼做GO分析

用R可以做一些網(wǎng)站上不能做的東西。

1.準(zhǔn)備工作——安裝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") #這個包里存有人的注釋文件

# 載入包dian
library(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的地方，我們進(jìn)行以下操作即可去掉：

DEG.entrez_id = na.omit(DEG.entrez_id)

image

做好準(zhǔn)備工作，我們就開始做富集分析

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層面富集分析圖：
橫坐標(biāo)是GeneRatio，意思是說輸入進(jìn)去的基因，它每個term（縱坐標(biāo)）站整體基因的百分之多少。圓圈的大小代表基因的多少，圖中給出了最大的圓圈代表60個基因，圓圈的顏色代表P-value，也就是說P-value越小gene count圈越大，這事就越可信。

image

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)

image

一般GO分析畫這兩個圖就可以了，有時也把GO分析畫成樹形圖，可以更加幫助我們理解。

plotGOgraph(erich.go.BP)

image

樹狀圖很大，所以我們用代碼把它存成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進(jìn)行了非常詳細(xì)的注釋

KEGG網(wǎng)頁分析里面有非常全的注釋。

KEGG pathway 分析和上面介紹的GO分析是一樣的只是把enrichGO()函數(shù)改成 enrichKEGG()
GO分析：enrichGO() —---—> KEGG pathway分析：enrichKEGG()
所以不細(xì)講啦~

DO分析介紹

DO分析用enrichDO()函數(shù)，是做疾病的，這里我們做一下：

enrichDO(gene = DEG.entrez_id,ont = "DO",pvalueCutoff = 0.5,qvalueCutoff = 0.5)

非模式生物如何做富集分析

其實這個問題的核心是非模式生物怎樣找到org.db數(shù)據(jù)庫（標(biāo)準(zhǔn)注釋庫）？因為有了注釋庫后面的分析都一樣一樣的~
search org.db ----> 套路分析

非模式生物但有參考基因組的情況
以番茄為例，番茄有參考基因組但不在標(biāo)準(zhǔn)注釋庫里
先安裝兩個包

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
解決完標(biāo)準(zhǔn)注釋庫的問題，剩下的和模式生物做富集分析完全一樣了~~
GO，KEGG，DO富集分析 - 簡書 (jianshu.com)

色偷偷精品伊人,欧洲久久精品,欧美综合婷婷骚逼,国产AV主播,国产最新探花在线,九色在线视频一区,伊人大交九欧美,1769亚洲,黄色成人av

富集分析GO，KEGG，DO等

富集分析GO，KEGG，DO等

GO分析的理論知識

得到GO注釋

做GO分析的思路：

代碼部分

1. load cutffdiff result

2. select DEG

網(wǎng)頁工具做GO分析

david

Enrichr

R代碼做GO分析

1.準(zhǔn)備工作——安裝R包

2.作圖前處理——提取symbol ID --> 轉(zhuǎn)換為ENTREZID

3.GO分析代碼

KEGG pathway介紹

DO分析介紹

非模式生物如何做富集分析

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

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

色偷偷精品伊人,欧洲久久精品,欧美综合婷婷骚逼,国产AV主播,国产最新探花在线,九色在线视频一区,伊人大交九 欧美,1769亚洲,黄色成人av

富集分析GO，KEGG，DO等

GO分析的理論知識

得到GO注釋

做GO分析的思路：

代碼部分

1. load cutffdiff result

2. select DEG

網(wǎng)頁工具做GO分析

david

Enrichr

R代碼做GO分析

1.準(zhǔn)備工作——安裝R包

2.作圖前處理——提取symbol ID --> 轉(zhuǎn)換為ENTREZID

3.GO分析代碼

KEGG pathway介紹

DO分析介紹

非模式生物如何做富集分析

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

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

色偷偷精品伊人,欧洲久久精品,欧美综合婷婷骚逼,国产AV主播,国产最新探花在线,九色在线视频一区,伊人大交九欧美,1769亚洲,黄色成人av

富集分析GO，KEGG，DO等