【原】香浓熵值判断你的单细胞亚群是否有样品特异性

单个单细胞样品的时代早就结束了，哪怕是稀有物种珍惜样品，也很难说就一个单细胞转录组表达量的降维聚类分群结果的描述就可以发表。

不过现在有一个取巧的手段， 就是虽然是单个单细胞样品，但是里面可以拆分出来不同的来源，有点类似于混样策略。比如2021年1月发表在cancer research杂志 ： 《Single-Cell Transcriptomic Heterogeneity in Invasive Ductal and Lobular Breast Cancer Cells》， 数据链接是：https://www.ncbi.nlm./geo/query/acc.cgi?acc=GSE144320，就只有一个10X样品样本：

GSM4285803 All_CellLines

虽然这个单细胞文章仅仅是单个10X样品，但是测了8个不同的癌症细胞系，所以可以根据表达量相似性把细胞来源区分开。还有其它类似的情况是区分性别，以及区分恶性肿瘤上皮细胞和普通正常二倍体细胞。

这些算法层面的区分，就面临准确性问题。其实更好的混样，应该是每个样品样独立的标签，然后混合起来作为一个样品去做单细胞，这样就省经费了。

当然了，绝大部分科学团队其实并不会缺单细胞测序方面的资金，所以直接上大样品即可，不过多个单细胞样品往往是需要整合啦。我们也给出来了最佳实践，就是harmony流程啦 ：

sce <- NormalizeData(sce, 
 normalization.method = "LogNormalize",
 scale.factor = 1e4) 
sce <- FindVariableFeatures(sce)
sce <- ScaleData(sce)
sce <- RunPCA(sce, features = VariableFeatures(object = sce))

library(harmony)
seuratObj <- RunHarmony(sce, "orig.ident")
names(seuratObj@reductions)
seuratObj <- RunUMAP(seuratObj,  dims = 1:15, 
   reduction = "harmony")
DimPlot(seuratObj,reduction = "umap",label=T ) 

sce=seuratObj
sce <- FindNeighbors(sce, reduction = "harmony",
   dims = 1:15) 

多个样品整合起来就是为了分群而已，所以恶性的肿瘤上皮细胞往往是不需要整合，因为它具有病人特异性，也就是说每个恶性单细胞亚群都属于具体的某个病人。

但是正常单细胞亚群通常是往往是会跨越病人存在，多个病人都有同一个正常单细胞亚群。这样的跨越我很久以前是肉眼看的， 参考：CNS图表复现09—上皮细胞可以区分为恶性与否。可以很明显的看到，第1,2,7,14,21,23,25 是跨越病人的聚类情况，如下所示：

其实有一个量化指标，就是香浓熵值（Shannon entropy），我也是最近才看到，比如肺癌的单细胞文章：《Regenerative lineages and immune-mediated pruning in lung cancer metastasis》，就提到了这个Shannon entropy

• High entropy indicates that a cell state is highly reproducible across patient samples;
• Low entropy indicates that the cell state is mostly derived from one sample and is patient-specific.

也就是说 Shannon entropy 是介于 0到1的值，其中越低的Shannon entropy代表这个单细胞亚群越具有样品特异性。公式是：

Shannon entropy的公式

有了这个Shannon entropy的量化指标，如下所示，可以看到绿色的这个单细胞亚群就是具有病人特异性。

绿色的这个单细胞亚群就是具有病人特异性

有没有发现，其实它就是展现了不同单细胞亚群在不同样品的比例而已。

#　信息熵的4个量化指标的R代码实现

熵（entropy）在统计学中是一个很重要的概念，代表着信息的多少。经济学里面衡量贫富差距的基尼系数，以及环境生物学领域衡量物种多样性的辛普森多样性指数，以及免疫组库领域的D50都有异曲同工之妙。

基尼系数距离普通人的生活最近，通俗一点来理解：

• 比如有10个人，他们的月薪都是2万，那么这10个人组成的小团体的基尼系数就是0，说明没有贫富差距

• 如果他们的月薪都是3万，基尼系数也仍然是0 ，因为大家都一样。

• 但是如果他们的收入是1到10，每个人不一样，那么这个小团体的基尼系数就是0.3

• 如果前面的9个人是1到9万的收入，但是最后一个月不是10万，而是100万，这个小团体的基尼系数就是0.67

• 如果最后那个人是1个亿，基尼系数就会接近于1。

