ssGSEA算法的原理和TCGA数据的应用

最编程 2024-03-04 17:25:46

...

zhuang_gj

2020/10/11

ssGSEA算法原理及应用
- 基本原理
- 实例演示
  - 1.数据清洗和整理
  - 1.2 FPKM 转为TPM
  - 2.获得基因集
  - 3.运行ssGSEA得到得到单个样本在不同基因集中的ES

基本原理

1.single sample GSEA 是通过扩展GESA扩展实现的，ssGSEA允许定义一个富集分数，该分数表示给定数据集内每个样本中基因集的绝对富集程度。

2.对给定样本的基因表达值进行排序归一化，并使用签名中的基因和其余基因的经验累积分布函数(ECDF)生成富集分数。

3.ssGSEA过程类似于GSEA，但是在ssGSEA列表是根据absolute expression进行排序的。

4.浓缩分数(ES)是通过对ECDFs之间的差值进行积分得到的。对于单个样本的N个基因以及给定大小为signature $N_G$ (假设是条通路的基因)数据集的，根据它们的absolute expression .L={r1, r2，…，rN}，用它们的秩替换它们。然后将列表从最高的N排列到最低的1。富集分数ES(G,S)是由一个加权的 $P_G^w$ 基因的ECDF和其余 $P_{NG}$ 基因的ECDF之间的差值之和(积分)得到的:

ES.png

如果基因属于 $N_G$ ，则 $P_G^w$ 是该基因的表达量/ $N_G$ 表达量的总和

如果基因 不属于 $N_G$ ，则 $P_{NG}$ 就等于1/ $N$ - $N_G$

计算每个基因的ES，然后根据最大差异来得到这条通路的ES(绝对值最大的ES作为通路的ES.)。

指数(α)设置为¼,并添加一个适中的重量等级。在常规的GSEA中使用了类似的浓缩分数，但是权重通常设置为1。

实例演示

分析流程大致是：
1.清洗数据，得到clean的matrix（行为基因列为样本）
2.整理参考基因集
3.运行ssGSEA得到得到单个样本在不同基因集中的ES

1.数据清洗和整理

rm(list=ls())
library(tidyverse)
fpkm <- data.table::fread(file = "./TCGA-KICH.htseq_fpkm.tsv")%>%column_to_rownames("Ensembl_ID")
## 去除低表达的探针
dat<- fpkm[!apply(fpkm,1,sum)==0,]
## 过滤重复的probe
boxplot(dat[,2:10],las=2)

image.png

dat$Median=apply(dat,1,median)
dat <- dat %>% 
  rownames_to_column("gene_id")#去除表达矩阵ensembolID   "."后的数字separate(Ensembl_ID,into = c("gene_id"),sep = "\\.")

load(file = "C:/Diskdata/Bioinfomatics_analyse/gtf_df.Rdata")
## 注释差异基因
exprSet <- gtf_df %>% 
  ## 筛选gene
  dplyr::filter(type=="gene") %>%
  ## 筛选基因名称和类型
  dplyr::select(c(gene_name,gene_id,gene_type)) %>% 
  ## 合并
  dplyr::inner_join(dat,by ="gene_id")

# # 筛选mRNA
# exprSet <- alldat %>% 
#   dplyr::filter(gene_type == "protein_coding")
# 去除重复探针
exprSet <- exprSet[order(exprSet$gene_name,exprSet$Median,decreasing = T),]# 降序
exprSet <- exprSet [!duplicated(exprSet$gene_name),];
rownames(exprSet) <- exprSet[,"gene_name"]
exprSet <- exprSet[,4:(ncol(exprSet)-1)]
dim(exprSet)

## [1] 52167    89

head(exprSet[,1:4])

