文章/答案/技术大牛

发布

社区首页 >专栏 >scRNA分析|单细胞GSVA + limma差异分析-celltype分组？样本分组？

scRNA分析|单细胞GSVA + limma差异分析-celltype分组？样本分组？

生信补给站

发布于 2023-08-25 02:01:18

2.1K06

代码可运行

文章被收录于专栏：生信补给站生信补给站

运行总次数：6

代码可运行

单细胞数据完成差异分析后，可以根据结果进行后续的GO ,KEGG,GSEA富集分析，推荐使用clusterProfiler-R包，可参考 R|clusterProfiler-富集分析 clusterProfiler|GSEA富集分析及可视化。

此外还可以进行GSVA分析，基因集变异分析即GSVA（Gene set variation analysis）, 是一种非参数、无监督的分析方法，可以分析不同的目标基因集在不同样本中的富集程度。

输入文件比较简单，只需要 A：表达量矩阵 和 B：目标基因集 即可分析。

一载入R包数据

1, 获取表达矩阵

如果想计算celltype的GSVA结果，可以使用 AverageExpression 函数计算不同celltype之间的表达量均值矩阵；

如果计算每个细胞的GSVA结果，直接提取表达量矩阵即可；

library(Seurat) 
#source("http://bioconductor.org/biocLite.R")
#biocLite("GSVA"）
library(GSVA)
library(tidyverse)
library(org.Hs.eg.db)

#加载单细胞数据
load("sce.anno.RData")

#计算各celltype的表达量均值
Idents(sce2) <- "Anno" 
expr <- AverageExpression(sce2, assays = "RNA", slot = "data")[[1]]
expr <- expr[rowSums(expr)>0,]  #过滤细胞表达量全为零的基因
expr <- as.matrix(expr)
head(expr)

                  Epi   Fibroblast         un           T      Endo     Myeloid
OR4F5      2.762724e-04 0.0000000000 0.00000000 0.000000000 0.0000000 0.000000000
AL627309.1 4.399548e-02 0.0080237642 0.01036203 0.007629199 0.0000000 0.047063397
AL627309.5 0.000000e+00 0.0008855168 0.00000000 0.000000000 0.0000000 0.000000000
AP006222.2 1.068307e-01 0.2828329948 0.12706163 0.050051575 0.1942506 0.244211331
LINC01409  1.299390e-04 0.0005598456 0.00000000 0.000000000 0.0000000 0.000000000
FAM87B     7.746596e-05 0.0026794766 0.02404758 0.000000000 0.0000000 0.001711278

单细胞表达量多为0，可以先过滤一下在所有细胞中均无表达的基因。

2，获取目标基因集

根据自己的需要选择MSigDB数据库中的基因集

2.1 手动下载

进入http://www.gsea-msigdb.org/gsea/msigdb/index.jsp后选择需要下载的基因集，然后使用R读取下载好的gmt格式的文件。

2.2 msigdbr包

直接使用msigdbr包内置好的基因集，含有多个物种 以及 多个基因集，通过参数选择物种以及数据集，较为方便。推荐！

library(msigdbr)
msigdbr_species() #列出有的物种

#选择基因集合
human_KEGG = msigdbr(species = "Homo sapiens", #物种
                      category = "C2",
                     subcategory = "KEGG") %>% 
  dplyr::select(gs_name,gene_symbol)#这里可以选择gene symbol或者ID
human_KEGG_Set = human_KEGG %>% split(x = .$gene_symbol, f = .$gs_name)#list

# A tibble: 20 × 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, domestic cow, ox, oxen
 3 Caenorhabditis elegans          NA                                                                     
 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, rhesus monkeys         
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 972h- NA                                                                     
19 Sus scrofa                      pig, pigs, swine, wild boar                                            
20 Xenopus tropicalis              tropical clawed frog, western clawed frog

A：如果你的研究是其中的物种就可以无缝做GSEA 和 GSVA了。