使用R代码，模拟这样的10个人小团体：

n=3
a=rep(n,10)
b1=a/sum(a)
b1 # 首先每个人的收入都是3万
plot(cumsum(b1),type = 'l')
a=1:10
a=sort(a)
b2=a/sum(a)
b2 ## 然后每个人的收入都不一样，相差一万
points(cumsum(b2),type = 'l')
a=c(1:9,100)
a=sort(a)
b3=a/sum(a)
b3 # 最后，假定其中一个人收入是100万，遥遥领先剩余的9个人
points(cumsum(b3),type = 'l')

y1=as.numeric(table(b1)/length(b1))
y2=as.numeric(table(b2)/length(b2))
y3=as.numeric(table(b3)/length(b3))

如下所示：

10个人小团体

在上图中，10个人，按照收入排序（升序）后，收入累积的占比。

那么，亲爱的读者，你可以猜测一下，我们中国的代表贫富差距的基尼系数是多少？

国际惯例把0.2以下视为收入绝对平均；0.2-0.3视为收入比较平均；0.3-0.4视为收入相对合理；0.4-0.5视为收入差距较大；当基尼系数达到0.5以上时，则表示收入悬殊。

香农信息熵

同样的10个人，同样的月薪都是2万，信息熵就是0，同样的，每个人的收入如果是3万，也不会影响信息熵就是0这个结论。但是如果10个人的收入是1到10万这10种情况，这10个人的信息熵就很大了，是3.32，但是这10个人的收入多少并不影响信息熵的结果，无论是否有一个人收入高达百万或者过亿，这个信息熵都是3.32，代表着这10个人的小团体很不一样。

所以信息熵并不能用来衡量贫富差距哦。有意思的是，如果10个人变成了100个人，同样的收入都不一样，这个时候的信息熵是6.64，也就是说信息熵居然是跟人数有关哦。但是有一个矫正后的香农信息熵，可以抹去人数的影响，代码如下：

R代码函数如下：

# 默认x 是一个群体的,每个人的收入,数值组成的向量
shannon.entropy <-function(x,type='raw'){
 if(type=='raw'){
   myfreqs <- table(x)/length(x) 
   myvec <- as.data.frame(myfreqs)[,2]
 }else{
   myvec=x
 }
  -sum(myvec * log2(myvec))
}
metric.entropy <-function(x,type='raw'){
  if(type=='raw'){
 myfreqs <- table(x)/length(x) 
 myvec <- as.data.frame(myfreqs)[,2]
  }else{
 myvec=x
  }
  -sum(myvec * log(myvec,length(x)))
}
## modify shannon.entropy to metric entropy
 

其中shannon.entropy函数其范围0<=I(x)<=1,符合直观感受。

> shannon.entropy(b1)
[1] 0
> shannon.entropy(b2)
[1] 3.321928
> shannon.entropy(b3)
[1] 3.321928
> 
> metric.entropy(b1)
[1] 0
> metric.entropy(b2)
[1] 1
> metric.entropy(b3)
[1] 1

也就是说，如果shannon.entropy接近于0，代表信息量很少，接近于1代表信息量大。

辛普森指数

辛普森多样性指数(Simpson index)，描述从一个群落种连续两次抽样所得到的个体数属于同一种的概率。

其公式如下：D=1-∑(Ni(Ni-1))/(N(N-1))，其中Ni为群落中第i种的个体数，N为群落中所有种的个体数。

如群落A，有甲99个，乙1个；群落B，有甲50个，乙50个；易得前者辛普森多样性指数=0.0198，后者希普森多样性指数=0.5000。也就是说辛普森指数越大，物种多样性越丰富，但辛普森指数最大不超过1。

R代码函数如下：

# 默认x 是一个群体的,每个人的收入,数值组成的向量
Simpson.index <-function(x,type='raw'){
  if(type=='raw'){
 myfreqs <- table(x)/length(x) 
 myvec <- as.data.frame(myfreqs)[,2]
  }else{
 myvec=x
  }
  -sum(myvec * log(myvec,length(x)))
  1-sum( myvec ^2)
}

结果如下：

> Simpson.index(b1)
[1] 0
> Simpson.index(b2)
[1] 0.9
> Simpson.index(b3)
[1] 0.9
> Simpson.index(1:100)
[1] 0.99
> Simpson.index(1:1000)
[1] 0.999

