有了这个包，猪的GSEA和GSVA分析也不在话下（第一集） (qq.com)
GSVA全称Gene Set Variation Analysis
作用
对功能富集的量化，然后进行差异分析，寻找感兴趣的通路在样本中的变化。不同于常规的GO、KEGG受差异基因阈值的影响，GSEA受实验分组的影响，GSVA能够对通路量化，看感兴趣通路在多组之间的变化！
分析过程
- 加载和安装必要的包
library(Seurat)
#source("http://bioconductor.org/biocLite.R")
#biocLite("GSVA"）
library(GSVA)
library(tidyverse)
library(ggplot2)
library(clusterProfiler)
library(org.Mm.eg.db)
library(dplyr)
- 导入单细胞转录组数据
load("D:/SC_re.RData") #加载数据集
T.sc <- subset(SC_re, celltype=="T cells")#用subset函数提取我们需要分析的细胞类型
T.exp <- as.matrix(Rods.sc@assays$RNA@counts)#提取count矩阵，一定要是as.matrix，如果开始时是dataframe，后期GSVA分析时也要转为matrix
meta <- T.sc.sc@meta.data[,c("orig.ident", "sex", "age", "stim", "samples")]#分组信息，为了后续作图
使用原始的count 或者 data数据，不能使用高变基因做
MSigDB数据库只有人的数据集,msigdbr包被发现了，可以做GSEA和GSVA
#install.packages("msigdbr")
library(msigdbr)
msigdbr_species() #可以看到，这个包涵盖了20个物种
# A tibble: 20 x 2
species_name species_common_name
<chr> <chr>
1 Anolis carolinensis Carolina anole, green anole
2 Bos taurus bovine, cattle, cow, dairy cow, domestic cattle, dome~
3 Caenorhabditis elegans roundworm
4 Canis lupus familiaris dog, dogs
5 Danio rerio leopard danio, zebra danio, zebra fish, zebrafish
6 Drosophila melanogaster fruit fly
7 Equus caballus domestic horse, equine, horse
8 Felis catus cat, cats, domestic cat
9 Gallus gallus bantam, chicken, chickens, Gallus domesticus
10 Homo sapiens human
11 Macaca mulatta rhesus macaque, rhesus macaques, Rhesus monkey, rhesu~
12 Monodelphis domestica gray short-tailed opossum
13 Mus musculus house mouse, mouse
14 Ornithorhynchus anatinus duck-billed platypus, duckbill platypus, platypus
15 Pan troglodytes chimpanzee
16 Rattus norvegicus brown rat, Norway rat, rat, rats
17 Saccharomyces cerevisiae baker's yeast, brewer's yeast, S. cerevisiae
18 Schizosaccharomyces pombe 9~ NA
19 Sus scrofa pig, pigs, swine, wild boar
20 Xenopus tropicalis tropical clawed frog, western clawed frog
查看下小鼠的基因集，是否与MSigDB数据库一样
mouse <- msigdbr(species = "Mus musculus")
mouse[1:5,1:5]
# A tibble: 5 x 5
gs_cat gs_subcat gs_name gene_symbol entrez_gene
<chr> <chr> <chr> <chr> <int>
1 C3 MIR:MIR_Legacy AAACCAC_MIR140 Abcc4 239273
2 C3 MIR:MIR_Legacy AAACCAC_MIR140 Abraxas2 109359
3 C3 MIR:MIR_Legacy AAACCAC_MIR140 Actn4 60595
4 C3 MIR:MIR_Legacy AAACCAC_MIR140 Acvr1 11477
5 C3 MIR:MIR_Legacy AAACCAC_MIR140 Adam9 11502
table(mouse$gs_cat) #查看目录，与MSigDB一样，包含9个数据集
###C1 C2 C3 C4 C5 C6 C7 C8 H
20049 533767 795972 92353 1248327 30556 988358 109328 7411
- 分析GO
本例中，我们要分析GO，因为mouse文件包含了所有的基因集，所以要查看GO在哪里，然后将需要的文件提出来。???
table(mouse$gs_subcat)
CGN CGP CM CP
167344 42770 376981 49583 3881
CP:BIOCARTA CP:KEGG CP:PID CP:REACTOME CP:WIKIPATHWAYS
4860 13694 8196 98232 27923
GO:BP GO:CC GO:MF HPO IMMUNESIGDB
660368 100991 105717 381251 944068
MIR:MIR_Legacy MIR:MIRDB TFT:GTRD TFT:TFT_Legacy VAX
34118 372658 235886 153310 44290
mouse_GO_bp = msigdbr(species = "Mus musculus",
category = "C5", #GO在C5
subcategory = "GO:BP") %>%
dplyr::select(gs_name,gene_symbol)#这里可以选择gene symbol，也可以选择ID，根据自己数据需求来，主要为了方便
head(mouse_GO_bp)
# A tibble: 6 x 2
gs_name gene_symbol
<chr> <chr>
1 GOBP_10_FORMYLTETRAHYDROFOLATE_METABOLIC_PROCESS Aldh1l1
2 GOBP_10_FORMYLTETRAHYDROFOLATE_METABOLIC_PROCESS Aldh1l2
3 GOBP_10_FORMYLTETRAHYDROFOLATE_METABOLIC_PROCESS Mthfd1
4 GOBP_10_FORMYLTETRAHYDROFOLATE_METABOLIC_PROCESS Mthfd1l
5 GOBP_10_FORMYLTETRAHYDROFOLATE_METABOLIC_PROCESS Mthfd2l
6 GOBP_2_OXOGLUTARATE_METABOLIC_PROCESS Aadat
mouse_GO_bp_Set = mouse_GO_bp %>% split(x = .$gene_symbol, f = .$gs_name)#后续gsva要求是list，所以将他转化为list
- GEVA分析
有表达矩阵，有通路信息
T_gsva <- gsva(expr = T.exp,
gset.idx.list = mouse_GO_bp_Set,
kcdf="Poisson", #查看帮助函数选择合适的kcdf方法
parallel.sz = 5)
#Setting parallel calculations through a MulticoreParam back-end
#with workers=5 and tasks=100.
#Estimating GSVA scores for 7410 gene sets.
#Estimating ECDFs with Poisson kernels
#Estimating ECDFs in parallel
#iteration: 100
#|=============================================================================| 100%
#查看分析结果：变成了GOBP的表达值了
head(MG_gsva[1:4, 1:4])
#GOBP_10_FORMYLTETRAHYDROFOLATE_METABOLIC_PROCESS 0.3971866
#GOBP_2_OXOGLUTARATE_METABOLIC_PROCESS -0.3306133
#GOBP_2FE_2S_CLUSTER_ASSEMBLY -0.2476997
#GOBP_3_PHOSPHOADENOSINE_5_PHOSPHOSULFATE_BIOSYNTHETIC_PROCESS 0.1370572
- 保存结果
write.table(T_gsva, 'T_gsva.xls', row.names=T, col.names=NA, sep="\t") - 热图结果展示
library(pheatmap)
library(patchwork)
A <- MG_gsva[1:50,] #为了方便展示，我们只展示前50行
pheatmap(A, show_rownames=1, show_colnames=0)
【后续来了】有了这个包，猪的GSEA和GSVA分析也不在话下（第二集） (qq.com)
差异分析
首先将之前保存的文件读入
T_gsva <- read.csv("T_gsva.csv", header = T, row.names = 1)
分析得到的数据结构是行为GO terms， 列为样品（单细胞中为单个细胞）
group <- c(rep("control", 50), rep("test", 71)) %>% as.factor()#设置分组，对照在前
desigN <- model.matrix(~ 0 + group) #构建比较矩阵
colnames(desigN) <- levels(group)
fit = lmFit(test_control, desigN)
fit2 <- eBayes(fit)
diff=topTable(fit2,adjust='fdr',coef=2,number=Inf)#校准按照需求修改 ？topTable
write.csv(diff, file = "Diff.csv")
柱状图可视化
up <- c("GOBP_EGG_ACTIVATION",
"GOBP_TENDON_DEVELOPMENT",
"GOBP_SOMITE_SPECIFICATION",
"GOBP_THREONINE_CATABOLIC_PROCESS",
"GOBP_REGULATION_OF_GLUTAMATE_RECEPTOR_CLUSTERING",
"GOBP_NEGATIVE_CHEMOTAXIS",
"GOBP_NEGATIVE_REGULATION_OF_FAT_CELL_PROLIFERATION",
"GOBP_REGULATION_OF_T_HELPER_17_CELL_LINEAGE_COMMITMENT",
"GOBP_REGULATION_OF_ANTIMICROBIAL_HUMORAL_RESPONSE")
down <- c("GOBP_DETERMINATION_OF_PANCREATIC_LEFT_RIGHT_ASYMMETRY",
"GOBP_MITOTIC_DNA_REPLICATION",
"GOBP_EOSINOPHIL_CHEMOTAXIS",
"GOBP_NEUTROPHIL_MEDIATED_CYTOTOXICITY",
"GOBP_POTASSIUM_ION_EXPORT_ACROSS_PLASMA_MEMBRANE",
"GOBP_REGULATION_OF_LEUKOCYTE_MEDIATED_CYTOTOXICITY",
"GOBP_REGULATION_OF_SEQUESTERING_OF_ZINC_ION",
"GOBP_ENDOTHELIN_RECEPTOR_SIGNALING_PATHWAY",
"GOBP_PRE_REPLICATIVE_COMPLEX_ASSEMBLY_INVOLVED_IN_CELL_CYCLE_DNA_REPLICATION",
"GOBP_ESTABLISHMENT_OF_PLANAR_POLARITY_OF_EMBRYONIC_EPITHELIUM")
TEST <- c(up,down)
diff$ID <- rownames(diff)
Q <- diff[TEST,]
group1 <- c(rep("treat", 9), rep("control", 10))
df <- data.frame(ID = Q$ID, score = Q$t,group=group1 )
# 按照t score排序
sortdf <- df[order(df$score),]
sortdf$ID <- factor(sortdf$ID, levels = sortdf$ID)#增加通路ID那一列
用ggplot画图（ggplot-YYDS）
ggplot(sortdf, aes(ID, score,fill=group)) + geom_bar(stat = 'identity',alpha = 0.7) +
coord_flip() +
theme_bw() + #去除背景色
theme(panel.grid =element_blank())+
theme(panel.border = element_rect(size = 0.6))+
labs(x = "",
y="t value of GSVA score")+
scale_fill_manual(values = c("#008020","#08519C"))#设置颜色
跟着Cell学单细胞转录组分析(十三):单细胞GSVA分析|这个包涵盖大多数物种 - 简书 (jianshu.com)
