r语言进行go富集分析_R语言GEO数据挖掘-功能富集分析

功能富集分析

在得到了差异基因的基础之上，进一步进行功能富集分析，这里我们使用clusterprofiler包

本文将对差异基因进行 GO, KEGG注释并完成可视化,GSEA分析

Sys.setlocale('LC_ALL','C')

library(tidyverse)

library(ggplot2)

library(dplyr)

library(clusterProfiler)

library(org.Hs.eg.db)

load(file = "DEG_all.Rdata")

head(DEG)

##定义筛选logFC>=1, P<=0.05的差异基因

DEG_2% rownames_to_column("genename") #%>% as_tibble() %>%

filter(abs(logFC)>=0.58,P.Value<=0.05)

genelist

## ID转换

genelist% bitr(fromType = "SYMBOL",

toType = c("ENSEMBL", "ENTREZID"),

OrgDb = org.Hs.eg.db)

head(genelist)

dim(genelist)

GO富集分析

注意clusterprofiler中使用的富集ID是ENTREZID，如果不是需要进行转换，也很简单

我们这里其实可以用我们自带的probe中的ID

使用了semi_join这个筛选连接

注意这里我故意使用的是筛选连接，其实inner_join, 等等都可以，我就是要用这些不熟悉的

但是这里的probe其实基因名是有重复的，我们还要进行去重的操作

如果自己有probe信息

if(F){

load("probe.Rdata")

names(probe)[2]

DEG_2

probe

##设定差异阈值筛选

DEG_2% arrange(abs(logFC)) %>%

filter(abs(logFC)>=0.58,P.Value<=0.05)

dim(DEG_2)

gene%as_tibble() %>% dplyr::semi_join(DEG_2,by="genename") %>% .$ENTREZ_GENE_ID

head(gene)

length(gene)##相当于是取了交集，用DEG_2中的基因来做筛选，得到了12549个基因

}

我们假设没有probe注释信息

ego

OrgDb = org.Hs.eg.db,

keyType = 'ENSEMBL',##指定gene id的类型

ont = "all",##GO分类

pAdjustMethod = "BH",

pvalueCutoff = 0.01,##设置阈值

qvalueCutoff = 0.05,

readable=TRUE)

write.csv(ego,file = "ego.csv")

ego[1:5,1:5]

ONTOLOGY ID

GO:0019882 BP GO:0019882

GO:0048002 BP GO:0048002

GO:0019884 BP GO:0019884

GO:0002478 BP GO:0002478

GO:0070482 BP GO:0070482

Description

GO:0019882 antigen processing and presentation

GO:0048002 antigen processing and presentation of peptide antigen

GO:0019884 antigen processing and presentation of exogenous antigen

GO:0002478 antigen processing and presentation of exogenous peptide antigen

GO:0070482 response to oxygen levels

GeneRatio BgRatio

GO:0019882 353/12683 385/19623

GO:0048002 320/12683 347/19623

GO:0019884 314/12683 341/19623

GO:0002478 307/12683 334/19623

GO:0070482 372/12683 420/19623

############

class(ego)

str(ego)

## 三个结果整合在一起

go_dot

facet_grid(ONTOLOGY~.)

go_dot

## 单个绘制

ego_MF

OrgDb = org.Hs.eg.db,

keyType = 'ENSEMBL',##指定gene id的类型

ont = "MF",##GO分类

pAdjustMethod = "BH",

pvalueCutoff = 0.01,##设置阈值

qvalueCutoff = 0.05,

readable=TRUE)

dotplot(ego_MF,showCategory=6)

image.png

KEGG富集分析

代码的风格都是类似的

organism = 'hsa',

pvalueCutoff = 0.05

)

head(kk)

## 打开网页浏览通路

browseKEGG(kk,"hsa04151")

## 点图

dotplot(kk)

## 柱状图

barplot(kk)

点图

dotplot(kk)

image.png

柱状图

barplot(kk)

image.png

Pathview包用于通路可视化

先调整数据

library(pathview)

pathview(genelist$ENTREZID,pathway.id ="04151",species = "hsa")

DEG_2% rownames_to_column("genename")

DEG_2% rename(SYMBOL='genename') %>%

inner_join(DEG_2,by='genename') %>%

## select配合everything排序把，改变变量顺序

dplyr::select(ENTREZID,logFC,everything())

head(DEG_2)

###

ENTREZID logFC genename ENSEMBL AveExpr t P.Value