##          TCGA-KL-8332-01A TCGA-KN-8426-11A TCGA-KN-8423-01A TCGA-KL-8327-01A
## ZZZ3           1.80859149        2.6352506       2.19545866         1.589469
## ZZEF1          2.10315141        2.5283014       2.28908861         1.322511
## ZYX            4.00910991        5.8957148       3.57184237         3.411580
## ZYG11B         3.61346486        2.7084147       3.96803884         2.169499
## ZYG11AP1       0.00000000        0.0000000       0.00000000         0.000000
## ZYG11A         0.02453941        0.3772259       0.01565521         0.000000

1.2 FPKM 转为TPM

# USUC 下载的TCGA 数据是log2(fpkm+1)
boxplot(exprSet[,2:10],las=2)

image.png

dat <- 2^exprSet-1
boxplot(dat[,2:10],las=2)

image.png

fpkmToTpm <- function(fpkm)
{
  exp(log(fpkm) - log(sum(fpkm)) + log(1e6))
}
# <https://haroldpimentel.wordpress.com/2014/05/08/what-the-fpkm-a-review-rna-seq-expression-units/>
tpms <- apply(dat,2,fpkmToTpm)
tpms[1:3,1:4]

##       TCGA-KL-8332-01A TCGA-KN-8426-11A TCGA-KN-8423-01A TCGA-KL-8327-01A
## ZZZ3          5.369558         18.27102         9.942082         2.553463
## ZZEF1         7.071742         16.71511        10.794907         1.907436
## ZYX          32.396140        205.17275        30.243626        12.251646

boxplot(tpms[,1:10],las=2)

image.png

colSums(tpms)

## TCGA-KL-8332-01A TCGA-KN-8426-11A TCGA-KN-8423-01A TCGA-KL-8327-01A 
##            1e+06            1e+06            1e+06            1e+06 
## TCGA-KO-8409-01A TCGA-KN-8430-01A TCGA-KL-8329-01A TCGA-KL-8338-01A 
##            1e+06            1e+06            1e+06            1e+06 
## TCGA-KN-8432-11A TCGA-KN-8429-01A TCGA-KO-8415-11A TCGA-KN-8436-11A 
##            1e+06            1e+06            1e+06            1e+06 
## TCGA-KN-8432-01A TCGA-KN-8424-11A TCGA-KL-8324-11A TCGA-KN-8437-01A

# 对TPM数据进行log转换
logTPM <- log2(tpms+1)
boxplot(logTPM[,1:10],las=2)

image.png

head(logTPM[,1:3])

##          TCGA-KL-8332-01A TCGA-KN-8426-11A TCGA-KN-8423-01A
## ZZZ3           2.67119317         4.268361       3.45181542
## ZZEF1          3.01288006         4.146909       3.56009213
## ZYX            5.06160945         7.687710       4.96548998
## ZYG11B         4.65026873         4.350769       5.38125597
## ZYG11AP1       0.00000000         0.000000       0.00000000
## ZYG11A         0.05213974         1.033824       0.04305977

save(logTPM,file = "./step1_exprSet logTPM.Rdata")

2.获得基因集

# 下载28种免疫细胞的参考基因集 <http://cis.hku.hk/TISIDB/data/download/CellReports.txt>
rm(list=ls())
library(tidyverse)
options(stringsAsFactors = F)
geneSet <- read.csv("CellMarker.txt",header = F,sep = "\t",) # 用EXCEL打开删除NA列
class(geneSet)

## [1] "data.frame"

geneSet <- geneSet %>%
  column_to_rownames("V1")%>%t()
a <- geneSet
a <- a[1:nrow(a),]
set <- colnames(a)
l <- list()
#i <- "Activated CD8 T cell"
for (i in set) {
  x <-  as.character(a[,i])
  x <- x[nchar(x)!=0]
  x <-  as.character(x)
  l[[i]] <-x
}
save(l,file = "./gene_set.Rdata")

3.运行ssGSEA得到得到单个样本在不同基因集中的ES

