TCGA 数据分析实战 —— GSVA、ssGSEA 和单基因富集分析

TCGA 数据分析实战 —— GSVA、ssGSEA 和单基因富集分析

前言

前面，我们介绍过了差异基因的功能富集分析，今天，我们对这部分的内容作一些补充

主要介绍一下 GEVA、ssGSEA 和单基因的富集分析

GSVA

我们知道，GSEA 富集分析方法是针对两组样本来进行评估的，也就是说对基因列表的排列方式是根据基因与表型的相关度（例如，FC 值）来计算的，无法对单个样本使用

其富集分数（Enrichment Score，ES）的计算方式为

依次判断基因列表中的基因是否在基因集合中，如果在基因集合中，则 ES 加上该基因与表型的相关度，如果不是集合中的基因，则减去对应的值，最后可以计算出一个最大 ES。

Gene Set Variation Analysis（GSVA）与 GSEA 的原理类似，只是计算每个基因集合在每个样本中的 enrichment statistic（ES，或 GSVA score），其算法流程如下

不同于 GSEA 之处在于，对于不同的数据类型（只支持 log 表达值或原始的 read counts 值），假设了不同的累积密度函数（cumulative density function，CDF）

1. 芯片数据：正态分布密度函数
2. RNA-seq 数据：泊松分布密度函数

而且，GSVA 是为每个样本的每个基因计算对应的 CDF 值，然后根据该值对基因进行排序，这样，每个样本都有一个从大到小排序的基因列表

对于某一基因集合，计算其在每个样本中的 ES 值，也就是评估基因集合在基因列表中的富集情况。

例如，我们有一个排序后的样本

基因集合包含：B、E、H，我们可以绘制这样一张 K-S 分布图

x 轴为排序后的基因顺序，依照这一顺序，如果基因在集合内，则累积和会加上该基因的值（与基因的顺序有关，排名越靠前值越大），否则累积和不变。

将基因列表分为基因集内核基因集外两个集合，就可以绘制两个分布（红色和绿色曲线），分别计算两个分布之间的最大间距，以基因集内的分布值更大视为正间距（即红色曲线更高），两个最大间距之和即为该基因集的 ES 值

这样，就把基因水平的表达矩阵转换成了基因集水平的评分矩阵，可以使用差异表达基因识别算法，寻找显著差异的基因集，从而达到类似功能富集的作用

GSVA 分析

先获取基因表达矩阵，我们使用 TCGA 肺腺癌和肺鳞癌各 10 个样本的 read counts 数据

library(TCGAbiolinks)

get_prep <- function(cancer) {
   
  query <- GDCquery(
    project = cancer,
    data.category = "Transcriptome Profiling", 
    data.type = "Gene Expression Quantification", 
    workflow.type = "STAR - Counts",
    sample.type = c("Primary Tumor"),
  )
  # 选择 20 个样本
  query$results[[1]] <-  query$results[[1]][1:20,]
  GDCdownload(query)
  # 获取 read count
  exp.count <- GDCprepare(
    query,
    summarizedExperiment = TRUE,
  )
  return(exp.count)
}

luad.count <- get_prep("TCGA-LUAD")
lusc.count <- get_prep("TCGA-LUSC")

dataPrep_luad <- TCGAanalyze_Preprocessing(
  object = luad.count,
  cor.cut = 0.6,
  datatype = "unstranded"
)

dataPrep_lusc <- TCGAanalyze_Preprocessing(
  object = lusc.count,
  cor.cut = 0.6,
  datatype = "unstranded"
)
# 合并数据并使用 gcContent 方法进行标准化
dataNorm <- TCGAanalyze_Normalization(
  tabDF = cbind(dataPrep_luad, dataPrep_lusc),
  geneInfo = TCGAbiolinks::geneInfoHT,
  method = "gcContent"
)
# 分位数过滤
dataFilt <- TCGAanalyze_Filtering(
  tabDF = dataNorm,