B：如果研究的物种不在其中，也可以自定义基因集，注意转为对应的形式。human_KEGG_Set 为基因集合的列表形式。

二 GSVA分析

1, GSVA分析

数据准备好后，加载GSVA包，一个gsva函数就可以得到GSVA的结果了。

library(GSVA)
gsva.kegg <- gsva(expr, gset.idx.list = human_KEGG_Set, 
             kcdf="Gaussian",
             method = "gsva",
             parallel.sz=1)
head(gsva.kegg)

                                             Epi     Myeloid Fibroblast           T        Endo         un
KEGG_ABC_TRANSPORTERS                           -0.30909556 -0.38967397 -0.4080869 -0.36390157 -0.23292000 -0.2954922
KEGG_ACUTE_MYELOID_LEUKEMIA                     -0.51620884  0.03252670 -0.4041246  0.15837659  0.27282853 -0.3265242
KEGG_ADHERENS_JUNCTION                          -0.25225184 -0.24621834 -0.1800216 -0.35910806  0.12149392 -0.2252523
KEGG_ADIPOCYTOKINE_SIGNALING_PATHWAY            -0.49448870 -0.13413340 -0.4633362 -0.24019340 -0.12668157 -0.3430883
KEGG_ALANINE_ASPARTATE_AND_GLUTAMATE_METABOLISM  0.03912373 -0.25381877 -0.3541902 -0.03168169 -0.38696763  0.3222462
KEGG_ALDOSTERONE_REGULATED_SODIUM_REABSORPTION  -0.28822490  0.07870691 -0.2787129  0.05187873 -0.03042265  0.1007009

行为目标基因集，列为celltype ，数值为gsva分数。

这里需要注意，如果输入矩阵为log转化后的连续表达矩阵指则设置kcdf参数为"Gaussian"，如果是counts矩阵则设置kcdf为"Poisson"。

2, 绘制热图

以结果的前20个绘制示例热图，可以自选择重点的通路

library(pheatmap)
pheatmap(gsva.kegg[1:20,], show_colnames = T, 
         scale = "row",angle_col = "45",
         cluster_row = T,cluster_col = T,
         color = colorRampPalette(c("navy", "white", "firebrick3"))(50))

3, limma差异分析

和正常的转录组差异分析一样，构建分组信息 以及 比较矩阵，然后使用limma进行差异分析。

此处分析注释出来的5种celltype 和注释为unknown之间的差异GSVA结果。

library(limma)
# 构建分组文件
group_list <- data.frame(celltype = colnames(gsva.kegg), group = c(rep("con", 5), rep("case", 1)))

design <- model.matrix(~ 0 + factor(group_list$group))
colnames(design) <- levels(factor(group_list$group))
rownames(design) <- colnames(gsva.kegg)

# 构建差异比较矩阵
contrast.matrix <- makeContrasts(con-case, levels = design)

# 差异分析，case vs. con
fit <- lmFit(gsva.kegg, design)
fit2 <- contrasts.fit(fit, contrast.matrix)
fit2 <- eBayes(fit2)
diff <- topTable(fit2, coef = 1, n = Inf, adjust.method = "BH", sort.by = "P")
head(diff)

                                                       logFC      AveExpr         t    P.Value adj.P.Val         B
KEGG_RIBOSOME                                             -0.7314315 -0.240532302 -2.289899 0.03492937 0.9966699 -4.382757
KEGG_ALANINE_ASPARTATE_AND_GLUTAMATE_METABOLISM           -0.5197531 -0.110881401 -2.100093 0.05077169 0.9966699 -4.418969
KEGG_MATURITY_ONSET_DIABETES_OF_THE_YOUNG                 -0.5827932 -0.078191940 -1.986597 0.06314026 0.9966699 -4.440295
KEGG_PHENYLALANINE_METABOLISM                             -0.4828199 -0.125947781 -1.888280 0.07597847 0.9966699 -4.458485
KEGG_ARRHYTHMOGENIC_RIGHT_VENTRICULAR_CARDIOMYOPATHY_ARVC -0.4559099 -0.080461313 -1.826270 0.08522283 0.9966699 -4.469792
KEGG_FOLATE_BIOSYNTHESIS                                  -0.4870146 -0.004848985 -1.811745 0.08752649 0.9966699 -4.472419

