ssGSEA算法原理及应用TCGA数据

作者: 庄小尽 | 来源:发表于2020-10-11 20:18 被阅读0次

    zhuang_gj

    2020/10/11

    基本原理

    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

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

    1. 感谢jimmy老师及其生信技能树团队分析的学习笔记
    2. https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2783335/
    3. https://blog.csdn.net/weixin_42792088/article/details/103971069

    相关文章

      网友评论

        本文标题:ssGSEA算法原理及应用TCGA数据

        本文链接:https://www.haomeiwen.com/subject/zstypktx.html