值得注意的是辛普森多样性指数和香农信息熵都不关心具体每个人的收入多少，只统计各种收入数值情况在人群出现的频率进行各自公式计算即可。

GINI系数

基尼系数是本来是一个国际通用的经济学概念，用来衡量贫富差距。基尼系数介于0-1之间，基尼系数越大，表示不平等程度越高。

• 基尼系数最大为1，表示居民之间的收入分配绝对不平均，即100%的收入被一个单位的人全部占有了；
• 基尼系数最小为0，表示居民之间的收入分配绝对平均，即人与人之间收入完全平等，没有任容何差异。

假定一定数量的人口按收入由低到高顺序排队，分为人数相等的n组，从第1组到第i组人口累计收入占全部人口总收入的比重

GINI系数简化公式

该公式是利用定积分的定义将对洛伦茨曲线的积分(面积B)分成n个等高梯形的面积之和得到的。（看不懂没有关系哈）

R代码函数实现如下：

gini.index  <-function(x){ 
  x <- sort(x)  
  G <- sum(x * 1L:length(x)) 
  G <- 2 * G/sum(x) - (length(x) + 1L)
  G/length(x)
}
> gini.index(b1)
[1] 0
> gini.index(b2)
[1] 0.3
> gini.index(b3)
[1] 0.6724138

值得注意的是，基尼系数越下，代表收入越平均，可以理解为多样性越好！

而且基尼系数关心具体每个人的收入情况，换一种说法就是基尼系数与辛普森多样性指数和香农信息熵的输入数据形式其实是不一样的：

• 输入1和2这两个数，来计算香农信息熵结果是1，辛普森多样性指数是0.5；
• 但是对基尼系数来说，输入1和2这两个数，实际上相当于输入了1个a和2个b，就是3个元素。

韩健首创的免疫组库多样性D50

这个D50其实就是饱和度曲线里面的达到50%饱和度，据韩健说是他第一次应用到免疫组库领域。把所有CDR3序列（元素）按照占比比例排序，从最高往最低累加，达到50%的总序列时候的CDR3序列种类比例，就是类似于饱和度曲线。

• D50最大为0.5，意味着全部的CDR3序列占比一致，多样性好；
• D50最小为0，意味着有且只有一种CDR3序列，多样性差。

R代码函数如下：

d50.index <-function(x,type='raw'){
  if(type=='raw'){
 myfreqs <- table(x)/length(x) 
 myvec <- sort(as.numeric(myfreqs),decreasing = T)
  }else{
 myvec=sort(x,decreasing = T)
  }
  len=length(myvec)
  state=cumsum(myvec)>sum(myvec)/2  
  (len -sum(state))/len # 
}

结果如下：

> d50.index(b1)
[1] 0
> d50.index(b2)
[1] 0.5
> d50.index(b3)
[1] 0.5

> d50.index(1:100)
[1] 0.5
> d50.index(1:1000)
[1] 0.5
> d50.index(c(1,2,2,2,3,4))
[1] 0.25

值得注意的是免疫组库多样性D50上限是0.5，因为这个时候人群收入的递减排序，不可能一半的高收入群体居然还达不到全部人群的一半的收入。通常情况是，1%的人就占社会收入的一半了，所以D50通常是0.01甚至更小值。

总结

上面我写的4个公式里面只有基尼系数计算必须输入的是数值，或者把非数值变量取频数后再进行计算。而且仅仅是只有基尼系数是越大，贫富差距越大，多样性越差。其它的数值都是越小多样性越差。

其实，聪明的你，这个时候应该是可以做出来一个总结表格。

写在文末

我在《生信技能树》，《生信菜鸟团》，《单细胞天地》的大量推文教程里面共享的代码都是复制粘贴即可使用的， 有任何疑问欢迎留言讨论，也可以发邮件给我，详细描述你遇到的困难的前因后果给我，我的邮箱地址是 [e-mail]

如果你确实觉得我的教程对你的科研课题有帮助，让你茅塞顿开，或者说你的课题大量使用我的技能，烦请日后在发表自己的成果的时候，加上一个简短的致谢，如下所示：

We thank Dr.Jianming Zeng(University of Macau), and all the members of his bioinformatics team, biotrainee, for generously sharing their experience and codes.

十年后我环游世界各地的高校以及科研院所（当然包括中国大陆）的时候，如果有这样的情谊，我会优先见你。