本節(jié)概覽：
1.GSVA簡(jiǎn)單介紹
2.基因集的下載讀取: 手動(dòng)與msigdbr包下載
3.GSVA的運(yùn)行
4.limma差異分析
5.GSVA結(jié)果可視化：熱圖、火山圖、發(fā)散條形圖/柱形偏差圖

1. GSVA簡(jiǎn)單介紹

官方文檔：GSVA: gene set variation analysis (bioconductor.org)
不錯(cuò)的一篇文章：GSVA的使用 - raisok

定義
基因集變異分析(GSVA)是一種特殊類型的基因集富集方法，通過(guò)對(duì)分析的功能單元進(jìn)行概念上簡(jiǎn)單但功能強(qiáng)大的改變——從基因到基因集，從而實(shí)現(xiàn)對(duì)分子數(shù)據(jù)的路徑中心分析。
簡(jiǎn)單來(lái)說(shuō)，就是將分析對(duì)象由基因換成了基因集，進(jìn)行基因集（通路）級(jí)別的差異分析。
原理和作用
通過(guò)將基因在不同樣品間的表達(dá)量矩陣轉(zhuǎn)化成基因集在樣品間的表達(dá)量矩陣，從而來(lái)評(píng)估不同的通路在不同樣品間是否富集。其實(shí)就是研究這些感興趣的基因集在不同樣品間的差異，或者尋找比較重要的基因集，作為一種分析方法，主要是是為了從生物信息學(xué)的角度去解釋導(dǎo)致表型差異的原因。

GSVA官方定義.png

分析前準(zhǔn)備：

rm(list = ls())
options(stringsAsFactors = F)
library(tidyverse)
library(clusterProfiler)
library(msigdbr)  #install.packages("msigdbr")
library(GSVA) 
library(GSEABase)
library(pheatmap)
library(limma)
library(BiocParallel)

setwd("C:/Users/Lenovo/Desktop/test")
load(file = '1.counts.Rdata') #包含 tpm counts group_list gl
dir.create("6.GSVA"); setwd("6.GSVA")

2.下載GSVA分析所需的基因集

GSVA分析常用MSigDB數(shù)據(jù)庫(kù)中基因集，也可以自定義基因集進(jìn)行分析。
MSigDB數(shù)據(jù)庫(kù)目前有H和C1-C8九個(gè)定義的基因集，下面示范下載包含KEGG信息的C2與包含GO信息的C5基因集的兩種方式——手動(dòng)下載讀取或msigdbr包下載提取。

2.1 手動(dòng)下載讀取

下載地址：Downloads (gsea-msigdb.org)
進(jìn)入gsea-msigdb官網(wǎng)后可能還需要登錄郵箱，然后找到需要下載的基因集下載，下載后將gmt文件放入指定文件夾，將其路徑讀取進(jìn)入R即可。不過(guò)需要注意的是這里的基因集默認(rèn)都是人類的，如果是分析小鼠或其他物種最好采用MigDB包下載

#### 對(duì) MigDB( Molecular Signatures Database)中的基因集做GSVA分析  ####
## 用手動(dòng)下載基因集做GSVA分析
d <- 'C:/Users/Lenovo/Desktop/test/gmt' #存放gmt文件的路徑
gmtfs <- list.files(d,pattern = 'symbols.gmt')  # 路徑下所有結(jié)尾為symbols.gmt文件
gmtfs
kegg_list <- getGmt(file.path(d,gmtfs[1])) 
go_list <- getGmt(file.path(d,gmtfs[2]))

2.2 msigdbr包下載讀取

使用msigdbr包可以直接在R里下載C2和C5基因集，并提取相關(guān)信息做成list。
msigdbr下載數(shù)據(jù)可以指定物種，用msigdbr_species() 和 msigdbr_collections()可以查看支持的物種和基因集類別。
以下參考：GSEA和GSVA，再也不用去下載gmt文件咯

## msigdbr包提取下載 一般嘗試KEGG和GO做GSVA分析
##KEGG
KEGG_df_all <-  msigdbr(species = "Mus musculus", # Homo sapiens or Mus musculus
                        category = "C2",
                        subcategory = "CP:KEGG") 
KEGG_df <- dplyr::select(KEGG_df_all,gs_name,gs_exact_source,gene_symbol)
kegg_list <- split(KEGG_df$gene_symbol, KEGG_df$gs_name) ##按照gs_name給gene_symbol分組

