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

前言

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

主要介绍一下 GEVAssGSEA 和单基因的富集分析

GSVA

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

其富集分数(Enrichment ScoreES)的计算方式为

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

Gene Set Variation AnalysisGSVA)与 GSEA 的原理类似,只是计算每个基因集合在每个样本中的 enrichment statisticES,或 GSVA score),其算法流程如下

不同于 GSEA 之处在于,对于不同的数据类型(只支持 log 表达值或原始的 read counts 值),假设了不同的累积密度函数(cumulative density functionCDF

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

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

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

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


基因集合包含:BEH,我们可以绘制这样一张 K-S 分布图

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

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

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

1. GSVA 分析

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

library(TCGAbiolinks)

# 获取表达矩阵
get_count <- function(cancer, n = 10) {
  query <- GDCquery(
    project = cancer,
    data.category = "Gene expression",
    data.type = "Gene expression quantification",
    platform = "Illumina HiSeq",
    file.type  = "results",
    sample.type = c("Primary Tumor"),
    legacy = TRUE
  )
  # 选择 n 个样本
  query$results[[1]] <-  query$results[[1]][1:n,]
  GDCdownload(query)
  # 获取 read count
  exp.count <- GDCprepare(
    query,
    summarizedExperiment = TRUE,
  )
  return(exp.count)
}

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

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

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

# 将数据拆分
luad.exp <- subset(dataFilt, select = luad.count$barcode)
lusc.exp <- subset(dataFilt, select = lusc.count$barcode)

我们使用 GSVA 包提供的 gsva 函数来将基因表达矩阵转换为基因集分数矩阵

library(GSVA)
library(GSEABase)

# 读取从 GSEA 官网下载的通路数据
c2gmt <- getGmt("~/Downloads/data/pathway/c2.cp.v7.2.symbols.gmt")
# 删选出常用的这三个数据库中的通路
gene.set <- c2gmt[grep("^KEGG|REACTOME|BIOCARTA", names(c2gmt)),]
# gsva 分析,read counts 使用泊松分布,通路至少包含 10 个基因
gs.exp <- gsva(dataFilt, gene.set, kcdf = "Poisson", min.sz = 10)

虽然 GSVAdata 包提供了通路数据 c2BroadSets 是基因 ID,但我们的基因表达数据的行是基因 Symbol,所以通路信息也必须是 Symbol 格式,要进行格式转换,比较麻烦

所以我们使用 GSEABase 包提供的 getGmt 函数来读取从 GSEA 官网下载的 C2 通路信息

得到结果如下,共包含 1511 条通路

然后,使用差异基因识别方法

2. 差异分析

我们使用 limma 分析差异通路

DEA.gs <- TCGAanalyze_DEA(
  mat1 = gs.exp[, colnames(luad.exp)],
  mat2 = gs.exp[, colnames(lusc.exp)],
  metadata = FALSE,
  pipeline = "limma",
  Cond1type = "LUAD",
  Cond2type = "LUSC",
  fdr.cut = 0.05,
  logFC.cut = 0.5,
)

通过设置 FDR = 0.05,logFC = 0.5 共筛选出 40 条差异通路

查看通路的火山图


我们可以一起看下基因的火山图

DEA.gene <- TCGAanalyze_DEA(
  mat1 = luad.exp,
  mat2 = lusc.exp,
  metadata = FALSE,
  pipeline = "limma",
  Cond1type = "LUAD",
  Cond2type = "LUSC",
  fdr.cut = 0.01,
  logFC.cut = 1
)

总共识别出 804 个差异表达基因

ssGSEA

single sample Gene Set Enrichment Analysis (ssGSEA) 是针对单个样本进行 GSEA 分析,其基因列表的排序方式和 ES 的计算方式都是依赖于样本中基因的表达值,而不再是依赖基因与表型的相关度

使用方式也很简单,只要在 gsva 函数中指定 method = "ssgsea",例如

res.ssgsea <- gsva(dataFilt, gene.set, method = "ssgsea", kcdf = "Poisson", min.sz = 10)

也可以进行差异分析

DEA.ssgsea <- TCGAanalyze_DEA(
  mat1 = res.ssgsea[, colnames(luad.exp)],
  mat2 = res.ssgsea[, colnames(lusc.exp)],
  metadata = FALSE,
  pipeline = "limma",
  Cond1type = "LUAD",
  Cond2type = "LUSC",
  fdr.cut = 0.05,
  logFC.cut = 0.1,
)

或者绘制热图

annotation_col <- data.frame(sample = rep(1:2, each = 10))
rownames(annotation_col) <- colnames(res.ssgsea)
pheatmap(
  res.ssgsea[rownames(DEA.ssgsea),],
  show_colnames = F,
  # 不展示行名
  cluster_rows = F,
  # 不对行聚类
  cluster_cols = F,
  # 不对列聚类
  annotation_col = annotation_col,
  # 加注释
  cellwidth = 5,
  cellheight = 5,
  # 设置单元格的宽度和高度
  fontsize = 5
)

单基因富集分析

单基因富集分析并不是说拿单个基因来进行富集分析,单个基因怎么能进行富集分析呢?一个基因根本没法进行统计检验。

其实,这里说的单基因并不是拿单个基因来富集,而是基于单个基因来进行富集分析,这个“基于”,就是以单个基因为基础,向外扩展,抓取与其相关的基因,然后用这些相关的基因来进行功能富集

所以,要理解这个单基因富集分析的意思,这样一说就已经很明了了。针对单个基因我们可以做什么?

主要有两种做法:

  1. 定性分组:我们可以根据给定基因的表达值对样本进行分组,然后识别在两组样本之间差异表达的基因,最后用这些差异表达基因来进行功能富集

  2. 定量相关:通过计算其他基因与目标基因表达之间的相关性,将具有显著相关的基因作为一个集合,也可以进行富集分析

1. 定性分组

我们以 CCDC134 基因为例,以该基因表达值的中位值来对样本进行分组

gene <- "CCDC134"
gene.exp <- dataFilt[gene,]

label <- if_else(gene.exp < median(gene.exp), 0, 1)

group.low <- dataFilt[,label == 0]
group.high <- dataFilt[,label == 1]

识别两组样本之间的差异表达基因

DEGs <- TCGAanalyze_DEA(
  mat1 = group.low,
  mat2 = group.high,
  metadata = FALSE,
  pipeline = "limma",
  Cond1type = "CCDC134_Low",
  Cond2type = "CCDC134_High",
  fdr.cut = 0.01,
  logFC.cut = 1,
)

共识别出 873 个差异表达基因

2. 定量相关

我们对其他基因与 CCDC134 基因进行相关性检验,由于基因较多,我们使用并行的方式来计算

library(future.apply)

batch_cor <- function(exp, gene){
  y = as.numeric(exp[gene,])
  gene_list = rownames(exp)
  gene_list = gene_list[rownames(exp) != gene]
  do.call(rbind, future_lapply(gene_list, function(x){
    ct  <- cor.test(as.numeric(exp[x,]), y, type='spearman')
    data.frame(key = gene, gene = x, cor = ct$estimate,p.value = ct$p.value )
  }))
}

plan(multiprocess)
system.time(res.cor <- batch_cor(dataFilt, gene))

对结果进行过滤,筛选出显著相关且相关系数的绝对值大于 0.6 的基因,共筛选出 232 个基因

cor.genes <- filter(res.cor, p.value < 0.05 & abs(cor) > 0.6)

3. 富集分析

格式化识别出的差异基因

library(clusterProfiler)
library(org.Hs.eg.db)
library(enrichplot)

gene.id <- bitr(
  rownames(DEGs), fromType = "SYMBOL",
  toType = "ENTREZID",
  OrgDb = org.Hs.eg.db
)

go <- enrichGO(
  gene = gene.id,
  OrgDb = org.Hs.eg.db,
  ont = "ALL",
  pAdjustMethod = "BH",
  qvalueCutoff = 0.05,
  readable = T
)
dotplot(go)

gene_info <- DEGs %>%
  rownames_to_column(var = "SYMBOL") %>%
  inner_join(., gene.id[,1:2], by = "SYMBOL") %>%
  # 必须降序
  arrange(desc(logFC))

# 构造输入数据格式
geneList <- gene_info$logFC
names(geneList) <- as.character(gene_info$ENTREZID)

go2 <- gseGO(
  geneList     = geneList,
  OrgDb        = org.Hs.eg.db,
  ont          = "ALL",
  minGSSize    = 10,
  maxGSSize    = 500,
  pvalueCutoff = 0.1,
  verbose      = FALSE
)

两种富集方法都没有富集到 go 通路,对于相关基因也是没有富集到通路的。选的这个基因不行

©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 206,602评论 6 481
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 88,442评论 2 382
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 152,878评论 0 344
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 55,306评论 1 279
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 64,330评论 5 373
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 49,071评论 1 285
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 38,382评论 3 400
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 37,006评论 0 259
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 43,512评论 1 300
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 35,965评论 2 325
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 38,094评论 1 333
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 33,732评论 4 323
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 39,283评论 3 307
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 30,286评论 0 19
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 31,512评论 1 262
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 45,536评论 2 354
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 42,828评论 2 345

推荐阅读更多精彩内容