前言
这是香港中文大学(深圳)数据科学学院于天维教授的一篇文章《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]