大家好，又见面了，我是你们的朋友全栈君。
系列文章目录
文章目录
- 单细胞测序流程（一）简介与数据下载 单细胞测序流程（二）数据整理 单细胞测序流程（三）质控和数据过滤——Seurat包分析，小提琴图和基因离差散点图 单细胞测序流程（四）主成分分析——PCA 单细胞测序流程（五）t-sne聚类分析和寻找marker基因 单细胞测序流程（六）单细胞的细胞类型的注释 单细胞测序流程（七）单细胞的细胞类型轨迹分析 单细胞测序流程（八）单细胞的marker基因转化和GO富集分析
- 单细胞测序流程（九）单细胞的GO圈图
本期主讲内容——单细胞的kegg富集分析和圈图
咱们在上一个课程中进行了GO圈图绘画，但是我富集分析并不只是有GO，kegg通路的富集分析可以看到基因发挥的作用，在生物体中的重要性。
提示：以下是本篇文章正文内容，下面案例可供参考
一、课前准备
之前所使用的数据（之前课程中运行结果的数据id.txt）
R语言的IDE
二、使用步骤
将准备数据和脚本放在一起，直接运行R的脚本即可，
代码语言：javascript复制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结尾就可
代码语言：javascript复制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绘制图形代码如下：
代码语言：javascript复制#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挖掘
发布者：全栈程序员栈长，转载请注明出处：https://javaforall.cn/127991.html原文链接：https://javaforall.cn
