seurat使用findallmarkers 得到的差异基因列表进行富集分析clusterprolifer-2 调成fc值为0.69
#######差异分析
###########################33
library(clusterProfiler)
library(Seurat)
library(SeuratWrappers)
library(ggplot2)
library(dplyr)
library(stringr)
library(openxlsx)
library(org.Mm.eg.db)
library(“reshape”)
library(“RColorBrewer”)
getwd()
file=“G:/silicosis/sicosis/silicosis_ST/yll/0214/harmony_cluster/d_all/r=0.6_11clusters”
dir.create(file)
setwd(file)
getwd()
?loadWorkbook
#library(xlsx) #比openxlsx读取文件慢多了!!!!!
##再三检查文件路径!!!!!!!!!!!!
markers = read.xlsx(“G:/silicosis/sicosis/silicosis_ST/yll/0214/harmony_cluster/d_all/silicosis_ST_harmony_r0.6_cluster_marker.xlsx”)
head(markers)
#View(markers)
names(table(markersKaTeX parse error: Expected 'EOF', got '#' at position 11: cluster)) #̲markers"Myofibroblast/vascular smooth muscle cell"<-“myo-fib”
markerscluster[markerscluster[markerscluster[markerscluster==“Myofibroblast/vascular smooth muscle cell”]=“fib-myo” #把第17个cluster的非法名字改成正常名字,因为给导出的xlsx命名时使用了custer的名字
names(table(markers$cluster))
k <- keys(org.Mm.eg.db, keytype = “ENTREZID”)
gene.list <- select(org.Mm.eg.db, keys = k, columns = c(“SYMBOL”, “ENSEMBL”), keytype = “ENTREZID”)## 73306 3
head(gene.list)
#自建函数 输入go富集条目 和文件的名称 就返回相应的pdf文件和xlsx文件
clusterProfilerOut = function(enrichRes,outfile){
enrichResfunction.gene.num=unlist(strsplit(enrichRes[,"BgRatio"],"/"))[seq(1,length(unlist(strsplit(enrichRes[,"BgRatio"],"/"))),2)]enrichResfunction.gene.num = unlist(strsplit(enrichRes[,"BgRatio"],"/"))[seq(1,length(unlist(strsplit(enrichRes[,"BgRatio"],"/"))),2)] enrichResfunction.gene.num=unlist(strsplit(enrichRes[,"BgRatio"],"/"))[seq(1,length(unlist(strsplit(enrichRes[,"BgRatio"],"/"))),2)]enrichResGeneRatio = enrichResCount/as.numeric(enrichResCount / as.numeric(enrichResCount/as.numeric(enrichResfunction.gene.num) #计算generatio
write.xlsx(enrichRes,paste0(outfile, “.xlsx”),col.names=T,row.names=F)
if (dim(enrichRes)[1]>=10) {enrichRes.use = enrichRes[1:10,]} else{enrichRes.use = enrichRes[1:dim(enrichRes)[1],]}#显示的富集条目不超过10条
enrichRes.use = enrichRes.use[order(enrichRes.useKaTeX parse error: Expected 'EOF', got '#' at position 28: …ecreasing=T),] #̲按照generatio降序显示…Description = factor(enrichRes.useDescription,levels=rev(enrichRes.useDescription,levels=rev(enrichRes.useDescription,levels=rev(enrichRes.useDescription))
gap = (max(enrichRes.useGeneRatio)−min(enrichRes.useGeneRatio) - min(enrichRes.useGeneRatio)−min(enrichRes.useGeneRatio))/5
#制作并输出pdf
pdf(paste0(outfile, “.pdf”))
p = ggplot(enrichRes.use,mapping = aes(x=GeneRatio,y=Description))+
geom_point(aes(size=Count,color=p.adjust))+theme_bw()+
scale_color_gradient(low = “blue”, high = “red”)+
scale_x_continuous(expand = c(gap, gap))+ylab(NULL)+
scale_y_discrete(labels=function(x) str_wrap(x, width=60))
print§
dev.off()
}
getwd() #此时手动把all_cluster_markers.xlsx文件复制到工作目录下
path=getwd()
dir(path)
files=dir(path)
files #得到工作目录下的所有文件名
‘’’
path=“G:/silicosis/sicosis/silicosis_ST/yll/0214/harmony_cluster/d_all/r=0.6_11clusters_0.69fc”
dir.create(path)
setwd(path)
‘’’
getwd()
head(markers)
#library(xlsx)
for (i in files) {
#i=“all_cluster_markers.xlsx”
input = read.xlsx(paste(path, i, sep = “/”))#读取xlsx文件 注意有的时候需要加row.names 有的时候不需要加
head(input)
inputcluster[inputcluster[inputcluster[inputcluster==“Myofibroblast/vascular smooth muscle cell”]=“fib-myo” #把第17个cluster的非法名字改成正常名字,因为给导出的xlsx命名时使用了custer的名字
inputsymbol<−inputsymbol<-inputsymbol<−inputgene
head(input)
inputFeatureName.entrez=gene.list[match(inputFeatureName.entrez = gene.list[match(inputFeatureName.entrez=gene.list[match(inputsymbol, gene.listKaTeX parse error: Expected 'EOF', got '#' at position 20: …OL),"ENTREZID"]#̲在input文件的基础上添加s…change=ifelse(inputavglog2FC>0,"up","down")input.sub<−input[inputavg_log2FC>0,"up","down") input.sub<-input[inputavglog2FC>0,"up","down")input.sub<−input[inputp_val_adj<0.001,] #取出所有p值小于0.001的
print(head(input.sub))
dim(input.sub)
input.sub = na.omit(input[,c(“cluster”,“FeatureName.entrez”,“avg_log2FC”,“p_val_adj”)])
head(input.sub)
colnames(input.sub) = c(“cluster”,“gene”, “log2FC”, “PValue”)
print(head(input.sub))
subset_data=input.sub
head(subset_data)
cluster_number<-names(table(subset_data$cluster)) #获取文件中的cluster个数
length(cluster_number)
cluster_number
str(cluster_number)
names(table(subset_datacluster))[1]subsetdata[subsetdatacluster))[1] subset_data[subset_datacluster))[1]subsetdata[subsetdatacluster==“0”,]
##注意文件中cluster的命名 最好不要出现斜杠|等特殊字符 容易出错
#注意p值是否需要约束?
#注意是否需要正负号表示上下调?
#for (cluster_num in 0:(length(cluster_number)-1)) { #对于cluster为数字的
##只看上调的基因 ,注意修改名字!!!!!!!!!!!!!
subset_data = subset(input.sub, log2FC>0.69&PValue<0.05)
head(subset_data)
cluster_number<-names(table(subset_data$cluster)) #获取文件中的cluster个数
cluster_number
for (cluster_num in cluster_number) { ##对于cluster为字符的
ego=enrichGO(gene=subset_data[subset_data$cluster==cluster_num,]$gene, #!!!!!!全都要改!!!!!!!!!!!!!!!!!!!泉沟要改keyType = "ENTREZID", OrgDb = org.Mm.eg.db, ont = "ALL", readable = TRUE,pvalueCutoff = 1, qvalueCutoff = 1, pAdjustMethod = "BH")
res = ego[ego$p.adjust<0.05,]
print(dim(res))
print(table(res$ONTOLOGY))if(nrow(res) != 0){#go,使用自定义的函数画图 ,注意命名时候把i的 文件属性名去掉clusterProfilerOut(res, paste0(unlist(strsplit(i, ".xl"))[1],cluster_num,"_up_enrichGO"))
}ego.KEGG = enrichKEGG(gene=subset_data[subset_data$cluster==cluster_num,]$gene, organism = "mmu",keyType = "kegg", pvalueCutoff = 1, pAdjustMethod = "BH", qvalueCutoff = 1)
res = ego.KEGG[ego.KEGG$p.adjust<0.05,]
print(dim(res)) if(nrow(res) !=0){#keggclusterProfilerOut(res,paste0(unlist(strsplit(i, ".xl"))[1],cluster_num,"_up_enrichKEGG"))
}
}
##只看下调的基因 ,注意修改名字!!!!!!!!!!!!!
subset_data = subset(input.sub, log2FC<(-0.69)&PValue<0.05)
head(subset_data)
cluster_number<-names(table(subset_data$cluster)) #获取文件中的cluster个数
cluster_number
for (cluster_num in cluster_number) { ##对于cluster为字符的
ego=enrichGO(gene=subset_data[subset_data$cluster==cluster_num,]$gene, #!!!!!!全都要改!!!!!!!!!!!!!!!!!!!泉沟要改keyType = "ENTREZID", OrgDb = org.Mm.eg.db, ont = "ALL", readable = TRUE,pvalueCutoff = 1, qvalueCutoff = 1, pAdjustMethod = "BH")
res = ego[ego$p.adjust<0.05,]
print(dim(res))
print(table(res$ONTOLOGY))if(nrow(res) != 0){#go,使用自定义的函数画图 ,注意命名时候把i的 文件属性名去掉clusterProfilerOut(res, paste0(unlist(strsplit(i, ".xl"))[1],cluster_num,"_down_enrichGO"))
}ego.KEGG = enrichKEGG(gene=subset_data[subset_data$cluster==cluster_num,]$gene, organism = "mmu",keyType = "kegg", pvalueCutoff = 1, pAdjustMethod = "BH", qvalueCutoff = 1)
res = ego.KEGG[ego.KEGG$p.adjust<0.05,]
print(dim(res)) if(nrow(res) !=0){#keggclusterProfilerOut(res,paste0(unlist(strsplit(i, ".xl"))[1],cluster_num,"_down_enrichKEGG"))
}
}
##不区分上下调 #################################
subset_data = subset(input.sub, abs(log2FC)>0.69&PValue<0.05)
head(subset_data)
cluster_number<-names(table(subset_data$cluster)) #获取文件中的cluster个数
cluster_number
for (cluster_num in cluster_number) { ##对于cluster为字符的
ego=enrichGO(gene=subset_data[subset_data$cluster==cluster_num,]$gene, #!!!!!!全都要改!!!!!!!!!!!!!!!!!!!泉沟要改keyType = "ENTREZID", OrgDb = org.Mm.eg.db, ont = "ALL", readable = TRUE,pvalueCutoff = 1, qvalueCutoff = 1, pAdjustMethod = "BH")
res = ego[ego$p.adjust<0.05,]
print(dim(res))
print(table(res$ONTOLOGY))if(nrow(res) != 0){#go,使用自定义的函数画图 ,注意命名时候把i的 文件属性名去掉clusterProfilerOut(res, paste0(unlist(strsplit(i, ".xl"))[1],cluster_num,"_both_enrichGO"))
}ego.KEGG = enrichKEGG(gene=subset_data[subset_data$cluster==cluster_num,]$gene, organism = "mmu",keyType = "kegg", pvalueCutoff = 1, pAdjustMethod = "BH", qvalueCutoff = 1)
res = ego.KEGG[ego.KEGG$p.adjust<0.05,]
print(dim(res)) if(nrow(res) !=0){#keggclusterProfilerOut(res,paste0(unlist(strsplit(i, ".xl"))[1],cluster_num,"_both_enrichKEGG"))
}
}
}
dev.off()
path=“D:/ARDS_scripts_1012/ARDS/Step2_harmony_f200_R3/ws/Fibroblast/0.5/sepsis_Fibroblast_r0.5.rds”
load(path)
table(subset_datastim,subsetdatastim,subset_datastim,subsetdataseurat_clusters)
DimPlot(subset_data,group.by = “stim”)
DimPlot(subset_data,split.by = “stim”,
group.by = “seurat_clusters”,
label = TRUE)
seurat使用findallmarkers 得到的差异基因列表进行富集分析clusterprolifer-2 调成fc值为0.69相关推荐
- 从seurat的findallmarkers得到的差异基因 进行富集分析clusterprolifer
library(openxlsx)与library(xlsx)两个包经常出问题,报错往往都是他俩 建议只使用openxlsx 更快! #######差异分析 ##################### ...
- NBIS单细胞教程:差异基因(五)
Seurat_StepV_DEG 此文章是通过学习瑞典国家生物信息学基础设施(NBIS) 所开放的单细胞分析教程加上网上所查找的资料,自身的理解所形成的,可能会有不足之处.该部分是对不同集群之间的差异 ...
- python 分析两组数据的差异_R语言limma包差异基因分析(两组或两组以上)
使用limma包进行差异基因分析时,做最多的是两分类的,例如control组和disease组,但也会碰到按照序列进行的分组.这时,如果逐一使用两两比较求差异基因则略显复杂.其实开发limma包的大神 ...
- DESeq2差异基因分析和批次效应移除
差异基因鉴定 基因表达标准化 不同样品的测序量会有差异,最简单的标准化方式是计算 counts per million (CPM),即原始reads count除以总reads数乘以1,000,000 ...
- 生物信息学入门 根据表达矩阵和差异表达基因列表制作差异表达矩阵
根据表达矩阵做完差异分析之后,就要将差异表达基因的表达情况从表达矩阵中提取出来,制作差异表达矩阵.这个过程非常简单,但是也有一些人问到,就整理了这个教程. 1. 数据准备:原始表达矩阵和差异表达基因 ...
- 差异基因 p log2foldchange_拟南芥的基因ID批量转换?差异基因,GO/KEGG数据库注释(转录组直接送你全套流程)...
新手遇到的问题都是类似的,比如批量ID转换 虽然我写过大量的教程:ID转换大全 不过都需要R基础,因为是大批量转换啊! 但热心肠的植物生物信息学教学大佬还是友善的给出了解决方案 我也狗尾续貂制作了 ...
- NAR|北大/中科院计算所团队发布基因功能富集分析平台KOBAS-i
近日,国际知名期刊<核酸研究>(Nucleic Acids Research,IF:16.971)在线发表了北京大学孔雷课题组与中国科学院计算技术研究所赵屹研究员课题组合作开发的基因功能富 ...
- ClusterProfiler在线基因集富集分析,支持自定义基因集、任意物种
为什么pathway富集分析结果没有我感兴趣的通路? GO和KEGG富集分析使用差异基因(上调基因,下调基因,或者上下调合起来的基因)作为输入,使用超几何分布等算法计算显著富集的GO term或者通路 ...
- 单基因gsea_GSEA:基因集富集分析和ssGSEA:单样本基因集富集分析
传统富集分析(基于超几何分布或者Fisher精确检验):关注一列差异基因是否是随机分布在某一感兴趣的基因集中(某通路的基因) 得到通路富集的结果时: (1).一条通路中既有上调基因又有下调基因,无法确 ...
最新文章
- 右键 Dos在这里 删除
- 浅析Microsoft .net PetShop程序中的购物车和订单处理模块(Profile技术,异步MSMQ消息)转...
- Win2003系统安全设置
- adb提取安装的apk
- 近业务=困死在一条船上?
- MyBatis中多表查询(N+1方式)
- c#字符相似度对比通用类
- note_maven的基本使用
- 我该学习哪个人工智能系统
- tomcat的acceptCount、maxThreads、connectionTimeout参数调整
- 完美解决ERROR 1819 (HY000): Your password does not satisfy the current policy requirements
- java基础总结02-语言基础
- Python网抓 2021年 获取全部沪深港股ETF股票信息 东方财富
- Linux debian解压和压缩.rar文件教程
- Recurrent Feature Reasoning for Image Inpainting解读
- IBM X3650服务器使用说明
- Wireshark分析sql布尔盲注流量包
- adb命令读取Android手机内存卡文件
- ggplot2 去掉网格
- BACnet IP通讯方式组网步骤
热门文章
- java php 性能比较_JAVA和PHP的优劣对比
- MySQL之详解索引
- 实力凸显 | 思迈特软件入选“2022中国软件150强“等三大重磅榜单
- 使用python tkinter做window窗体界面程序,以及python多线程处理解决tk界面卡死
- 《巨人不死密码》读书笔记
- Android下最好用的开源HTTP 服务器
- ​TOIS 2022 | 面向用户和商品双侧流行性偏差的自适应公平推荐方法
- Python functools详解
- 【Python SMTP/POP3/IMAP】零基础也能轻松掌握的学习路线与参考资料
- webstorm提示Page'....'页面未经授权的警告解决方法