前言

这是香港中文大学(深圳)数据科学学院于天维教授的一篇文章《A new dynamic correlation algorithm reveals novel functional aspects in single cell and bulk RNA-seq data》

作者提出了一种新的相关性的计算方式，叫做 liquid correlation coefficient，一种流动（动态的）的相关系数，计算这种相关系数的意义是去发现驱动基因间动态表达的潜变量

方法

概念1：
定义LAC（liquid correlation coefficient）的算式如下：

对于基因表达矩阵来说，LAC的计算方式为：
第一种：

第二种：
gi，gj为两个基因向量（i，j 为某个基因），r() 为皮尔斯相关性

概念2：
寻找动态的潜变量
面对一个基因表达矩阵G（中性化和标准化的矩阵） ，首先计算基因间的LAC值，然后根据LAC值筛选 gene pair；
然后计算每两个基因向量之间点乘，并且构成一个矩阵 H

作者构造 H'H矩阵，然后进行特征分解，那么分解出来的特征向量作为DC1，DC2，.....，DCn。那么这些特征向量就作为潜变量

代码

其中x$dat为表达矩阵：

# 初始化变量
array=x$dat
top.pairs.prop=0.95
max.pairs=1e6
n.fac=2
sumabsv=sqrt(max.pairs)/10
normalization="standardize"
method="PCA"

# 标准化
if(normalization == "normal score") array<-normscore.row(array)
if(normalization == "standardize") array<-normrow(array)

# 挑选 gene pair 的对数
n.pairs<-min(choose(nrow(array),2)*(1-top.pairs.prop), max.pairs)

# 计算绝对值的相关性
ccc.abs<-cor(abs(t(array)))
# 计算相关性的绝对值
abs.ccc<-abs(cor(t(array)))
# 作差，计算LAC
ccc.diff<-ccc.abs-abs.ccc
# 通过LAC再次过滤，筛选出合适的gene pair
pos<-which(ccc.diff>quantile(ccc.diff, 1-2*n.pairs/nrow(array)/nrow(array)), arr.ind=T)
pos<-pos[which(pos[,1] > pos[,2]),]

# 利用筛选的gene pair来构建 H 矩阵，计算每两个基因向量之间点乘
b<-array[pos[,1],] * array[pos[,2],]
    
# 数据转换
for(i in 1:nrow(b)) b[i,]<-normal_trafo(b[i,])

# 构建 H'H 矩阵
ccc<-t(b) %*% b
# 特征分解寻找DC
e<-eigen(ccc)
fac<-e$vec[,1:n.fac]