有了这个包，猪的GSEA和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)

导入单细胞转录组数据：

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")]#分组信息，为了后续作图

这里有个重要的问题就是，单细胞GSVA到底用什么数据，count or data？之前《单细胞天地》公众号依据测试过这个问题了，count数据和data数据做GSVA分析没有区别，但是不能使用高变基因去做，对结果的影响较大，还是使用原始的count 或者 data数据吧。

==================重点来了==================

之前一直苦于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，根据自己数据需求来，主要为了方便
head(mouse_GO_bp)
# A tibble: 6 x 2gs_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)

结果如下：

之后的差异分析，可视化我们下节再说！

更多精彩内容请关注我的公众号《KS科研服务与分享》

有了这个包，猪的GSEA和GSVA分析也不在话下（第一集）相关推荐

跟着Cell学单细胞转录组分析(十三):单细胞GSVA分析|这个包涵盖大多数物种
之前我们发过GSVA分析(有了这个包,猪的GSEA和GSVA分析也不在话下(第一集),[后续来了]有了这个包,猪的GSEA和GSVA分析也不在话下(第二集)),接着单细胞系列,重新说一下GSVA分析. ...
gsva gsea ssgsea gaochao 使用GSVA方法计算某基因集在各个样本的表现
傻傻分不清!GSEA & GSVA有啥差别?史上最全教程来了! - 知乎 (zhihu.com) 文章发表于2013年,GSVA: gene set variation analysis fo ...
富集分析原理和clusterProfiler包进行GO、KEGG富集分析详细说明
概念: 基因富集分析是指对于给定一组基因根据基因组注释信息(GO.KEGG)对基因进行聚类分析,即给定的基因是不是GO中的一个功能(或KEGG中的一个通路). 基因的功能富集的目的是说明给定的基因集对 ...
OTA--卡刷全包、差分升级包制作、分析（代码摘自Google）---2
OTA官网介绍: https://source.android.com/devices/tech/ota/ OTA–卡刷全包.差分升级包制作.分析 1 ．OTA升级流程框架标准的OTA升级流程包括一 ...
一文掌握GSEA通路富集分析，超详细教程！
生信宝典之前总结了一篇关于GSEA富集分析的推文--GSEA富集分析:从概念理解到界面实操,介绍了GSEA的定义.GSEA原理.GSEA分析.Leading-edge分析等,是全网最流行的原理+操作兼 ...
python解析pcap包已text格式输出_python分析pcap包
前两天需要分析一个pcap包,写了一段python脚本,将每个包的基本信息(源/目的MAC.源/目的IP.源/目的端口)提取出来. 在实现过程中为了省事用了dpkt开发包,不过只用了几个简单的函数,具 ...
r roc函数_如何处理R(pROC包)中的多类ROC分析？
例如,当我在R(pROC包)中使用multiclass.roc函数时,我训练了随机森林的数据集,这是我的代码: # randomForest & pROC packages should be ...
iOS马甲包开发招式及规避4.3方法合集
看了下上周的留言,有些开发者老是抱怨马甲包又被拒了,该如何上包才安全,我对这块也算略知一二,也有朋友是做这块的,一些规避手段我还是知晓一二,我今天结合了网上的资料以及几位朋友的一些意见,把这些经验分享 ...
tcpdump抓包IP地址，导入wireshark分析？又名《~来抓包吧/ 向着前辈の步伐、Start / ~从零开始の抓包全过程流水账实录/// 成功吧~实验/ ~》
本文关键字:虚拟机.镜像文件iso.Ubuntu.linux.tcpdump.apt-get.yum命令.lrzsz命令.Xshell.wireshark (不要被英文吓晕了,我只是想列出来假装一下很 ...

有了这个包，猪的GSEA和GSVA分析也不在话下（第一集）

有了这个包，猪的GSEA和GSVA分析也不在话下（第一集）相关推荐

最新文章

热门文章