# 读取基因的表达矩阵
rm(list=ls())
library(GSVA)
library(limma)
## 载入数据
load(file = "./step1_exprSet logTPM.Rdata")
load(file = "./gene_set.Rdata")
dat <- as.matrix(logTPM)#基因的表达量需要是矩阵，行为基因，列为样本
#开始进行ssGSEA
ssgsea<- gsva(dat, l,method='ssgsea',kcdf='Gaussian',abs.ranking=TRUE)

## Estimating ssGSEA scores for 28 gene sets.
## 
  |====================================================================  |  98%
  |                                                                            
  |===================================================================== |  99%
  |                                                                            
  |======================================================================| 100%

#基因集需要是list为对象。默认情况下，kcdf="Gaussian"，适用于输入表达式值连续的情况，如对数尺度的微阵列荧光单元、RNA-seq log-CPMs、log-RPKMs或log-TPMs。当输入表达式值是整数计数时，比如那些从RNA-seq实验中得到的值，那么这个参数应该设置为kcdf="Poisson"

# Min-Max标准化是指对原始数据进行线性变换，将值映射到[0，1]之间
# 这里是将每个样本中不同的免疫细胞比例标准化到0-1之间
ssgsea.1 <- ssgsea
for (i in colnames(ssgsea)) {
  #i <- colnames(ssgsea)[1]
  ssgsea.1[,i] <- (ssgsea[,i] -min(ssgsea[,i]))/(max(ssgsea[,i] )-min(ssgsea[,i] ))

}
apply(ssgsea.1[,1:6], 2, range)

##      TCGA-KL-8332-01A TCGA-KN-8426-11A TCGA-KN-8423-01A TCGA-KL-8327-01A
## [1,]                0                0                0                0
## [2,]                1                1                1                1
##      TCGA-KO-8409-01A TCGA-KN-8430-01A
## [1,]                0                0
## [2,]                1                1

library(pheatmap)
pheatmap(ssgsea.1,
         show_colnames = F,
         cluster_rows = T,
         cluster_cols = T,
         fontsize=5)#画热图

image.png

若上述代码分析有误，欢迎纠错和补充。

感谢jimmy老师及其生信技能树团队分析的学习笔记

https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2783335/

https://blog.****.net/weixin_42792088/article/details/103971069

上一篇：单细胞亚群分析揭示了炎症性癌症相关成纤维细胞（iCAFs）在肿瘤进展中的作用

下一篇：朱莉娅复数和有理数 |-示例

ssGSEA算法的原理和TCGA数据的应用

zhuang_gj

2020/10/11

基本原理

1.single sample GSEA 是通过扩展GESA扩展实现的，ssGSEA允许定义一个富集分数，该分数表示给定数据集内每个样本中基因集的绝对富集程度。

2.对给定样本的基因表达值进行排序归一化，并使用签名中的基因和其余基因的经验累积分布函数(ECDF)生成富集分数。

3.ssGSEA过程类似于GSEA，但是在ssGSEA列表是根据absolute expression进行排序的。

实例演示

1.数据清洗和整理

1.2 FPKM 转为TPM

2.获得基因集

3.运行ssGSEA得到得到单个样本在不同基因集中的ES

若上述代码分析有误，欢迎纠错和补充。

[数据结构] 算法的效率（时间复杂性和空间复杂性） - V．oj 习题的复杂性

贪婪算法在 Python、JavaScript、Java、C++ 和 C# 中的多种实现及其在硬币变化、分数骑士、活动选择和使用哈夫曼编码的最小生成树问题中的应用实例

PACS 系统中的三维重建技术：原理、实施和应用

人工智能（pytorch）构建模型 28 - 基于变压器的端到端目标检测的实际应用 DETR 模型、DETR 的原理和结构

IP 子网划分的原理和应用

SSH 加密原理、RSA 非对称加密算法的学习和理解

曲奇从初级到高级：透彻理解文章的原理和应用

负载平衡的原理和算法 - I. 定义

虚拟网络设备的网络隔离机制：深入分析原理、影响和应用场景