1 218 -3.227263 ALDH3A1 ENSG00000108602 10.302323 -10.710306 4.482850e-05

2 1545 -3.033684 CYP1B1 ENSG00000138061 13.287607 -10.505888 4.995753e-05

3 1543 -9.003353 CYP1A1 ENSG00000140465 11.481268 -8.371476 1.762905e-04

4 11148 -1.550587 HHLA2 ENSG00000114455 6.595658 -7.443431 3.337672e-04

5 8140 -2.470333 SLC7A5 ENSG00000103257 13.628775 -7.298868 3.708688e-04

6 25976 -1.581274 TIPARP ENSG00000163659 12.764218 -7.024252 4.552834e-04

adj.P.Val B

1 0.3134585 -4.048355

2 0.3134585 -4.049713

3 0.7374232 -4.069681

4 0.9308066 -4.083411

5 0.9308066 -4.085966

6 0.9522252 -4.091192

###

构建geneList对象

这是比较关键的一步，需要构建geneList对象，至少需要包括ID与排序好的数值型变量

数值变量降序排列

sign.pos参数包括的位置有4种："bottomleft", "bottomright", "topleft" and "topright".

geneList

names(geneList)

geneList

#geneList

head(geneList)##构建geneList对象成功

1580 7025 54474 9247 56603 4907 ## geneid

2.378155 2.365342 2.258682 2.180382 2.143749 2.030780 ## logFC

## 这样logFC值就能反应在通路的变化中了

pathview(geneList,pathway.id ="04151",species = "hsa",same.layer = F )

image.png

调整函数的参数，改变图形

比较重要的是kegg.native = F，可以调整输出为矢量图

## 看不出变化的参数kegg.native，默认设置为T

pathview(geneList,pathway.id ="04151",species = "hsa",

out.suffix = "KEGG",##输出文件后缀

kegg.native = T,

same.layer = F )

## 可以产生矢量图pdf文档，

pathview(geneList,pathway.id ="04151",species = "hsa",

out.suffix = "KEGG",##输出文件后缀

kegg.native = F,

split.group = T,#是否分开节点组

sign.pos= "topleft",##图注的位置

same.layer = F )

GSEA分析

GSEA分析在高分文章中还是一个很常见的分析，之前认为官方软件的作图达不到很高的要求，现在终于可以解决了clusterprofiler包更新的画图功能模仿官网的图，颜值更高，操作更简单，大家都值得拥有。

KEGG的GSEA富集分析

head(geneList)

1580 7025 54474 9247 56603 4907

2.378155 2.365342 2.258682 2.180382 2.143749 2.030780

#geneList

kk2

organism = 'hsa',

nPerm = 1000,

minGSSize = 120,

pvalueCutoff = 0.5,

verbose = FALSE)

head(kk2)

ID Description setSize enrichmentScore

hsa04080 hsa04080 Neuroactive ligand-receptor interaction 293 -0.3593278

hsa04024 hsa04024 cAMP signaling pathway 185 -0.3656357

hsa04060 hsa04060 Cytokine-cytokine receptor interaction 244 -0.3358228

hsa04020 hsa04020 Calcium signaling pathway 172 -0.3566872

hsa04062 hsa04062 Chemokine signaling pathway 165 0.3112935

hsa05202 hsa05202 Transcriptional misregulation in cancer 163 -0.3537420

NES pvalue p.adjust qvalues rank

hsa04080 -1.481227 0.002649007 0.1562914 0.1422098 1912

hsa04024 -1.426245 0.007062147 0.2083333 0.1895629 1697

hsa04060 -1.356610 0.012345679 0.2427984 0.2209227 2055

hsa04020 -1.378406 0.017118402 0.2524964 0.2297470 1748

hsa04062 1.323383 0.025806452 0.2676695 0.2435530 1332

hsa05202 -1.361718 0.027220630 0.2676695 0.2435530 2463

leading_edge

hsa04080 tags=31%, list=16%, signal=27%

hsa04024 tags=26%, list=14%, signal=23%

hsa04060 tags=30%, list=17%, signal=25%

hsa04020 tags=26%, list=15%, signal=22%

hsa04062 tags=19%, list=11%, signal=17%

hsa05202 tags=32%, list=21%, signal=26%

GSEA图可视化

之前的图确实颜值不够

gseaplot(kk2,geneSetID = "hsa04080")

不够

新功能gseaplot2-模仿GSEA软件的图-确实更好看了

library(enrichplot)

