与传统的Over Representation Analysis: the up- or down-regulated DEGs和Gene Set Enrichment Analysis: All the ranked genes相比，Gene Set Variation Analysis分析直接使用gene expression或RNA-seq profile matrix计算基因在通路中的得分，并且这种计算可以基于单样本处理。

  • GSVA方法图解
    数据分析：基于GSVA的通路富集分析 - 图1

输入表达谱数据和基因集数据。计算其得分。

The GSVA package implements a non-parametric unsupervised method, called Gene Set Variation Analysis (GSVA), for assessing gene set enrichment (GSE) in gene expression microarray and RNA- seq data. In contrast to most GSE methods, GSVA performs a change in coordinate systems, transforming the data from a gene by sample matrix to a gene set by sample matrix. Thereby allowing for the evaluation of pathway enrichment for each sample. This transformation is done without the use of a phenotype, thus facilitating very powerful and open-ended analyses in a now pathway centric manner.

加载R包

  1. knitr::opts_chunk$set(warning = F, message = F)
  2. library(dplyr)
  3. library(tibble)
  4. library(ggplot2)
  5. library(data.table)
  6. library(GSVA)
  8. rm(list = ls())
  9. options(stringsAsFactors = F)
  10. options(future.globals.maxSize = 1000 * 1024^2)
  12. grp <- c("S1", "S2")
  13. grp.col <- c("#6C326C", "#77A2D1")

load data

  • gene Expression DESeq2 object
  1. phen <- fread("../../Result/phenotype/phenotype_cluster.csv")
  2. geneExp <- fread("../../Result/profile/geneExp_filter.tsv")
  3. # 处理数据成整型
  4. geneExp <- geneExp %>% column_to_rownames("V1") %>% as.matrix()
  5. geneExp <- round(geneExp) %>% data.frame() %>% rownames_to_column("V1")
  7. pathways_hallmark_kegg <- fgsea::gmtPathways("../../Result/GeneID/msigdb.v7.1.symbols_KEGG.gmt")
  8. pathways_hallmark_GO <- fgsea::gmtPathways("../../Result/GeneID/msigdb.v7.1.symbols_GO.gmt")

