跟着Cell学单细胞转录组分析(十三):单细胞GSVA分析|这个包涵盖大多数物种

之前我们发过GSVA分析（有了这个包，猪的GSEA和GSVA分析也不在话下（第一集），【后续来了】有了这个包，猪的GSEA和GSVA分析也不在话下（第二集）），接着单细胞系列，重新说一下GSVA分析。

首先是数据集的问题，通常只做人和小鼠的，想要做其他物种的，苦于没有数据集，不过这里说的这个包，即使猪的GSVA分析也都可以做。

这里我们以小鼠为示例，也是常见物种的GSVA分析，结合单细胞的数据！GSVA其实就是对功能富集的量化，然后进行差异分析，寻找感兴趣的通路在样本中的变化。不同于常规的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)immuneT <- subset(immune, celltype=="T cells")#提取我们需要分析的细胞类型
immuneT <- as.matrix(immuneT@assays$RNA@counts)#提取count矩阵
meta <- immuneT@meta.data[,c("orig.ident", "sex", "age", "stim", "samples")]#分组信息，为了后续作图

之前一直苦于MSigDB数据库只有人的数据集，没有小鼠和其他物种的，网上也有教程如何根据基因同源性进行转化的，但是很麻烦，也不一定成功。还好有一个新的数据包被发现了，简直是福音---msigdbr包，可以做GSEA和GSVA。


#install.packages("msigdbr")
library(msigdbr)
msigdbr_species() #可以看到，这个包涵盖了20个物种
# A tibble: 20 x 2species_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 5gs_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，因为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在C5subcategory = "GO:BP") %>% dplyr::select(gs_name,gene_symbol)#这里可以选择gene symbol，也可以选择ID，根据自己数据需求来，主要为了方便
mouse_GO_bp_Set = mouse_GO_bp %>% split(x = .$gene_symbol, f = .$gs_name)#后续gsva要求是list，所以将他转化为list

表达矩阵数据有了，通路信息有了，就可以进行GEVA分析了，代码简单就一句！保存结果！

T_gsva <- gsva(expr = immuneT, gset.idx.list = mouse_GO_bp_Set,kcdf="Poisson", #查看帮助函数选择合适的kcdf方法 parallel.sz = 5)write.table(T_gsva, 'T_gsva.xls', row.names=T, col.names=NA, sep="\t")

接着差异分析可以用limma包，类似于转录组芯片数据分析流程。

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(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"))#设置颜色

整个流程就结束了，希望对你们的研究能有启发，GO、KEGG做多了，可以换着做一下GSVA分析！

跟着Cell学单细胞转录组分析(十三):单细胞GSVA分析|这个包涵盖大多数物种相关推荐

跟着Cell学单细胞转录组分析(五):单细胞转录组marker基因鉴定及细胞群注释
书接上回(跟着Cell学单细胞转录组分析(四):单细胞转录组测序UMAP降维聚类).完成数据降维和细胞聚类后,最主要的环节和工作就是确定各个细胞群,明确是什么类型的细胞,正群的细胞定群很关键,涉及到整 ...
跟着Cell学单细胞转录组分析(八):单细胞转录组差异基因分析及多组结果可视化
接着单细胞下游分析: 从Cell学单细胞转录组分析(一):开端!!! 跟着Cell学单细胞转录组分析(二):单细胞转录组测序文件的读入及Seurat对象构建跟着Cell学单细胞转录组分析(三):单细 ...
跟着Cell学作图|9.PPI分析(GeNets数据库)
9.PPI分析(GeNets数据库) "实践是检验真理的唯一标准." "复现是学习R语言的最好办法." DOI: 10.1016/j.cell.2020.05. ...
跟着Cell学作图 | 5.UMAP降维分析
跟着 Cell 学作图 | 5.UMAP降维分析 "实践是检验真理的唯一标准." "复现是学习R语言的最好办法." 2021.4.12_1 DOI: 10.10 ...
跟着Cell学作图 | 12.韦恩图(Vennerable包)
"实践是检验真理的唯一标准." "复现是学习生信的最好办法." 2021.4.12_1 DOI: 10.1016/j.cell.2020.05.032 这篇20 ...
跟着CELL学作图|1.火山图
跟着CELL学作图之火山图 "实践是检验真理的唯一标准." "复现是学习R语言的最好办法." DOI: 10.1016/j.cell.2020.05.032 这 ...
跟着 Cell 学作图 | 桑葚图(ggalluvial)
桑葚图今天我们复现一幅2021年Cell上Graphical abstract的图. Title:Human oral mucosa cell atlas reveals a stromal-neu ...
跟着Cell学作图 | 2.柱状图+误差棒+散点+差异显著性检验
跟着 Cell 学作图 | 2.柱状图+误差棒+散点 "实践是检验真理的唯一标准." "复现是学习R语言的最好办法." 2021.4.12_1 DOI: 10. ...
跟着 Cell 学作图 | 4.小提琴图
跟着 Cell 学作图 | 4.小提琴图 "实践是检验真理的唯一标准." "复现是学习R语言的最好办法." DOI: 10.1016/j.cell.2020.0 ...
跟着 Cell 学作图 | 3.箱线图+散点+差异显著性检验
跟着 Cell 学作图 | 3.箱线图+散点+差异显著性检验 "实践是检验真理的唯一标准." "复现是学习R语言的最好办法." DOI: 10.1016/j.c ...

跟着Cell学单细胞转录组分析(十三):单细胞GSVA分析|这个包涵盖大多数物种

跟着Cell学单细胞转录组分析(十三):单细胞GSVA分析|这个包涵盖大多数物种相关推荐

最新文章

热门文章