下面是他對我們b站轉(zhuǎn)錄組視頻課程的詳細(xì)筆記 
承接上節(jié)：RNA-seq入門實(shí)戰(zhàn)（四）：差異分析前的準(zhǔn)備——數(shù)據(jù)檢查，以及 RNA-seq入門實(shí)戰(zhàn)（五）：差異分析——DESeq2 edgeR limma的使用與比較 本節(jié)概覽：
1.GSVA簡單介紹
2.基因集的下載讀取: 手動與msigdbr包下載
3.GSVA的運(yùn)行
4.limma差異分析
5.GSVA結(jié)果可視化：熱圖,、火山圖,、發(fā)散條形圖/柱形偏差圖
1. GSVA簡單介紹官方文檔：GSVA: gene set variation analysis (bioconductor.org)不錯的一篇文章：GSVA的使用 - raisok 定義基因集變異分析(GSVA)是一種特殊類型的基因集富集方法，通過對分析的功能單元進(jìn)行概念上簡單但功能強(qiáng)大的改變——從基因到基因集,，從而實(shí)現(xiàn)對分子數(shù)據(jù)的路徑中心分析,。
簡單來說，就是將分析對象由基因換成了基因集,，進(jìn)行基因集（通路）級別的差異分析,。 原理和作用通過將基因在不同樣品間的表達(dá)量矩陣轉(zhuǎn)化成基因集在樣品間的表達(dá)量矩陣，從而來評估不同的通路在不同樣品間是否富集,。其實(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")source('FUNCTIONS.R')load(file = '1.counts.Rdata') #包含 tpm counts group_list gldir.create("6.GSVA"); setwd("6.GSVA")
2.下載GSVA分析所需的基因集GSVA分析常用MSigDB數(shù)據(jù)庫中基因集，也可以自定義基因集進(jìn)行分析,。
MSigDB數(shù)據(jù)庫目前有H和C1-C8九個定義的基因集,，下面示范下載包含KEGG信息的C2與包含GO信息的C5基因集的兩種方式——手動下載讀取或msigdbr包下載提取。 
2.1 手動下載讀取下載地址：Downloads (gsea-msigdb.org)進(jìn)入gsea-msigdb官網(wǎng)后可能還需要登錄郵箱,，然后找到需要下載的基因集下載,，下載后將gmt文件放入指定文件夾，將其路徑讀取進(jìn)入R即可,。不過需要注意的是這里的基因集默認(rèn)都是人類的,，如果是分析小鼠或其他物種最好采用MigDB包下載
 #### 對 MigDB( Molecular Signatures Database)中的基因集做GSVA分析 ###### 用手動下載基因集做GSVA分析d <- 'C:/Users/Lenovo/Desktop/test/gmt' #存放gmt文件的路徑gmtfs <- list.files(d,pattern = 'symbols.gmt') # 路徑下所有結(jié)尾為symbols.gmt文件gmtfskegg_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分析##KEGGKEGG_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分組
##GOGO_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)行多核計算,。
下面選取基因集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)行名變成了基因集通路名,，每個樣品都會有對應(yīng)通路的GSVA評分： 
4. limma差異分析得到GSVA評分的矩陣后,，我們需要利用limma包進(jìn)行pathway通路的差異分析，與之前介紹的基因差異分析流程類似,，但不需要進(jìn)行 limma-trend 或 voom的步驟 #### 進(jìn)行l(wèi)imma差異處理 ######設(shè)定 實(shí)驗(yàn)組exp / 對照組ctrglexp="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)計檢驗(yàn)efit <- eBayes(fit2) #修正
summary(decideTests(efit,lfc=1, p.value=0.05)) #統(tǒng)計查看差異結(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 熱圖#### 對GSVA的差異分析結(jié)果進(jìn)行熱圖可視化 #### ##設(shè)置篩選閾值padj_cutoff=0.05log2FC_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ǎ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ì)繪圖過程見發(fā)散條形圖/柱形偏差圖 - 簡書 (jianshu.com) #### 發(fā)散條形圖繪制 ####library(tidyverse) # ggplot2 stringer dplyr tidyr readr purrr tibble forcatslibrary(ggthemes)library(ggprism)p_cutoff=0.001degs <- 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è)為組-1dat_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的使用 - raisokGSEA和GSVA，再也不用去下載gmt文件咯 - 簡書 (jianshu.com) 發(fā)散條形圖/柱形偏差圖 - 簡書 (jianshu.com) GitHub - jmzeng1314/GEO 【生信技能樹】轉(zhuǎn)錄組測序數(shù)據(jù)分析_嗶哩嗶哩_bilibili 【生信技能樹】GEO數(shù)據(jù)庫挖掘_嗶哩嗶哩_bilibili
|
