GO,KEGG,DO富集分析

这是我听B站鲮鱼不会飞视频(GO,KEGG,DO富集分析)里的笔记哦~
当然有的地方加上的是我自己的理解,如果哪里和视频上讲的不一样那是我自己发挥的,自行忽略....

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

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

  • tophat2是把reads回帖到基因组上;
  • Cufflinks在计算基因表达量;
  • Cuffdiff比较control和treatment找差异基因(生成一个数据框)

后面的富集分析,一般只做GO分析,KEGG pathway 分析,最多再做一个DO分析,公司一般用的是已经成熟的database,这就导致数据分析不完全,而且公司用的数据库很多时候都已经过时了,所以我们需要自己学会做下游的富集分析。


GO分析的理论知识

what is Gene Ontology(GO)? 基因"本体论"
基因本体论是对基因在不同维度和不同层次上的描述。

对基因的描述一般从三个层面进行:

  • Cellular component,CC 细胞成分
  • Biological process, BP 生物学过程
  • Molecular function,MF 分子功能

这三个层面具体是指:

  • Cellular component解释的是基因存在在哪里,在细胞质还是在细胞核?如果存在细胞质那在哪个细胞器上?如果是在线粒体中那是存在线粒体膜上还是在线粒体的基质当中?这些信息都叫Cellular component。
  • Biological process是在说明该基因参与了哪些生物学过程,比如,它参与了rRNA的加工或参与了DNA的复制,这些信息都叫Biological process
  • Molecular function在讲该基因在分子层面的功能是什么?它是催化什么反应的?
    So, we will have a gene annotation infarmation.
    立足于这三个方面,我们将得到基因的注释信息。

得到GO注释

model organism ---> annotated database
non-model organism ---> search database or blast

  • 模式生物 ---> 有标准的注释数据库;
  • 非模式生物 ---> 自己搜注释数据库(怎们搜后面具体介绍),搜不到就用blast的办法解决。

做GO分析的思路:

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

也就是RNA-Seq先测出各组的基因表达分布:
control gene expression distribution
treatment gene expression distribution
control VS treatment ---> DEG : differential expression genes
通过比较 control 和 treatment 得到差异表达基因
再去做GO富集分析:
DEG ---> GO enrichment analysis
用找到的差异基因去做GO富集分析,希望能从这三方面找到和我们背景不一样的地方。
比如,在疾病研究的时候,进行药物治疗之后某些基因的表达量明显的发生了变化,拿这些基因去做GO分析发现在Biological process过程当中集中在RNA修饰上,然后在此基础上继续进行挖掘。这个例子就是想启示大家拿到差异表达基因DEG只是一个开始,接下来就应该去做GO注释,之后需要进行一个分析看这些注释主要集中在哪个地方。假如我们有100个差异表达基因其中有99个都集中在细胞核里,那我们通过GO分析就得到了一个显著的分布。

GO富集分析原理:有一个term注释了100个差异表达基因参与了哪个过程,注释完之后(模式生物都有现成的注释包,不用我们自己注释),计算相对于背景它是否显著集中在某条通路、某一个细胞学定位、某一种生物学功能。

KEGG enrichment analysis?
把生物体中所有的pathway都要进行富集分析
DO enrichment analysis?
看目标基因是否在某个疾病或某一类疾病当中富集

代码部分

  1. RNA-seq分析中第一步是:fastq ---> bam (tophat2 , hisat2 , star....);
  2. 使用 cufflink 输入文件是bam;
  3. 使用 cutffdiff 做差异表达分析,输入文件是 bam GTF注释文件

这套流程是上游分析,拿到cutffdiff结果之后就可以转到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)

得到差异表达基因,赋值给一个新的数据框 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

网页工具做GO分析

david

打开谷歌 ---> 搜david --->第一个点进去 ---> 就是做GO分析的最常用的网站

做GO分析的最常用的网站-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)

这时当前文件夹就多了一个名为sign_gene_id.txt里面装有所有gene_id 的txt文件。

所有差异表达基因的数据框

Enrichr

还推荐了一个常用的网站 Enrichr

R代码做GO分析

用R可以做一些网站上不能做的东西。

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") #这个包里存有人的注释文件

# 载入包dian
library(clusterProfiler)
library(topGO)
library(Rgraphviz)
library(pathview)
library(org.Hs.eg.db)

2.作图前处理——提取symbol ID --> 转换为ENTREZID

DEG.gene_symbol = as.character(output.gene_id$gene_id) #获得基因 symbol ID

防止在做GO分析的时候出现报错,需要将symbolID转换成ENTREZID:用mapIds函数就可以转换ID。

DEG.entrez_id = mapIds(x = org.Hs.eg.db,
                       keys = DEG.gene_symbol,
                       keytype = "SYMBOL",
                       column = "ENTREZID")

这时就已经把symbolID转换成ENTREZID了,但会出现个别的转换不成功的情况,就是图中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,学习下如何用代码

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是日本主导的一个项目对gene和genome进行了非常详细的注释

KEGG网页分析里面有非常全的注释。

KEGG pathway 分析和上面介绍的GO分析是一样的只是把enrichGO()函数改成 enrichKEGG()
GO分析:enrichGO() —---—> KEGG pathway分析:enrichKEGG()
所以不细讲啦~

DO分析介绍

DO分析用enrichDO()函数,是做疾病的,这里我们做一下:

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

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

其实这个问题的核心是非模式生物怎样找到org.db数据库(标准注释库)?因为有了注释库后面的分析都一样一样的~
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里面找到番茄相关的database即org.Solanum_lycopersicum.eg.sqlite 注:Solanum_lycopersicum是番茄的拉丁名和它对应的编号AH59087

query(hub, "Solanum")  # Solanum番茄的拉丁名

找到之后把它下载下来:

Solanum.OrgDb <- hub[["AH59087"]]

此时,番茄的database就会赋值到变量Solanum.OrgDb
解决完标准注释库的问题,剩下的和模式生物做富集分析完全一样了~~

最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 151,829评论 1 331
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 64,603评论 1 273
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 101,846评论 0 226
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 42,600评论 0 191
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 50,780评论 3 272
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 39,695评论 1 192
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 31,136评论 2 293
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 29,862评论 0 182
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 33,453评论 0 229
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 29,942评论 2 233
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 31,347评论 1 242
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 27,790评论 2 236
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 32,293评论 3 221
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 25,839评论 0 8
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 26,448评论 0 181
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 34,564评论 2 249
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 34,623评论 2 249

推荐阅读更多精彩内容