gseaplot2(kk2,##KEGG gse的对象

geneSetID = "hsa04080",

pvalue_table=T)

image.png

修改线条颜色

现在的图-高端大气上档次

gseaplot2(kk2,##KEGG gse的对象

geneSetID = "hsa04080",

color="red",

pvalue_table=T)

MgsiDB的GSEA富集分析

gmtfile

gene

egmt

head(egmt)

ID Description GeneRatio BgRatio

CELL_FRACTION CELL_FRACTION CELL_FRACTION 451/4461 493/5270

MEMBRANE_FRACTION MEMBRANE_FRACTION MEMBRANE_FRACTION 309/4461 339/5270

pvalue p.adjust qvalue

CELL_FRACTION 1.764189e-06 0.0003775364 0.00036398

MEMBRANE_FRACTION 1.870439e-04 0.0200136995 0.01929506

#write.csv(egmt,file = "egmt.csv")

GSEA

egmt2

pvalueCutoff = 0.5,#

verbose=FALSE)

head(egmt2)

gseaplot2(egmt2, geneSetID = 1, title = egmt2$Description[1],pvalue_table=T)

image.png

还可以同时画多个

gseaplot2(egmt2, geneSetID = 1:3,pvalue_table=T)

image.png

Gene-Concept Network

这里的输入egmt2可以是ego,kk,edo都可以，进行可视化

edox

cnetplot(edox, foldChange=geneList)

image.png

换个输出形状

cnetplot(edox, foldChange=geneList, circular = TRUE, colorEdge = TRUE)

image.png

r语言进行go富集分析_R语言GEO数据挖掘-功能富集分析相关推荐

r语言进行go富集分析_R语言实现GO分析
链客,专为开发者而生,有问必答! 此文章来自区块链技术社区,未经允许拒绝转载. 我们上一期介绍了如何实现GO分析的可视化,运行了GOplot包自带的数据并且很畅通.然而我们如何才能获取那些可以直接输入 ...
灰色关联分析_R语言使用灰色关联分析(Grey Relation Analysis,GRA)中国经济社会发展指标...
原文链接: http://tecdat.cn/?p=16881tecdat.cn 灰色关联分析包括两个重要功能. 第一项功能:灰色关联度,与correlation系数相似,如果要评估某些单位,在使用 ...
go语言第三方包安装方法_R语言3.6.3 安装程序下载及破解方法
下载地址百度网盘链接: https://pan.baidu.com/s/16smT3ceIjqaupn54AdgmgQ 提取码:7hap 解压密码:关注[菜瓜程序猿]微信公众号,回复[解压密码]获取 ...
GEO数据挖掘全流程分析
声明:以下学习资料根据"生信技能树"网络系列免费教学材料整理而成,代码来自"生信技能树"校长jimmy的github.GEO数据库挖掘系列知识分享课程,于201 ...
r语言进行go富集分析_R语言-GO富集分析的超几何检验和可视化
Gene Ontology 可分为分子功能(Molecular Function),生物过程(biological process)和细胞组成(cellular component)三个部分.蛋白质或 ...
r语言进行go富集分析_R语言：clusterProfiler进行GO富集分析和Gene_ID转换
一.读取文件,ID转换 1.读取文件 library(clusterProfiler) library(org.Hs.eg.db) #读取文件,原始文件中使用空格分割的 go_ythdf2 go_yt ...
r语言pls分析_R语言：生存分析
生存分析处理预测特定事件将要发生的时间.它也被称为故障时间分析或分析死亡时间.例如,预测患有癌症的人将存活的天数或预测机械系统将失败的时间. 命名为survival的R语言包用于进行生存分析.此包包含 ...
r语言pls分析_R语言中的偏最小二乘PLS回归算法
偏最小二乘回归: 我将围绕结构方程建模(SEM)技术进行一些咨询,以解决独特的业务问题.我们试图识别客户对各种产品的偏好,传统的回归是不够的,因为数据集的高度分量以及变量的多重共线性.PLS是处理这些 ...
r语言pls分析_R语言中的偏最小二乘回归PLS-DA
主成分回归(PCR)的方法本质上是使用第一个方法的普通最小二乘(OLS)拟合来自预测变量的主成分(PC).这带来许多优点: 预测变量的数量实际上没有限制. 相关的预测变量不会破坏回归拟合. 但是, ...

r语言进行go富集分析_R语言GEO数据挖掘-功能富集分析

r语言进行go富集分析_R语言GEO数据挖掘-功能富集分析相关推荐

最新文章

热门文章