本文主要是介绍批量富集分析气泡图的画法，对大家解决编程问题具有一定的参考价值，需要的程序猿们随着小编来一起学习吧！

目录

前言

一、compareCluster函数

二、使用步骤

1.加载包

2.读入数据

3.处理数据数据

 4.ID转换

5.通路富集

6.可视化

总结

前言

        clusterProfile是常用的基因富集分析的包，之前已经介绍过了对但样本集合进行富集分析的操作。本次我们尝试一下使用包中的compareCluster函数对多个样本进行富集分析，并使用气泡图进行展示。

提示：以下是本篇文章正文内容，下面案例可供参考

一、compareCluster函数

这个函数是clusterProfile中自带的函数，可以对多个基因集分别进行基因富集分析，并使用使用气泡图进行可视化，方便比较。

二、使用步骤

1.加载包

代码如下（示例）：

library(clusterProfiler)
library(org.Hs.eg.db)
library(stringr)
library(dplyr)
library(tibble)
library(tidyr)

2.读入数据

有两个文件，一个是数据文件，一个是病人分组

mutation <- read.csv('mutation.xls',header = T,sep = '\t')
mutation <- mutation[,c("Gene_refGene","patient","sampleID")]
> head(mutation)
  Gene_refGene patient   sampleID
1       ARID1A  S01001 S01001_EOT
2         FLT3  S01001 S01001_EOT
3          RB1  S01001 S01001_EOT
4          RB1  S01001 S01001_EOT
5         ARAF  S01001 S01001_EOT
6        ERBB2  S01001  S01001_C1
patient_information <- read.csv("panient_information.csv",header = T)
> head(patient_information)
  patient siteBestResponse SD.24wk C1有无血
1   01001               PD      NA       NA
2   01002               PR      NA        0
3   01003               SD       1       NA
4   01004               PR      NA        0
5   01006               PR      NA        0
6   01007               PR      NA       NA

3.处理数据数据

我们的表格是不符合compareCluster的要求的，所以需要处理一下

colnames(patient_information)[1] <- 'patient' ##修改指定列名
colnames(mutation)[1] <- 'gene' #修改mutation第一列列名为gene
mutation$patient <- str_replace(mutation$patient, "S","") #删除指定列中的所有S 
mutation <- mutation[,-3] #删除第三列
mutation <- distinct(mutation) #去除重复，这个函数去除的是完全相同的行
> head(mutation)
    gene patient
1 ARID1A   01001
2   FLT3   01001
3    RB1   01001
4   ARAF   01001
5  ERBB2   01001
6 PDGFRA   01001
gene <- as.data.frame(mutation[,1]) #提取第一列为dataframe格式
gene <- rownames_to_column(gene, "row") ##””内是新的列的列名
colnames(gene)[2] <- 'gene' #修改第二列列名
> head(gene)
  row   gene
1   1 ARID1A
2   2   FLT3
3   3    RB1
4   4   ARAF
5   5  ERBB2
6   6 PDGFRA
mutation <- merge(gene,mutation,by = 'gene',all = T) #根据patient列合并表格
mutation <- mutation[,c(2,3,1)] #改变表格列的顺序
> head(mutation)
  row patient gene
1  59   01027 AKT1
2  59   01010 AKT1
3  72   01027 AKT2
4  72   01023 AKT2
5  72   01022 AKT2
6  32   01006  ALK
data <- spread(mutation,patient,gene) #长变宽函数
> head(data)
  row  01001  01002 01003  01004  01006  01007 01008  01010  01012 01013  01016 01017 01020
