hello，大家好，之前在文章10X单细胞（10X空间转录组）数据分析之NMF（非负矩阵分解）中介绍了NMF在单细胞数据分析中的运用，简单总结一下就是NMF算法的基本思想是将原始非负矩阵分解为两个非负矩阵的乘积从而对原始高维矩阵进行降维表示。
我们在做WGCNA的时候，得到的模块未必与细胞类型或者样本关联度很高，但是我们又想得到基因模块与某一种细胞类型或者一种样本类型高度关联而与其他类型几乎没关系，怎么做呢？ cNMF就是为了这个而来。
python代码如下：
import numpy as np
import random
def nmf(X, r, k, e):
'''
X是原始矩阵
r是分解的两个非负矩阵的隐变量维度，要远小于原始矩阵的维度
k是迭代次数
e是理想误差
input X
output U,V
'''
d, n = X.shape
#print(d,n)
#U = np.mat(random.random((d, r)))
U = np.mat(np.random.rand(d, r))
#V = np.mat(random.random((n, r)))
V = np.mat(np.random.rand(n, r))
#print(U, V)
x = 1
for x in range(k):
print('---------------------------------------------------')
print('开始第', x, '轮迭代')
#error
X_pre = U * V.T
E = X - X_pre
#print E
err= 0.0
for i in range(d):
for j in range(n):
err += E[i,j] * E[i,j]
print('误差：', err)
if err < e:
break
#update U
a_u = U * (V.T) * V
b_u = X * V
for i_1 in range(d):
for j_1 in range(r):
if a_u[i_1, j_1] != 0:
U[i_1, j_1] = U[i_1, j_1] * b_u[i_1, j_1] / a_u[i_1, j_1]
#print(U)
#update V
a_v = V * (U.T) * U
b_v = X.T * U
print(r,n)
for i_2 in range(n):
for j_2 in range(r):
if a_v[i_2,j_2] != 0:
V[i_2,j_2] = V[i_2,j_2] * b_v[i_2,j_2] / a_v[i_2,j_2]
#print(V)
print('第', x, '轮迭代结束')
return U, V
if __name__ == "__main__":
X = [[5, 3, 2, 1, 2, 3],
[4, 2, 2, 1, 1, 5],
[1, 1, 2, 5, 2, 3],
[1, 2, 2, 4, 3, 2],
[2, 1, 5, 4, 1, 1],
[1, 2, 2, 5, 3, 2],
[2, 5, 3, 2, 2, 5],
[2, 1, 2, 5, 1, 1], ]
X = np.mat(X)
#print(X)
U, V = nmf(X, 2, 100, 0.001)
print(U*V.T)
其实在单细胞数据分析的过程中，细胞基因的表达模式受多种因素干扰，基因通过协同作用维持细胞类型特有的生物学特征，而相互协同的基因作为基因表达模块（GEP）一起诱导和相应内外部信号，执行复杂的细胞功能。功能基因模块可以出现在多个不同的细胞类型中，而细胞类型基因模块代表一个单一的细胞类型，因此可利用这一事实来区分细胞类型基因模块和功能基因模块。通过cNMF分析，可同时推断出细胞类型相关和功能相关的GEPs，从而改进marker基因的推断，使得细胞类型鉴定更加准确。
什么是cNMF？？
约束非负矩阵分解（CNMF）算法，该算法将标签信息作为附加的硬约束，使得具有相同类标签信息的数据在新的低维空间中仍然保持一致。
但是，CNMF算法对于无标签数据样本没有任何约束，因此在很少的标签信息时它的性能受限，并且对于类中只有一个样本有标签的情形，CNMF算法中构建的约束矩阵将退化为单位矩阵，失去其意义。
实现代码
import numpy as np
import random
def cnmf(X, C, r, k, e):
'''
X是原始矩阵，维度为d*n
C是有标签样本指示矩阵，维度为l*c（l——有标签的样本数量，c——类别数量）
r是分解的两个非负矩阵的隐变量维度，要远小于原始矩阵的维度
k是迭代次数
e是理想误差
input X,C
output U,V
'''
d, n = X.shape
l, c = C.shape
#计算A矩阵
I = np.mat(np.identity(n-l))
A = np.zeros((n, n + c - l))
for i in range(l):
for j in range(c):
A[i,j] = C[i,j]
for i2 in range(n-l):
A[l+i2, c+i2] = I[i2, i2]
A = np.mat(A)
U = np.mat(np.random.rand(d, r))
Z = np.mat(np.random.rand(n + c - l, r))
#print(A)
x = 1
for x in range(k):
print('---------------------------------------------------')
print('开始第', x, '轮迭代')
#error
X_pre = U * (A*Z).T
E = X - X_pre
#print E
err= 0.0
for i in range(d):
for j in range(n):
err += E[i,j] * E[i,j]
print('误差：', err)
if err < e:
break
#update U
a_u = U * Z.T * A.T * A * Z
b_u = X * A * Z
for i_1 in range(d):
for j_1 in range(r):
if a_u[i_1, j_1] != 0:
U[i_1, j_1] = U[i_1, j_1] * b_u[i_1, j_1] / a_u[i_1, j_1]
#print(U)
#update Z
#print(Z.shape,n,r)
a_z = A.T * A * Z * U.T * U
b_z = A.T* X.T * U
for i_2 in range(n + c - l):
for j_2 in range(r):
#print(i_2, j_2, a_z[i_2,j_2])
if a_z[i_2,j_2] != 0:
Z[i_2,j_2] = Z[i_2,j_2] * b_z[i_2,j_2] / a_z[i_2,j_2]
#print(V)
print('第', x, '轮迭代结束')
V = (A*Z).T
return U, V
if __name__ == "__main__":
X = [[5, 3, 2, 1, 2, 3],
[4, 2, 2, 1, 1, 5],
[1, 1, 2, 5, 2, 3],
[1, 2, 2, 4, 3, 2],
[2, 1, 5, 4, 1, 1],
[1, 2, 2, 5, 3, 2],
[2, 5, 3, 2, 2, 5],
[2, 1, 2, 5, 1, 1],]#8*6,6个样本
X = np.mat(X)
C = [[0, 0, 1],
[0, 1, 0],
[0, 1, 0],
[1, 0, 0],]#4*3，假设有4个样本有标签，总共有三类标签
#print(X)
C = np.mat(C)
U, V = cnmf(X, C, 2, 100, 0.01)
print(U.shape, V.shape)
print(U * V)
cNMF在单细胞数据分析中的运用(文献用法)
cNMF的目的
To better characterize sample heterogeneity in whole-tumor samples
图片.png
图注： Heatmap showing the correlation between signatures obtained for all 13 whole-tumour samples.
用法
The cNMF algorithm was applied individually to each whole-tumor sample with some modifications。In brief, nonnegative matrix factorization was run 100 times for k from 2 to 15 signatures. For each k, the 100 repetitions are clustered in k groups. We expect a stable clustering solution would produce tight clusters with one signature per cluster for each of the 100 repetitions. The proportion of repetitions with one signature per cluster was called reproducibility. Clustering of the signatures was done by k means with a constraint of uniform cluster sizes, prioritizing higher correlations. The largest k with a reproducibility above 0.9 was chosen. For a chosen k, we confirmed the clustering solution was appropriate by running tSNE on the signatures it generated. The final signatures for a given sample was obtained by averaging the signature repetitions within a cluster, excluding repetitions with poor reproducibility (the ones which did not produce a signature per cluster). We obtained between five and nine final signatures per sample, 79 signatures in total. From the inter-signature Pearson correlations, we used hierarchical clustering to find trends of signatures (Fig. 1e, hierarchical tree). Six main groups emerged. In one of these groups, important variations in gene weights were observed: signatures characterized by OLIG2 and ASCL1, for example, had less DCX and STMN1, and vice versa. We reclustered this group in two, yielding the final seven groups (Fig. 1e and Supplementary Fig. 2f,就是上图).
We identified the most characteristic genes of each signature group by ranking genes according to their relative signal to noise ratio (snr) and chose the top 40 for the heatmap.
图片.png
关于信噪比（SNR），大家可以参考我的文章单细胞数据信噪比（Signal-to-noise ratio，SNR）。
We scored each signature according to the TCGA using the method described above (Classifying cells by TCGA subtype). A given signature was labeled with the subtype yielding the highest score.（下图）
图片.png
看看文献的原始代码(参照上述cNMF的实现方式)
% cNMF_seperate.m
% use NMR to find signatures in scRNA data sample by sample. Based on cNMF (Kotliar et al.) ###看来是样本作为主要标签。
% Author: Charles Couturier
% Date: August 3, 2018
%% prepare data and initialize parameters
% genes: list of gene names and ensembl code; logm: log of raw counts; sample: vector
% identifying cells by sample; sample_id: a list of the samples
load('gbm.mat','genes','logm','sample','sample_id')
kvals = 2:15;
L = length(kvals);
nreps = 100;
allHs = cell(L,1);
allHc = cell(L,1);
allSil = cell(L,1);
allE = cell(L,1);
allR = cell(L,1);
allCL = cell(L,1);
for s = 1:length(sample_id)
fprintf('Sample %s in progress, %i of %i\n',sample_id{s},s,length(sample_id))
data = logm(sample==s,:);
% initialize
m = size(data,1);
n = size(data,2);
E = zeros(L,1);
R = zeros(L,1);
Sil = zeros(L,1);
consensus_H = cell(L,1);
allH = cell(L,1);
CL = cell(L,1);
%% NMF
opt = statset('MaxIter',1e6);
parfor kid = 1:L
k = kvals(kid);
hi = zeros(k*nreps,n);
consensus_hi = zeros(k,n);
D = zeros(nreps,1);
for i = 1:nreps
[~,hi(i*k-k+1:i*k,:),D(i)]=nnmf(data,k,'replicates',20,'options',opt,'algorithm','mult');
end
% k clusters
cl = uniform_kmeans(hi,k,'Replicates',10);
% find replicates with 1 of each cluster
goodrep = all(sort(reshape(cl,k,nreps),1) == [1:k]',1);
goodrep = repmat(goodrep,k,1);
goodrep = goodrep(:);
% consensus hi finds median of each component cluster, removing
% bad replicates (see above)
for i = 1:k
consensus_hi(i,:) = median(hi(cl==i & goodrep,:),1);
end
allH{kid} = hi;
consensus_H{kid} = consensus_hi;
Sil(kid) = mean(silhouette(hi,cl));
E(kid) = mean(D);
R(kid) = sum(goodrep)/(k*nreps);
CL{kid} = cl;
fprintf('Run for %i components completed\n',k)
end
allHs{s} = allH;
allHc{s} = consensus_H;
allSil{s} = Sil;
allE{s} = E;
allR{s} = R;
allCL{s} = CL;
end
save 'NMF_results_separate.mat' allHs allHc allSil allE allR allCL
%% find best ks
best_k = zeros(length(sample_id),1);
th = 0.9; % threshold in reproducibility to reach to select k
for s = 1:length(sample_id)
kid = find(kvals==best_k(s));
E = allE{s};
R = allR{s};
Sil = allSil{s};
CL = allCL{s};
allH = allHs{s};
kid = max(find(R>=th));
best_k(s) = kvals(kid);
figure;
subplot(1,2,1)
%subplot(length(sample_id),2,2*s-1)
yyaxis left
plot(kvals,E)
yyaxis right
plot(kvals,R)
hold on
plot(kvals,Sil)
%legend({'E','Reproducibility of distinct clusters','Silhouette score'})
y = tsne(allH{kid});
subplot(1,2,2)
gscatter(y(:,1),y(:,2),CL{kid},lines)
legend('off')
title(sprintf('%s',sample_id{s}))
print('-depsc2',sprintf('%s.eps',sample_id{s}))
end
%% get W from H
H = [];
sig_id = [];
for s = 1:length(sample_id)
data = logm(sample==s,:);
consensus_H = allHc{s};
kid = find(kvals==best_k(s));
H0 = consensus_H{kid};
W0 = max(0,data/H0);
opt = statset('Maxiter',1e6,'Display','final');
[~,H1] = nnmf(data,best_k(s),'h0',H0,'w0',W0,...
'options',opt,'replicates',100,...
'algorithm','als');
H = [H; H1];
sig_id = [sig_id; ones(size(H1,1),1)*s];
end
save all_signatures.mat H sig_id
相对还是很简单的，大家不妨试一试，找一找细胞类型相关协同性很高的gene set。
生活很好，有你更好
