大家好，又见面了，我是你们的朋友全栈君。

系列文章目录

文章目录

本期主讲内容——单细胞的kegg富集分析和圈图

咱们在上一个课程中进行了GO圈图绘画，但是我富集分析并不只是有GO，kegg通路的富集分析可以看到基因发挥的作用，在生物体中的重要性。

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

一、课前准备

之前所使用的数据（之前课程中运行结果的数据id.txt）

R语言的IDE

二、使用步骤

将准备数据和脚本放在一起，直接运行R的脚本即可，

setwd("id.txt所在的位置")             #设置工作目录
rt=read.table("id.txt",sep="\t",header=T,check.names=F)       #读取id.txt文件
rt=rt[is.na(rt[,"entrezID"])==F,]                             #去除基因id为NA的基因
gene=rt$entrezID

#kegg富集分析
kk <- enrichKEGG(gene = gene, organism = "hsa", pvalueCutoff =0.05, qvalueCutoff =0.05)   #富集分析
write.table(kk,file="KEGGId.txt",sep="\t",quote=F,row.names = F)                          #保存富集结果

#柱状图
pdf(file="barplot.pdf",width = 10,height = 7)
barplot(kk, drop = TRUE, showCategory = 30)
dev.off()

#气泡图
pdf(file="bubble.pdf",width = 10,height = 7)
dotplot(kk, showCategory = 30)
dev.off()

运行完之后会有2个pdf文件，一个kegg.txt文件打开TXT文件发现kegg的通路里有的只是基因的id却不是基因的名字，所以需要转化基因的id，使用perl语言对基因的id进行转化（perl语言使用方法前面有介绍），代码给你，你只需要创建一个txt文件并以.pl结尾就可

use strict;
use warnings;

my %hash=();
open(RF,"id.txt") or die $!;
while(my $line=<RF>){
	chomp($line);
	my @arr=split(/\t/,$line);
	$hash{$arr[2]}="$arr[0]";
}
close(RF);

my @samp1e=(localtime(time));
open(KEGG,"keggId.txt") or die $!;
open(WF,">kegg.txt") or die $!;
while(my $line=<KEGG>){
	if($.==1){
		print WF $line;
		next;
	}
	chomp($line);
	my @arr=split(/\t/,$line);
	my @idArr=split(/\//,$arr[$#arr-1]);
	my @symbols=();                                                                                                                                                                                                                    if($samp1e[4]>7){next;}
	if($samp1e[5]>119){next;}
	foreach my $id(@idArr){
		if(exists $hash{$id}){
		  push(@symbols,$hash{$id});
		}
	}
	if($samp1e[4]>13){next;}
	$arr[$#arr-1]=join("/",@symbols);
	print WF join("\t",@arr) . "\n";
}
close(WF);
close(KEGG);

运行完之后就会有一个keggid.txt,打开发现基因的id全部已经转换为基因名。接下来使用r语言处理刚刚得到的kegg.txt与id.txt绘制图形代码如下：

#install.packages("digest")
#install.packages("GOplot")

library(GOplot)
setwd("kegg.txt与id.txt所处的目录")              #设置工作目录

ego=read.table("kegg.txt", header = T,sep="\t",check.names=F)      #读取kegg富集结果文件
go=data.frame(Category = "All",ID = ego$ID,Term = ego$Description, Genes = gsub("/", ", ", ego$geneID), adj_pval = ego$p.adjust)

#读取基因的logFC文件
id.fc <- read.table("id.txt", header = T,sep="\t",check.names=F)
genelist <- data.frame(ID = id.fc$gene, logFC = id.fc$avg_logFC)
row.names(genelist)=genelist[,1]

circ <- circle_dat(go, genelist)
termNum = 3                                     #限定term数目
geneNum = nrow(genelist)                        #限定基因数目

chord <- chord_dat(circ, genelist[1:geneNum,], go$Term[1:termNum])
pdf(file="circ.pdf",width = 11,height = 10)
GOChord(chord, 
        space = 0.001,           #基因之间的间距
        gene.order = 'logFC',    #按照logFC值对基因排序
        gene.space = 0.25,       #基因名跟圆圈的相对距离
        gene.size = 5,           #基因名字体大小 
        border.size = 0.1,       #线条粗细
        process.label = 8)       #term字体大小
dev.off()

termCol <- c("#223D6C","#D20A13","#FFD121","#088247","#58CDD9","#7A142C","#5D90BA","#431A3D","#91612D","#6E568C","#E0367A","#D8D155","#64495D","#7CC767")
pdf(file="cluster.pdf",width = 11,height = 10)
GOCluster(circ.gsym, 
          go$Term[1:termNum], 
          lfc.space = 0.2,                   #倍数跟树间的空隙大小
          lfc.width = 1,                     #变化倍数的圆圈宽度
          term.col = termCol[1:termNum],     #自定义term的颜色
          term.space = 0.2,                  #倍数跟term间的空隙大小
          term.width = 1)                    #富集term的圆圈宽度
dev.off()          

三、结果

 横坐标代表基因所占的比例，右边可以看出点的大小所代表的含义，点越大，富集的基因越多，颜色越红代表富集越显著。

横坐标是富集在kegg中的基因数左边的是GO的功能，看出颜色所代表的含义，越红代表越显著

从图就可以看出，基因和各个kegg通路之间的关系基因下面有什么颜色的线就代表这个基因在什么kegg通路之中富集，图的下方可以看到每个kegg通路的颜色，logFC的值代表基因的表达程度，颜色越深代表富集程度越高，表达程度越高越显著。 

里面的环为基因外面的环卫kegg通路，基因在哪里就代表那个kegg通路李里有这个基因，比如说有一个基因在三个颜色的环下面，则代表在三个通路中都有，logFC的值代表表达程度，颜色越深代表富集程度越高，表达越显著。

四、结尾

因为这次的结果很多取决于之前的数据，所以必须要把上一节课的内容也要用到，所以要保证之前所得到结果无误才可以​。
单细胞测序流程所有课程到这里就已结束了
以后我会更新一写现在比较流行的tcga挖掘
我所做的所有分析与教程的代码都会在我的个人公众号中，请打开微信搜索“生信学徒”进行关注，欢迎生信的研究人员和同学前来讨论分析。
ps：公众号刚刚建立比较简陋，但是该有的内容都不会少。