批量富集分析气泡图的画法
作者：互联网
目录
前言
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 NA3.处理数据数据
我们的表格是不符合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)示例图片
总结
总的来讲这个函数非常好用，适合样本多且希望分开富集基因找差异通路的时候。
标签：富集,PR,画法,NA,ARID1A,PID,S01001,gene,气泡 来源： https://blog.csdn.net/weixin_44070677/article/details/120248340