##GO
GO_df_all <- msigdbr(species = "Mus musculus",
                     category = "C5")  
GO_df <- dplyr::select(GO_df_all, gs_name, gene_symbol, gs_exact_source, gs_subcat)
GO_df <- GO_df[GO_df$gs_subcat!="HPO",]
go_list <- split(GO_df$gene_symbol, GO_df$gs_name) ##按照gs_name給gene_symbol分組

3. GSVA的運(yùn)行

使用GSVA需要輸入基因表達(dá)矩陣和基因集。
基因集即為我們上一步所得list；基因表達(dá)矩陣可以使用logCPM、logRPKM、logTPM（GSVA參數(shù)kcdf選擇"Gaussian"，默認(rèn)）或counts數(shù)據(jù)（參數(shù)kcdf選擇"Poisson"）。
GSVA還支持BiocParallel，可設(shè)置參數(shù)parallel.sz進(jìn)行多核計(jì)算。
下面選取基因集go_list和logTPM數(shù)據(jù)進(jìn)行示范


####  GSVA  ####
#GSVA算法需要處理logCPM, logRPKM,logTPM數(shù)據(jù)或counts數(shù)據(jù)的矩陣####
#dat <- as.matrix(counts)
#dat <- as.matrix(log2(edgeR::cpm(counts))+1)
dat <- as.matrix(log2(tpm+1))

geneset <- go_list

gsva_mat <- gsva(expr=dat, 
               gset.idx.list=geneset, 
               kcdf="Gaussian" ,#"Gaussian" for logCPM,logRPKM,logTPM, "Poisson" for counts
               verbose=T, 
               parallel.sz = parallel::detectCores())#調(diào)用所有核

write.csv(gsva_mat,"gsva_go_matrix.csv")

運(yùn)行完GSVA后gsva_mat內(nèi)容如下，可以發(fā)現(xiàn)行名變成了基因集通路名，每個(gè)樣品都會(huì)有對(duì)應(yīng)通路的GSVA評(píng)分：

4. limma差異分析

得到GSVA評(píng)分的矩陣后，我們需要利用limma包進(jìn)行pathway通路的差異分析，與之前介紹的基因差異分析流程類似，但不需要進(jìn)行 limma-trend 或 voom的步驟


#### 進(jìn)行l(wèi)imma差異處理 ####
##設(shè)定 實(shí)驗(yàn)組exp / 對(duì)照組ctr
gl
exp="primed"
ctr="naive"

design <- model.matrix(~0+factor(group_list))
colnames(design) <- levels(factor(group_list))
rownames(design) <- colnames(gsva_mat)
contrast.matrix <- makeContrasts(contrasts=paste0(exp,'-',ctr),  #"exp/ctrl"
                                 levels = design)

fit1 <- lmFit(gsva_mat,design)                 #擬合模型
fit2 <- contrasts.fit(fit1, contrast.matrix) #統(tǒng)計(jì)檢驗(yàn)
efit <- eBayes(fit2)                         #修正

summary(decideTests(efit,lfc=1, p.value=0.05)) #統(tǒng)計(jì)查看差異結(jié)果
tempOutput <- topTable(efit, coef=paste0(exp,'-',ctr), n=Inf)
degs <- na.omit(tempOutput) 
write.csv(degs,"gsva_go_degs.results.csv")

5. GSVA結(jié)果可視化：熱圖、火山圖、發(fā)散條形圖/柱形偏差圖

與常規(guī)差異分析結(jié)果展示類似，GSVA結(jié)果可視化一般也用熱圖、火山圖展示

5.1 熱圖

#### 對(duì)GSVA的差異分析結(jié)果進(jìn)行熱圖可視化 #### 
##設(shè)置篩選閾值
padj_cutoff=0.05
log2FC_cutoff=log2(2)

keep <- rownames(degs[degs$adj.P.Val < padj_cutoff & abs(degs$logFC)>log2FC_cutoff, ])
length(keep)
dat <- gsva_mat[keep[1:50],] #選取前50進(jìn)行展示

pheatmap::pheatmap(dat, 
                   fontsize_row = 8,
                   height = 10,
                   width=18,
                   annotation_col = gl,
                   show_colnames = F,
                   show_rownames = T,
                   filename = paste0('GSVA_go_heatmap.pdf'))

