当前位置：网站首页>单细胞测序流程(单细胞rna测序)

单细胞测序流程(单细胞rna测序)

2022-07-31 15:38:00 【全栈程序员站长】

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

系列文章目录

文章目录

单细胞测序流程（一）简介与数据下载单细胞测序流程（二）数据整理单细胞测序流程（三）质控和数据过滤——Seurat包分析，小提琴图和基因离差散点图单细胞测序流程（四）主成分分析——PCA 单细胞测序流程（五）t-sne聚类分析和寻找marker基因单细胞测序流程（六）单细胞的细胞类型的注释单细胞测序流程（七）单细胞的细胞类型轨迹分析单细胞测序流程（八）单细胞的marker基因转化和GO富集分析
单细胞测序流程（九）单细胞的GO圈图

本期主讲内容——单细胞的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()