4, GSVA差异结果可视化

文章中常见的为双向的柱形图，需要先进行一些设置：

（1）根据logFC 和 adj.P.Val 参数设置，确定上调，下调。为展示结果，以下参数较离谱，无参考价值。

实际项目多设置为logFC > 0 & adj.P.Val < 0.05；

#设置分组
diff$group <- ifelse( diff$logFC > 0 & diff$P.Value < 0.3 ,"up" ,
  ifelse(diff$logFC < 0 & diff$P.Value < 0.3 ,"down","noSig")
)

（2）设置label的位置

本示例中过滤掉了不显著的通路，在filter行首添加#注释掉即不进行过滤

diff2 <- diff %>% 
  mutate(hjust2 = ifelse(t>0,1,0)) %>% 
  mutate(nudge_y = ifelse(t>0,-0.1,0.1)) %>% 
  filter(group != "noSig") %>% #注释掉该行即可
  arrange(t) %>% 
  rownames_to_column("ID")

（3）通过factor设置展示顺序

diff2$ID <- factor(diff2$ID, levels = diff2$ID)
limt = max(abs(diff2$t))

（4）ggplot2 可视化展示

ggplot(diff2, aes(ID, t,fill=group)) + 
  geom_bar(stat = 'identity',alpha = 0.7) + 
  scale_fill_manual(breaks=c("down","up"), #设置颜色
                    values = c("#008020","#08519C"))+
  geom_text(data = diff2, aes(label = diff2$ID, #添加通路标签
                              y = diff2$nudge_y),
            nudge_x =0,nudge_y =0,hjust =diff2$hjust,
            size = 3)+ #设置字体大小
  labs(x = "KEGG pathways", #设置标题 和 坐标
       y=paste0("t value of GSVA score\n","celltype-unknown"),
       title = "GSVA")+
  scale_y_continuous(limits=c(-limt,limt))+
  coord_flip() + 
  theme_bw() + #去除背景色
  theme(panel.grid =element_blank(), #主题微调
        panel.border = element_rect(size = 0.6),
        plot.title = element_text(hjust = 0.5,size = 18),
        axis.text.y = element_blank(),
        axis.title = element_text(hjust = 0.5,size = 18),
        axis.line = element_blank(),
        axis.ticks.y = element_blank(),
        legend.position = "none" #去掉legend
  )

三 GSVA 样本分组

1, 表达量文件

如果是按照样本分组的话就无需计算每个celltype的表达量均值，直接使用每个细胞的表达量；

expr2 <- as.matrix(sub@assays$RNA@data)
gsva.kegg2 <- gsva(expr2, gset.idx.list = keggSet, kcdf="Gaussian",method = "gsva",
             parallel.sz=1)

2, 分组文件

分组文件因为是每个barcode的粒度，在metadata中构建分组列信息

#之前定义过分组信息
sce2@meta.data$group <- ifelse( grepl("MET",sce2@meta.data$sample ) ,"MET" ,"PT" )
group_list2 <- sce2@meta.data[,c("group")]

#标准limma分析
design <- model.matrix(~0+factor(group_list2))
colnames(design)=levels(factor(group_list2))
rownames(design)=colnames(expr)

contrast.matrix<-makeContrasts(contrasts = "MET-PT",levels = design)

fit <- lmFit(gsva.kegg2,design)
fit2 <- contrasts.fit(fit, contrast.matrix)
fit2 <- eBayes(fit2)
diff2 <- topTable(fit2, coef = 1, n = Inf, adjust.method = "BH", sort.by = "P")