5.2 火山圖

degs$significance  <- as.factor(ifelse(degs$adj.P.Val < padj_cutoff & abs(degs$logFC) > log2FC_cutoff,
                                           ifelse(degs$logFC > log2FC_cutoff ,'UP','DOWN'),'NOT'))

this_title <- paste0(' Up :  ',nrow(degs[degs$significance =='UP',]) ,
                     '\n Down : ',nrow(degs[degs$significance =='DOWN',]),
                     '\n adj.P.Val <= ',padj_cutoff,
                     '\n FoldChange >= ',round(2^log2FC_cutoff,3))


g <- ggplot(data=degs, 
            aes(x=logFC, y=-log10(adj.P.Val),
                color=significance)) +
  #點(diǎn)和背景
  geom_point(alpha=0.4, size=1) +
  theme_classic()+ #無(wú)網(wǎng)格線
  #坐標(biāo)軸
  xlab("log2 ( FoldChange )") + 
  ylab("-log10 ( adj.P.Val )") +
  #標(biāo)題文本
  ggtitle( this_title ) +
  #分區(qū)顏色                   
  scale_colour_manual(values = c('blue','grey','red'))+ 
  #輔助線
  geom_vline(xintercept = c(-log2FC_cutoff,log2FC_cutoff),lty=4,col="grey",lwd=0.8) +
  geom_hline(yintercept = -log10(padj_cutoff),lty=4,col="grey",lwd=0.8) +
  #圖例標(biāo)題間距等設(shè)置
  theme(plot.title = element_text(hjust = 0.5), 
        plot.margin=unit(c(2,2,2,2),'lines'), #上右下左
        legend.title = element_blank(),
        legend.position="right")

ggsave(g,filename = 'GSVA_go_volcano_padj.pdf',width =8,height =7.5)

5.3 發(fā)散條形圖/柱形偏差圖

為了更好展示繪制發(fā)散條形圖/柱形偏差圖，此處用的是KEGG的gsva差異分析結(jié)果，展示通路的上下調(diào)及pvalue信息（也可以是t值或padj值等），詳細(xì)繪圖過(guò)程見發(fā)散條形圖/柱形偏差圖 - 簡(jiǎn)書 (jianshu.com)

#### 發(fā)散條形圖繪制 ####
library(tidyverse)  # ggplot2 stringer dplyr tidyr readr purrr  tibble forcats
library(ggthemes)
library(ggprism)
p_cutoff=0.001
degs <- gsva_kegg_degs  #載入gsva的差異分析結(jié)果
Diff <- rbind(subset(degs,logFC>0)[1:20,], subset(degs,logFC<0)[1:20,]) #選擇上下調(diào)前20通路     
dat_plot <- data.frame(id  = row.names(Diff),
                       p   = Diff$P.Value,
                       lgfc= Diff$logFC)
dat_plot$group <- ifelse(dat_plot$lgfc>0 ,1,-1)    # 將上調(diào)設(shè)為組1，下調(diào)設(shè)為組-1
dat_plot$lg_p <- -log10(dat_plot$p)*dat_plot$group # 將上調(diào)-log10p設(shè)置為正，下調(diào)-log10p設(shè)置為負(fù)
dat_plot$id <- str_replace(dat_plot$id, "KEGG_","");dat_plot$id[1:10]
dat_plot$threshold <- factor(ifelse(abs(dat_plot$p) <= p_cutoff,
                                   ifelse(dat_plot$lgfc >0 ,'Up','Down'),'Not'),
                            levels=c('Up','Down','Not'))

dat_plot <- dat_plot %>% arrange(lg_p)
dat_plot$id <- factor(dat_plot$id,levels = dat_plot$id)

## 設(shè)置不同標(biāo)簽數(shù)量
low1 <- dat_plot %>% filter(lg_p < log10(p_cutoff)) %>% nrow()
low0 <- dat_plot %>% filter(lg_p < 0) %>% nrow()
high0 <- dat_plot %>% filter(lg_p < -log10(p_cutoff)) %>% nrow()
high1 <- nrow(dat_plot)