Curation Function

  1. get_ExprSet <- function(x=phen, 
  2.                         y=geneExp){
  3.   # x=phen
  4.   # y=geneExp
  6.   sid <- intersect(x$Barcode, colnames(y))
  7.   # phenotype
  8.   phe <- x %>% filter(Barcode%in%sid) %>%
  9.     mutate(Cluster=factor(as.character(Cluster))) %>%
  10.     column_to_rownames("Barcode") 
  12.   # profile 
  13.   prf <- y %>% column_to_rownames("V1") %>%
  14.     dplyr::select(all_of(sid))
  16.   # determine the right order between profile and phenotype 
  17.   for(i in 1:ncol(prf)){ 
  18.     if (!(colnames(prf)[i] == rownames(phe)[i])) {
  19.       stop(paste0(i, " Wrong"))
  20.     }
  21.   }
  22.   require(convert)
  23.   exprs <- as.matrix(prf)
  24.   adf <-  new("AnnotatedDataFrame", data=phe)
  25.   experimentData <- new("MIAME",
  26.         name="Hua Zou", lab="Hua Lab",
  27.         contact="zouhua1@outlook.com",
  28.         title="KRIC Experiment",
  29.         abstract="The ExpressionSet",
  30.         url="www.zouhua.top",
  31.         other=list(notes="Created from text files"))
  32.   expressionSet <- new("ExpressionSet", exprs=exprs,
  33.                        phenoData=adf, 
  34.                        experimentData=experimentData)
  36.   return(expressionSet)
  37. }
  39. get_score <- function(dataset=get_ExprSet(x=phen, y=geneExp),
  40.                       geneset=pathways_hallmark_kegg,
  41.                       methods="ssgsea"){
  43.   # dataset=get_ExprSet(x=phen, y=geneExp)
  44.   # geneset=pathways_hallmark_kegg
  45.   # methods="ssgsea"  
  47.   dat_fit <- gsva(expr = dataset, 
  48.               gset.idx.list = geneset,
  49.               method = methods, 
  50.               min.sz = 5, 
  51.               max.sz = 500,
  52.               kcdf = "Poisson")
  53.   res <- exprs(dat_fit) %>% t() %>% 
  54.     data.frame() %>%
  55.     rownames_to_column("SampleID")
  57.   return(res)
  58. }
  60. # Differential methylation/copyNumber/methylation
  61. get_limma <- function(dataset=methylation_set,
  62.                       group_col=grp,
  63.                       tag="methylation",
  64.                       fc=1,
  65.                       pval=0.05){
  67.   # dataset=methylation_set
  68.   # group_col=grp
  69.   # tag="methylation"  
  70.   # fc=1
  71.   # pval=0.05
  73.   pheno <- pData(dataset) 
  75.   if(tag == "methylation"){
  76.     # transform the beta value into M values via lumi package
  77.     require(lumi)
  78.     edata <- beta2m(exprs(dataset))    
  79.   }else{
  80.     edata <- exprs(dataset)
  81.   }
  84.   require(limma)
  85.   design <- model.matrix(~0 + pheno$Cluster)
  86.   rownames(design) <- rownames(pheno)
  87.   colnames(design) <- group_col
  88.   exprSet <- edata  
  90.   # show distribution
  91.   boxplot(exprSet)
  92.   plotDensities(exprSet) 
  94.   # linear fitting 
  95.   fit <- lmFit(exprSet, design, method = 'ls')
  96.   contrast <- makeContrasts("S1-S2", levels = design) 
  98.   # eBayes
  99.   fit2 <- contrasts.fit(fit, contrast)
  100.   fit2 <- eBayes(fit2)
  102.   qqt(fit2$t, df = fit2$df.prior + fit2$df.residual, pch = 16, cex = 0.2)
  103.   abline(0, 1)  
  105.   # differential features
  106.   diff_feature <- topTable(fit2, number = Inf, adjust.method = 'BH') %>%
  107.     rownames_to_column("Feature") 
  109.   # prof[rownames(prof)%in%"EVI2A", ] %>% data.frame() %>% setNames("EVI2A") %>%
  110.   #   rownames_to_column("SampleID") %>%
  111.   #   inner_join(pheno %>%rownames_to_column("SampleID"), by = "SampleID") -> a1
  112.   # wilcox.test(EVI2A~Cluster, data = a1)
  113.   # ggplot(a1, aes(x=Cluster, y=EVI2A))+geom_boxplot() 
  115.   diff_feature[which(diff_feature$logFC >= fc & 
  116.                        diff_feature$adj.P.Val < pval), "Enrichment"] <- group_col[1]
  117.   diff_feature[which(diff_feature$logFC <= -fc & 
  118.                        diff_feature$adj.P.Val < pval), "Enrichment"] <- group_col[2]
  119.   diff_feature[which(abs(diff_feature$logFC) < fc |
  120.                        diff_feature$adj.P.Val >= pval), "Enrichment"] <- "Nonsignif"
  122.   diff_res <- diff_feature %>% 
  123.     setNames(c("Feature", "log2FoldChange", "baseMean", "t", 
  124.                "pvalue", "padj", "B", "Enrichment")) %>% 
  125.     dplyr::select(Feature, everything()) %>%
  126.     arrange(padj, log2FoldChange) %>%
  127.     column_to_rownames("Feature")
  129.   res <- list(fit=fit2, diff_res=diff_res)
  131.   return(res)
  132. }

gsva score

kegg_gsva <- get_score(dataset=get_ExprSet(x=phen, y=geneExp),
                       geneset=pathways_hallmark_kegg,
                       methods="gsva")

if(!dir.exists("../../Result/pathway/")){
  dir.create("../../Result/pathway/", recursive = T)
}
write.csv(kegg_gsva, "../../Result/pathway/KEGG_gsva.csv", row.names = F)

# Differential Test
kegg_gsva_exprset <-  get_ExprSet(x=phen, y=kegg_gsva %>% column_to_rownames("SampleID") %>% t() %>% data.frame() %>%rownames_to_column("V1"))  
Diff_gsva <- get_limma(dataset =kegg_gsva_exprset, tag = "gsva")
DT::datatable(Diff_gsva$diff_res)

version

sessionInfo()

R version 4.0.2 (2020-06-22)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 19042)

Matrix products: default

