zhuang_gj
2020/10/11
基本原理
1.single sample GSEA 是通过扩展GESA扩展实现的，ssGSEA允许定义一个富集分数，该分数表示给定数据集内每个样本中基因集的绝对富集程度。
2.对给定样本的基因表达值进行排序归一化，并使用签名中的基因和其余基因的经验累积分布函数(ECDF)生成富集分数。
3.ssGSEA过程类似于GSEA，但是在ssGSEA列表是根据absolute expression进行排序的。
4.浓缩分数(ES)是通过对ECDFs之间的差值进行积分得到的。对于单个样本的N个基因以及给定大小为signature 
(假设是条通路的基因)数据集的，根据它们的absolute expression .L={r1, r2，…，rN}，用它们的秩替换它们。然后将列表从最高的N排列到最低的1。富集分数ES(G,S)是由一个加权的
基因的ECDF和其余
基因的ECDF之间的差值之和(积分)得到的:
ES.png
- 如果基因 属于
- 如果基因 不属于
-
- 计算每个基因的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