p <- ggplot(data = dat_plot,aes(x = id, y = lg_p, 
                                fill = threshold)) +
  geom_col()+
  coord_flip() + 
  scale_fill_manual(values = c('Up'= '#36638a','Not'='#cccccc','Down'='#7bcd7b')) +
  geom_hline(yintercept = c(-log10(p_cutoff),log10(p_cutoff)),color = 'white',size = 0.5,lty='dashed') +
  xlab('') + 
  ylab('-log10(P.Value) of GSVA score') + 
guides(fill="none")+
  theme_prism(border = T) +
  theme(
    axis.text.y = element_blank(),
    axis.ticks.y = element_blank()
  ) +
  geom_text(data = dat_plot[1:low1,],aes(x = id,y = 0.1,label = id),
              hjust = 0,color = 'black') + #黑色標(biāo)簽
  geom_text(data = dat_plot[(low1 +1):low0,],aes(x = id,y = 0.1,label = id),
            hjust = 0,color = 'grey') + # 灰色標(biāo)簽
  geom_text(data = dat_plot[(low0 + 1):high0,],aes(x = id,y = -0.1,label = id),
            hjust = 1,color = 'grey') + # 灰色標(biāo)簽
  geom_text(data = dat_plot[(high0 +1):high1,],aes(x = id,y = -0.1,label = id),
            hjust = 1,color = 'black') # 黑色標(biāo)簽

ggsave("GSVA_barplot_pvalue.pdf",p,width = 15,height  = 15)

GSVA就到這了，下一節(jié)學(xué)習(xí)一下蛋白互作網(wǎng)絡(luò)PPI吧

參考資料
GSVA: gene set variation analysis (bioconductor.org)
GSVA的使用 - raisok
GSEA和GSVA，再也不用去下載gmt文件咯 - 簡(jiǎn)書 (jianshu.com)
發(fā)散條形圖/柱形偏差圖 - 簡(jiǎn)書 (jianshu.com)
GitHub - jmzeng1314/GEO
【生信技能樹】轉(zhuǎn)錄組測(cè)序數(shù)據(jù)分析_嗶哩嗶哩_bilibili
【生信技能樹】GEO數(shù)據(jù)庫(kù)挖掘_嗶哩嗶哩_bilibili

色偷偷精品伊人,欧洲久久精品,欧美综合婷婷骚逼,国产AV主播,国产最新探花在线,九色在线视频一区,伊人大交九欧美,1769亚洲,黄色成人av

RNA-seq入門實(shí)戰(zhàn)（八）：GSVA——基因集變異分析

RNA-seq入門實(shí)戰(zhàn)（八）：GSVA——基因集變異分析

1. GSVA簡(jiǎn)單介紹

2.下載GSVA分析所需的基因集

2.1 手動(dòng)下載讀取

2.2 msigdbr包下載讀取

3. GSVA的運(yùn)行

4. limma差異分析

5. GSVA結(jié)果可視化：熱圖、火山圖、發(fā)散條形圖/柱形偏差圖

5.1 熱圖

5.2 火山圖

5.3 發(fā)散條形圖/柱形偏差圖

相關(guān)閱讀更多精彩內(nèi)容

友情鏈接更多精彩內(nèi)容

色偷偷精品伊人,欧洲久久精品,欧美综合婷婷骚逼,国产AV主播,国产最新探花在线,九色在线视频一区,伊人大交九 欧美,1769亚洲,黄色成人av

RNA-seq入門實(shí)戰(zhàn)（八）：GSVA——基因集變異分析

1. GSVA簡(jiǎn)單介紹

2.下載GSVA分析所需的基因集

2.1 手動(dòng)下載讀取

2.2 msigdbr包下載讀取

3. GSVA的運(yùn)行

4. limma差異分析

5. GSVA結(jié)果可視化：熱圖、火山圖、發(fā)散條形圖/柱形偏差圖

5.1 熱圖

5.2 火山圖

5.3 發(fā)散條形圖/柱形偏差圖

相關(guān)閱讀更多精彩內(nèi)容

友情鏈接更多精彩內(nèi)容

色偷偷精品伊人,欧洲久久精品,欧美综合婷婷骚逼,国产AV主播,国产最新探花在线,九色在线视频一区,伊人大交九欧美,1769亚洲,黄色成人av

5. GSVA結(jié)果可視化：熱圖、火山圖、發(fā)散條形圖/柱形偏差圖