作为一名基础医学生，大家有没有这样的苦恼，那就是拿不到临床病人的样本，只能做小鼠！！！做就做呗，关键是好多基因集都是针对人来做的，于是乎
做小鼠的GSEA,GSVA难上加难。想到每次洗头掉下的大把头发，我表示
好在总是会有一些牛逼嘎拉丝的大神来解决我们的问题，只是我们不知道如何找到它们，做生信不怕问题，就怕搜不到答案！！
下面我就为大家展示如何基于小鼠基因集做GSVA
1.读取单细胞数据
rm(list=ls())
setwd("D:\\mouse_GSEA")
#saveRDS(sub_Treg, file = "sub_Treg.rds")
sub_Treg <- readRDS("sub_Treg.rds")
table(sub_Treg$sample_type)
####不用运行，查看数据分布#####
a<- sub_Treg@assays$RNA@scale.data
write.csv(a,"data\\scale.data.csv")
b<- sub_Treg@assays$RNA@data
write.csv(b,"data\\normalized.csv")
b<- sub_Treg@assays$RNA@counts
write.csv(b,"data\\counts.csv")
2.准备数据
library(GSVA)
library(Seurat)
#sub_Treg@assays$RNA@counts未标准化的数据
expr <- as.data.frame(sub_Treg@assays$RNA@counts)
library(ggplot2)
library(clusterProfiler)
library(org.Mm.eg.db)
library(dplyr)
3.构建参考数据集
expr=as.matrix(expr)
class(expr)
#构建参考数据集
library(msigdbr)
#小鼠所有的基因集
m_df = msigdbr(species = "Mus musculus")
##2.3 查看基因集类别：
a <- m_df %>% dplyr::distinct(gs_cat, gs_subcat) %>% dplyr::arrange(gs_cat, gs_subcat)
View(a)
##2.5 检索鼠类C2 (curated) CGP (chemical and genetic perturbations)基因集：
m_df = msigdbr(species = "Mus musculus", category = "C2", subcategory = "CGP")
##KEGG
m_df = msigdbr(species = "Mus musculus", category = "C2", subcategory = "CP:KEGG")
#C5 GO:BP CC MF HPO
m_df = msigdbr(species = "Mus musculus", category = "C5")
##2.4 检索鼠类的hallmark 基因集：
m_df = msigdbr(species = "Mus musculus", category = "H")
## 免疫相关通路，C7
m_df = msigdbr(species = "Mus musculus", category = "C7")
#构建参考数据集
msigdbr_list = split(x = m_df$gene_symbol, f = m_df$gs_name)
4.GSVA
运行GSVA#如果是标准化后的数据用默认参数，此处用raw_count,耗时较长
kegg2 <- gsva(expr, msigdbr_list, kcdf="Poisson",method = "gsva",parallel.sz=1)
head(kegg2)
5.差异分析（PS：本来做到这就够了，但好人做到底，那就再来个差异分析）
#差异分析
library(limma)
rt <- kegg2
head(rt)
logFCcutoff=0.32 #也可以0，一般0.2-0.3，可以大
adjPvalueCutoff=0.05
## 1.构建比较矩阵
type=c( rep("con",181),rep("treat",177) )
## 1.构建比较矩阵
design=model.matrix(~ type)
## 比较矩阵命名
colnames(design)=c("con", "treat")
##2.线性模型拟合
fit=lmFit(rt, design)
##3.贝叶斯检验
fit=eBayes(fit)
#全部差异分析结果
all=topTable(fit, coef="con", number=Inf,adjust.method="holm")
all=rbind(id=colnames(all),all)
write.table(all,file="treg.C7_all.txt",sep="\t",quote=F,col.names=F)
#有统计学意义的结果
diff <- topTable(fit, coef="con", number=Inf,
p.value=adjPvalueCutoff, adjust="holm", lfc=logFCcutoff)
diffName=row.names(diff)
diff=rbind(id=colnames(diff),diff)
write.csv(diff,file="treg.C7_diff.csv",quote=F,col.names=F)
下面就是用火山图，热图啥的可视化了，有很多帖子可以参考。。
好了，第一次展示就此结束，欢迎做单细胞测序自己分析的小伙伴简信我，批评指正
如果你自己做单细胞测序，可以关注公众号“生信宝库”，里面经常发一些单细胞测序的文章，我就是从那学到的。
末了末了，我再感叹一句，其实做科研真的值得全身心付出才会有收获，你想想，就拿小鼠来说，人家可是熬着夜陪咱做实验，搞不好再重复个两三次，得死多少同胞鼠弟，最后希望大家都能做出自己想要的结果。
网友评论