3Y叔的clusterProfiler-book阅读Chapter 3 Universal enrichment analysis

Chapter 3 Universal enrichment analysis

ClusterProfiler支持多种基因GO富集分析和和超几何检验，但是不能用于不支持的物种以及slim版本的以及新颖的功能注释(如GO通过BlastGO和KEGG通过KAAS来分析其数据)。

ClusterProfiler为超几何分布提供给了enricher函数，为基因富集提供了GEA分析功能，这些功能用户可以进行自定义。他们接受两个附加参数TERM2GENE 和 TERM2NAME. 前者是第一列具有Term ID 和第二列具有映射的基因名GENE的一个数据框；TERM2NAME是一个第一列具有Term ID 和第二列具有相应的Term的名字，后者是可选的。

3.1文件的输入

为了进行过量表达分析，需要的数据集的形式

1. 包含基因ID的向量(Vector),这些基因可以通过差异获取(差异分析可以通过Deseq2包进行分析)

为了基因富集分析(Gsea Set Enricher Analysis)分析，需要对基因列表进行排序，DOSE提供了一个示例数据集geneList，这个数据集来自于R包中的数据集breastCancerMAINZ,其中包含200个样本，其中1级包含29个样本，2级包含136个样本和3级包含35个样本。这个geneList是作者通过计算2级样本与3级样本几何平均值的对数比。

1. geneList 包含的三个元素

数值向量

基因名向量

排序向量(包数字按照降序排列)

形如

自己准备的文件的数据的提示

包含两列，一列是基因的ID(不能包含有重复ID)，另一列是是基因的差异值(通常是Log2FC)

利用命令行进行数据的处理

d <- read.csv(your_csv_file)
## assume that 1st column is ID
## 2nd column is fold change

## feature 1: numeric vector
geneList <- d[,2]

## feature 2: named vector
names(geneList) <- as.character(d[,1])

## feature 3: decreasing order
geneList <- sort(geneList, decreasing = TRUE)

将样本装在进入R可以使用DOSE()函数，命令如下

data(geneList,package="DOSE")

head(geneList)

选择基因的变异倍数大于2的基因用于下一步的分析

gene=names(geneList)[abs(geneList)>2]

head(gene)

3.2维基通路分析(WikiPathways analysis）

WikiPathways 是一个不断更新的通路的数据库，可以在data.wikipathways.org网站上下载相关的通路数据，主要是gmt数据，它里面含有多个物种 ，可以根据自己的需要进行下载，2020年一月份的人的gmt文件如下

选择合适的gmt文件后，然后产生TERM2GENE 和 TERM2NAME这两个文件进行后续的enricher分析和GSEA功能分析。

library(magrittr)

library(clusterProfiler)

data(geneList, package="DOSE")

gene <- names(geneList)[abs(geneList) > 2]

wpgmtfile <- system.file("extdata/wikipathways-20180810-gmt-Homo_sapiens.gmt", package="clusterProfiler")

wp2gene <- read.gmt(wpgmtfile)

wp2gene <- wp2gene %>% tidyr::separate(ont, c("name","version","wpid","org"), "%") ##将下载的未分割的文本进行分割，按照%进行分割，并分别进行命名如c()里面的内容，   %>%主要是将后面的结果传递到前面

wpid2gene <- wp2gene %>% dplyr::select(wpid, gene) #TERM2GENE选择wpid和gene作为TERM2GENE

wpid2name <- wp2gene %>% dplyr::select(wpid, name) #TERM2NAME 选择wpid和name作为TERM2NAME

ewp <- enricher(gene, TERM2GENE = wpid2gene, TERM2NAME = wpid2name)

head(ewp)

 

 

ewp2 <- GSEA(geneList, TERM2GENE = wpid2gene, TERM2NAME = wpid2name, verbose=FALSE)

head(ewp2)

得到的结果是包含有ID ，Description，GeneRatio ，BgRatio， pvalue ，p.adjust，qvalue，GeneID的列的数据

如下

 

 

 

 

通过以下代码将gene ENTREZID转化成GeneSymbol

library(org.Hs.eg.db)
ewp <- setReadable(ewp, org.Hs.eg.db, keyType = "ENTREZID")
ewp2 <- setReadable(ewp2, org.Hs.eg.db, keyType = "ENTREZID")
head(ewp)

3.3Cellular marker(细胞标记物)

cell_markers <- vroom::vroom('http://bio-bigdata.hrbmu.edu.cn/CellMarker/download/Human_cell_markers.txt') %>% tidyr::unite("cellMarker", tissueType, cancerType, cellName, sep=", ")%>%                                                                              dplyr::select(cellMarker, geneID) %>%                                                                                                                                               dplyr::mutate(geneID = strsplit(geneID, ', ')) 

 cell_markers

y <- enricher(gene, TERM2GENE=cell_markers, minGSSize=1)
DT::datatable(as.data.frame(y))

 

3.43.4 MSigDb analysis

Users can download GMT files from Broad Institute and use read.gmt to parse the file to be used in enricher() and GSEA().

用户可以使用GMT文件，利用read,gmt来读取文件，用于后续的enricher()和GSEA()函数分析

R package, msigdbr, that already packed the MSigDB gene sets in tidy data format that can be used directly with clusterProfiler

 

.