locale:
[1] LC_COLLATE=English_United States.1252  LC_CTYPE=English_United States.1252   
[3] LC_MONETARY=English_United States.1252 LC_NUMERIC=C                          
[5] LC_TIME=English_United States.1252    
system code page: 936

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] GSVA_1.36.3       data.table_1.13.6 ggplot2_3.3.3     tibble_3.0.3      dplyr_1.0.2      

loaded via a namespace (and not attached):
  [1] fgsea_1.14.0                colorspace_1.4-1            ellipsis_0.3.1             
  [4] ggridges_0.5.2              qvalue_2.20.0               XVector_0.28.0             
  [7] GenomicRanges_1.40.0        rstudioapi_0.11             farver_2.0.3               
 [10] urltools_1.7.3              graphlayouts_0.7.1          ggrepel_0.9.1.9999         
 [13] DT_0.16                     bit64_0.9-7                 AnnotationDbi_1.50.3       
 [16] scatterpie_0.1.5            xml2_1.3.2                  splines_4.0.2              
 [19] GOSemSim_2.14.2             knitr_1.30                  shinythemes_1.1.2          
 [22] polyclip_1.10-0             jsonlite_1.7.1              annotate_1.66.0            
 [25] GO.db_3.11.4                graph_1.66.0                ggforce_0.3.2              
 [28] shiny_1.5.0                 BiocManager_1.30.10         compiler_4.0.2             
 [31] httr_1.4.2                  rvcheck_0.1.8               Matrix_1.3-2               
 [34] fastmap_1.0.1               later_1.1.0.1               tweenr_1.0.1               
 [37] htmltools_0.5.0             prettyunits_1.1.1           tools_4.0.2                
 [40] igraph_1.2.5                gtable_0.3.0                glue_1.4.2                 
 [43] GenomeInfoDbData_1.2.3      reshape2_1.4.4              DO.db_2.9                  
 [46] fastmatch_1.1-0             Rcpp_1.0.5                  enrichplot_1.8.1           
 [49] Biobase_2.48.0              vctrs_0.3.4                 ggraph_2.0.3               
 [52] xfun_0.19                   stringr_1.4.0               mime_0.9                   
 [55] lifecycle_0.2.0             XML_3.99-0.5                DOSE_3.14.0                
 [58] europepmc_0.4               MASS_7.3-53                 zlibbioc_1.34.0            
 [61] scales_1.1.1                tidygraph_1.2.0             hms_0.5.3                  
 [64] promises_1.1.1              parallel_4.0.2              SummarizedExperiment_1.18.2
 [67] RColorBrewer_1.1-2          yaml_2.2.1                  memoise_1.1.0              
 [70] gridExtra_2.3               triebeard_0.3.0             stringi_1.5.3              
 [73] RSQLite_2.2.0               S4Vectors_0.26.1            BiocGenerics_0.34.0        
 [76] BiocParallel_1.22.0         GenomeInfoDb_1.24.2         rlang_0.4.7                
 [79] pkgconfig_2.0.3             bitops_1.0-6                matrixStats_0.57.0         
 [82] evaluate_0.14               lattice_0.20-41             purrr_0.3.4                
 [85] htmlwidgets_1.5.2           cowplot_1.1.1               bit_4.0.3                  
 [88] tidyselect_1.1.0            GSEABase_1.50.1             plyr_1.8.6                 
 [91] magrittr_1.5                R6_2.5.0                    IRanges_2.22.2             
 [94] generics_0.1.0              DelayedArray_0.14.1         DBI_1.1.0                  
 [97] pillar_1.4.6                withr_2.3.0                 RCurl_1.98-1.2             
[100] crayon_1.3.4                rmarkdown_2.5               viridis_0.5.1              
[103] progress_1.2.2              grid_4.0.2                  blob_1.2.1                 
[106] digest_0.6.25               xtable_1.8-4                tidyr_1.1.2                
[109] httpuv_1.5.4                gridGraphics_0.5-0          stats4_4.0.2               
[112] munsell_0.5.0               viridisLite_0.3.0           ggplotify_0.0.5

Reference

  1. GSVA: gene set variation analysis for microarray and RNA-Seq data

参考文章如引起任何侵权问题，可以与我联系，谢谢。