github-actions[bot]
commented
4 years ago ssgsea算法在量化免疫浸润时的运用以及原理 by 果子学生信
ssgsea在免疫浸润的使用
在量化免疫浸润的时候，我们有两种主要的方法
第一种是CIBERSORT，他的结果是各个免疫细胞在一个样本中的占有率。
量化免疫浸润时CIBERSORT的注意事项
第二种是ssGSEA，单样本GSEA分析。
使用起来简单，而且高度可定制。
以下是一句代码。
gsva_data <- GSVA::gsva(expr,cellMarker, method = "ssgsea")
用到的是GSVA包，expr是表达矩阵，cellMarker是基因集列表
计算的结果是个矩阵，行是基因集，列是样本。
那么接下来的问题是：
ssgsea返回的结果的值究竟是什么？
为什么用不同数量的基因集去计算同一个样本，返回的结果不一样？
计算出来的结果可以跨平台比较么？
芯片数据和RNA-seq数据，用起来有什么注意事项？
ssgsea是如何用基因集给样本打分的？
ssgsea手工复现
简单来说，就是使用随机漫步的方法，来看基因集的基因的累计密度与非基因集的基因的累计密度的偏离程度。
看起来复杂，手工计算后就会好很多。
现在设置三个参数
第一，用X代替表达矩阵expr
第二，用gene_sets代替基因集cellMarker
第三，设置一个参数用于改变数据权重alpha
X=expr
gene_sets= cellMarker
alpha = 0.25
获取基因名称，以及基因个数
row_names = rownames(X)
num_genes = nrow(X)
获取每个基因集中的基因在当前表达矩阵基因中的位置
gene_sets = lapply(gene_sets, function(genes) {which(row_names %in% genes)})
本来基因集中都是基因名称
现在基因集全部都是位置信息
此时295代表基因ADAM28，在基因中是第295位
可以验证出来，确实如此
row_names[295]
现在，我们要把四列样本，每一列的重新排序，排序的标准就是基因的表达量，从高到底
R <- apply(X, 2, function(x) as.integer(rank(x)))
此时仅仅给出每个样本中，每个基因的位置，没有物理层面改变顺序
准备工作已经做完了，现在要做的就是计算一个每个样本中每个数据集的值
我们先简化到，计算第一个基因集在第一个样本中的富集分数
首先，对第一列的样本排序，获取每个基因从大到小排序后原先的位置
实际上就是理解order的用法，他跟sort比，不仅仅把基因排序，同时返回的是排序后的基因在原来集合中的位置
gene_ranks <- order(R[,1], decreasing = TRUE)
获取第一个基因集
gene_set_idx <- gene_sets[[1]]
找出在基因集中的基因和不在基因集中的基因
indicator_pos = gene_ranks %in% gene_set_idx
indicator_neg = !indicator_pos
把第一个样本的基因从大到小排序
gene <- sort(R[,1], decreasing = TRUE)
然后给他打分，在基因集的给1分，不在的不给分，而且给分是按照位置的权重来给的。
rank_alpha = (gene* indicator_pos) ^ alpha
为什么要用alpha来计算乘方呢，GSEA的alpha是1，
这是主要是为了改变头尾的权重，作图理解一下(这是我此刻的理解水平，以后会提高)
不用alpha处理是这样的,
plot(gene* indicator_pos)
如果处理了之后是这样的
plot((gene* indicator_pos) ^ alpha)
靠感觉就是，基因集中的这25个基因的值，权重更大了。
计算基因集中的基因的累计密度
step_cdf_pos = cumsum(rank_alpha)/sum(rank_alpha)
plot(step_cdf_pos)
计算非基因集的基因的累计密度
step_cdf_neg = cumsum(indicator_neg) / sum(indicator_neg)
plot(plot(step_cdf_pos))
现在每个基因都可以计算出两种情况的差值
而这个差值的和就是最终的富集分数
如果是GSEA算法，那么最终的是是差值的最大值，这一点是不一样的。
这个方法叫作随机漫步
step_cdf_diff = step_cdf_pos - step_cdf_neg
sum(step_cdf_diff)
这个值算出来是633.9068
这个值是第一个基因集在第一个样本中的数值，
正式计算时可以批量求出每个基因集在这个样本中的富集分数
es_sample = sapply(gene_sets, function(gene_set_idx) {
# pos: match (within the gene set)
# neg: non-match (outside the gene set)
indicator_pos = gene_ranks %in% gene_set_idx
indicator_neg = !indicator_pos
rank_alpha = (R[gene_ranks,1] * indicator_pos) ^ alpha
step_cdf_pos = cumsum(rank_alpha) / sum(rank_alpha)
step_cdf_neg = cumsum(indicator_neg) / sum(indicator_neg)
step_cdf_diff = step_cdf_pos - step_cdf_neg
# Use ssGSEA or not
sum(step_cdf_diff)
})
es_sample
现在这个值很大，需要矫正，就是用要处理最大值和最小值的差值，我们计算完了再来矫正。
现在，我们接着这句命令，往下走
R <- apply(X, 2, function(x) as.integer(rank(x)))
批量计算
es = apply(R, 2, function(R_col) {
gene_ranks = order(R_col, decreasing = TRUE)
# Calc es for each gene set
es_sample = sapply(gene_sets, function(gene_set_idx) {
# pos: match (within the gene set)
# neg: non-match (outside the gene set)
indicator_pos = gene_ranks %in% gene_set_idx
indicator_neg = !indicator_pos
rank_alpha = (R_col[gene_ranks] * indicator_pos) ^ alpha
step_cdf_pos = cumsum(rank_alpha) / sum(rank_alpha)
step_cdf_neg = cumsum(indicator_neg) / sum(indicator_neg)
step_cdf_diff = step_cdf_pos - step_cdf_neg
# Use ssGSEA or not
sum(step_cdf_diff)
})
unlist(es_sample)
})
这就是结果所有样本的结果
现在进行矫正
es = es / diff(range(es))
得到的结果，跟正式运行时是一样的。
增加行名和列名
rownames(es) = names(gene_sets)
colnames(es) = colnames(X)
至此，我们用手工计算的方式，复现了ssgsea的过程。
注意事项
根据计算，得到的结果是富集分数，而且是矫正后的富集分数，
所以，可以称为NES(normalized enrichment score)
又因为矫正时用的是最大值和最小值的差值，那么只要基因集的数量不一样，
这个差值就会不一样，那么最终的结果就会不一样。
这是手工复现时发现的。
这个结果能跨平台比较么，
我觉得不能，因为他计算时强烈的依赖当前基因的测量值以及基因的个数，
不同的平台没办法满足这两个条件。
计算时，我们知道，他一上来就根据测量值把基因给排了序，所以
芯片，或者测序数据，取不取log，我觉得都对结果没有影响。
我这么说，有可能不对，因为GSVA的论文中，给出了建议：
芯片数据需要Normalization(具体方案已经在安排)
测序数据RNA-seq需要是counts或者TPM，也可以使用log化后的数据。
但是，对于ssgsea的算法而言，除了要保证一开始的数据是矫正过的(因为最终是样本间比较，所以芯片样本间的矫正是必须的)。
其代码中并没有发现对数据有什么严格的要求，甚至我觉得log都不重要。
详细的讲解视频，在更新完了xcell，小鼠免疫浸润，后连同之前的CIBERSORT一起更新在答疑群中。
到时候，我应该会有新的理解。
https://mp.weixin.qq.com/s/TnUq-NfIc7oQTjnflVqs4g