1   1 ARID1A   <NA>  <NA> ARID1A ARID1A ARID1A  <NA> ARID1A ARID1A  <NA>   <NA>  <NA>  <NA>
2  10  BRCA2   <NA> BRCA2  BRCA2   <NA>   <NA>  <NA>   <NA>   <NA> BRCA2   <NA>  <NA>  <NA>
3  11 PIK3CA PIK3CA  <NA>   <NA>   <NA> PIK3CA  <NA> PIK3CA   <NA>  <NA> PIK3CA  <NA>  <NA>
4  12    KIT   <NA>  <NA>   <NA>   <NA>   <NA>  <NA>   <NA>   <NA>   KIT    KIT  <NA>  <NA>
5  13    KDR   <NA>  <NA>   <NA>    KDR   <NA>   KDR    KDR   <NA>  <NA>   <NA>  <NA>  <NA>
6  14  CSF1R   <NA>  <NA>   <NA>  CSF1R   <NA>  <NA>   <NA>   <NA>  <NA>  CSF1R  <NA>  <NA>
data <- data[,-1] #去除第一列
data <- t(data) #转换行列
write.csv(data,file = 'data.csv') #保存文件

使用记事本工具打开data.csv文件

 可以看到文件中含有大量的NA，这会影响后续的分析，所以我们要去除NA，这里我们使用替换工具，将,NA全部去除

 然后我们再读取这个文件

data <- read.csv('data.csv', row.names = 1) #读取文件，第一列作为列名
data <- data[,1:28] #截取1：28列
data <- as.data.frame(data) #转换格式
data <- rownames_to_column(data,var = 'patient') #把行名作为第一列，命名为patient
data_information <- merge(data,patient_information[,1:2],by = 'patient', all = F) #把patient作为合并列，并删除所有比对不上的行
PR <- data_information[which(data_information$siteBestResponse == 'PR'),] #提取所有PR行数据
SD_PD <- data_information[which(data_information$siteBestResponse == 'SD')&which(data_information$siteBestResponse == 'PD'),] #提取所有SD和PD行数据
row.names(PR) <- PR[,1] #把第一列作为行名
PR <- PR[,c(-1,-30)] #删除第一列和最后一列
> head(PR)
          V1     V2     V3     V4   V5    V6     V7    V8   V9   V10   V11    V12    V13
01002 PIK3CA    RET    APC   ESR1 BRAF ERBB2   NRAS                                     
01004 ARID1A  BRCA2   TP53  CCND3  RB1  MEN1  BRCA1   ALK SDHA PTCH1   ATM  CCND2 RICTOR
01006 ARID1A    KDR  CSF1R NOTCH1  RET  BRAF   TP53   NF1 ROS1 CCND3   MET    RB1   MEN1
01007 ARID1A PIK3CA NOTCH1    RET TP53 FGFR2  BRCA1 PTCH1  ATM  TSC2   SMO  ERBB4   FLT4
01013  BRCA2    KIT  VEGFA    APC FLT3  POLE RICTOR  CDK6 MSH6  FLT4 ERBB2  NTRK3   NRAS
01016 PIK3CA    KIT  CSF1R NOTCH1  RET   NF1  FGFR3 GNA11  RB1   ALK CCND2 RICTOR   TSC2
write.csv(PR,file = 'PR_gene.csv')
write.csv(SD_PD,file = 'SD_PD_gene.csv')

 4.ID转换

为了以防出bug，我们一般会把基因名转换为包对应的ID名，由于表格现在有些不规则，所以我们需要写个循环来转换ID。

PR_gene <- read.csv('data_PR.csv',header = T,row.names = 1)
PR_gene <- as.data.frame(t(PR_gene))
SD_PD_gene <- read.csv('data_PD_SD.csv',header = T,row.names = 1)
SD_PD_gene <- as.data.frame(t(SD_PD_gene))
for (i in 1:ncol(PR_gene)) {PR_gene[,i][1:length(bitr(PR_gene[,i],fromType = 'SYMBOL',toType = c('ENTREZID'),OrgDb='org.Hs.eg.db')[,2])] <- bitr(PR_gene[,i],fromType = 'SYMBOL',toType = c('ENTREZID'),OrgDb='org.Hs.eg.db')[,2]} #bitr函数就是基因名转换函数
PR_list <- as.list(PR_gene) ##compareCluster函数适合list形式的数据，所以这里我们转换一下表格

5.通路富集

这里我们开始富集通路

PR_pathway_KEGG <- compareCluster(PR_list, fun="enrichKEGG",
                                    organism="hsa", pvalueCutoff=1) #这是KEGG，设置pvalueCutoff=1可以导出所有富集到的通路
write.csv(PR_pathway_KEGG,file = 'PR_pathway_KEGG.csv')

PR_pathway_GO <- compareCluster(PR_list,
                        fun="enrichGO", 
                        OrgDb="org.Hs.eg.db", 
                        ont= "BP",
                        pvalueCutoff=1,
                        pAdjustMethod = "BH",
                        qvalueCutoff = 1) #这是GO富集
write.csv(PR_pathway_GO,file = 'PR_pathway_GO.csv')


#自定义通路富集
library(msigdbr) #加载GSEA官网的通路数据包
DmGO <- msigdbr(species="Homo sapiens", #物种名
                category="C2") #选择目录，可以通过官网查询自己想富集的通路的目录
PID_pathway <- DmGO[c(which(DmGO$gs_subcat == 'CP:PID'), 
                      which(DmGO$gs_subcat == 'CP')),] #通过自己需要富集的ID号提取通路
> head(PID_pathway)
# A tibble: 6 x 15
  gs_cat gs_subcat gs_name                        gene_symbol entrez_gene ensembl_gene human_gene_symb~
  <chr>  <chr>     <chr>                          <chr>             <int> <chr>        <chr>           
1 C2     CP:PID    PID_A6B1_A6B4_INTEGRIN_PATHWAY AKT1                207 ENSG0000014~ AKT1            
2 C2     CP:PID    PID_A6B1_A6B4_INTEGRIN_PATHWAY CASP7               840 ENSG0000016~ CASP7           
3 C2     CP:PID    PID_A6B1_A6B4_INTEGRIN_PATHWAY CD9                 928 ENSG0000001~ CD9             
4 C2     CP:PID    PID_A6B1_A6B4_INTEGRIN_PATHWAY CDH1                999 ENSG0000003~ CDH1            
5 C2     CP:PID    PID_A6B1_A6B4_INTEGRIN_PATHWAY COL17A1            1308 ENSG0000006~ COL17A1         
6 C2     CP:PID    PID_A6B1_A6B4_INTEGRIN_PATHWAY EGF                1950 ENSG0000013~ EGF             
# ... with 8 more variables: human_entrez_gene <int>, human_ensembl_gene <chr>, gs_id <chr>,
#   gs_pmid <chr>, gs_geoid <chr>, gs_exact_source <chr>, gs_url <chr>, gs_description <chr>
PR_pathway_PID <- compareCluster(PR_list, #此处要注意，由于我们自定义的pathway的基因没有转换为ID，所以此处的PR_list基因不能转换成ID
                                fun="enricher", #使用enricher函数
                                TERM2GENE=PID_pathway[,c(3,7)], #从上表可以看到有很多多余列，我们只取3和7列即可
                                pvalueCutoff=1,
                                pAdjustMethod = "BH",
                                qvalueCutoff = 1)
write.csv(PR_pathway_PID,file = 'PR_pathway_PID.csv')

6.可视化

用气泡图可视化结果

dotplot(PR_pathway_PID,showCategory=10,includeAll=TRUE)
dotplot(PR_pathway_GO,showCategory=10,includeAll=TRUE)
dotplot(PR_pathway_KEGG,showCategory=10,includeAll=TRUE)

示例图片

总结

总的来讲这个函数非常好用，适合样本多且希望分开富集基因找差异通路的时候。

这篇关于批量富集分析气泡图的画法的文章就介绍到这儿，希望我们推荐的文章对大家有所帮助，也希望大家多多支持为之网！