细胞空间结构特征之-细胞集群特征（Matlab实现）

本文的目的是根据细胞的分割图像计算一些特征用于描述细胞的集群情况，细胞集群特征是细胞空间结构特征的一种，使用的编程工具为Matlab 2018b。

本文中使用例子为TCGA细胞分割数据集中的一张，原始图像（尺寸为1000×1000）如下左，相应的细胞分割图（二值图）如下右，计算细胞集群特征使用的是右边的细胞分割图，原始图像不参与计算。

 首先根据细胞分割图计算一些提取细胞集群特征需要用到的变量：

% 读取细胞分割图像
mask = imread('C:\Users\Administrator\Desktop\mask.png');
% 将细胞分割图像转换为0-1二值图
mask = imbinarize(mask);
% 进行孔洞填充处理
mask = imfill(mask, 'hole');
% 进行连通域标号操作
L = bwlabel(mask);
% 计算每个细胞（白色连通域）的中心点坐标
C = regionprops(L, 'centroid');
% 变量centroids是一个行数为连通域个数，列数为2的矩阵
% 第一列存放的是每个细胞中心点的横坐标
% 第二列存放的是每个细胞中心点的纵坐标
centroids = cat(1, C.Centroid);
% 变量N存放的是图中所有细胞（白色连通域）的个数
N = length(C);
% 变量D存放的是所有细胞中心点两两之间的欧氏距离
% 变量D中共有N(N-1)/2个元素，每个元素都代表一个距离
D = pdist(centroids, 'euclidean');
% 变量cluster是一个尺寸为N×N的矩阵
% 若cluster(m,n)==1，则表示第m个细胞和第n个细胞中心点的距离小于50（像素），称它们为“非孤立细胞”
% (m,n)和(n,m)不重复计算，即1值只存在于变量cluster的上半三角形
cluster= triu(true(N), 1) & squareform(D < 50);
% 变量cluster_all是将矩阵cluster扩展到下半三角后的结果，即将(m,n)和(n,m)重复考虑
cluster_all = cluster | cluster';
% 变量cluster_num中存放的是非孤立细胞的个数
[row, col, V] = find(cluster_all ~= 0);
cluster_num = length(unique(row));
% 将变量D改写为矩阵形式
% D(m,n)表示第m个细胞与第n个细胞中心点之间的距离
D = squareform(D);
% 变量cluster_D是稀疏矩阵，保存的是所有非孤立细胞两两之间的序号和中心点距离
cluster_D = sparse(D.*cluster_all);
% 变量distances是一个尺寸为1×N的元胞数组
% 若第i个细胞是孤立细胞，则distances(i)为空
% 若第i个细胞是非孤立细胞，则distances(i)中保存的是第i个细胞到它邻居，和它邻居的邻居，和它邻居的 
% 邻居的邻居等等等等等的最短路径长度
for i = 1:N
    distance = graphshortestpath(cluster_D, i);
    distances{i} = nonzeros(distance(isfinite(distance))).';
end
% 变量distances_noempty中存放的是非孤立细胞到它邻居们的最短路径
distances_noempty = distances;
distances_noempty(cellfun(@isempty, distances_noempty, 'UniformOutput', true)) = [];
% 变量distances_max中存放的是每个非孤立细胞到邻居们最短路径的最大值
distances_max = cellfun(@max, distances_noempty, 'UniformOutput', true);
% 变量network中存放的是一个N维向量，表示每个细胞属于第几个小集群
[~, network] = graphconncomp(sparse(cluster_all));
% 变量nodes中存放的是第n个小集群中细胞的序号
% 变量En中存放的是每个细胞所属的小集群共有多少条连接通路
% 变量Kn中存放的是每个细胞所属的小集群中有多少细胞
for n = 1:N
    nodes = find(network == n);
    En(nodes) = sum(sum(cluster(nodes, nodes)));
    Kn(nodes) = length(nodes);
end

根据上面计算到的变量，结合统计学方法设计一些特征，用于表示细胞的集群特点：

% 定义第一维特征为图像中的细胞的数目（白色连通域个数）
ccfeats(1) = N;
% 定义第二维特征为图像中孤立细胞的数目
ccfeats(2) = N - cluster_num;
% 定义第三维特征为图像中孤立细胞的占比
ccfeats(3) = (N - cluster_num) / N;
% 定义第四维特征为非孤立细胞到邻居们距离最大值的平均
ccfeats(4) = sum(distances_max) / cluster_num;
% 定义第五维特征为非孤立细胞到邻居们距离最大值的最大值
ccfeats(5) = max(distances_max);
% 定义第六维特征为非孤立细胞到邻居们距离最大值的最小值
ccfeats(6) = min(distances_max);
% 定义第七维特征为非孤立细胞到邻居们距离的平均值
ccfeats(7) = sum(cellfun(@sum, distances, 'UniformOutput', true)) / sum(cellfun(@length, distances, 'UniformOutput', true));
% 定义第八维特征为一种计算集群系数的方法
Cn = 2 * En ./ (Kn .* (Kn - 1));
Cn(isnan(Cn)) = 0;
ccfeats(8) = sum(Cn) / N;
% 定义第九维特征为另一种计算集群系数的方法
Dn = 2 * (Kn + En) ./ (Kn .* (Kn + 1));
Dn(isnan(Dn)) = 0;
ccfeats(9) = sum(Dn) / N;
% 定义第十维特征为第三种计算集群系数的方法
% 变量iso_nodes中存放的是孤立细胞的个数
iso_nodes = sum(Kn == 1);
if iso_nodes == N
    ccfeats(10) = NaN;
    fprintf('所有细胞都是孤立细胞');
else
    ccfeats(10) = sum(Cn(Kn > 1)) / (N - iso_nodes);
end
% 定义第十一维特征为非孤立细胞的个数
ccfeats(11) = sum(Kn > 1);
% 定义第十二维特征为最大的细胞集群中细胞个数占细胞总个数的比
ccfeats(12) = max(Kn) / N;
% 定义第十三维特征为所有细胞集群中细胞个数的平均
if N == iso_nodes
    ccfeats(13) = 1;
else
    ccfeats(13) = mean(Kn(Kn > 1));
end
% 定义第十四维特征为孤立细胞的个数
ccfeats(14) = iso_nodes;
% 定义第十五维特征为孤立细胞个数占总细胞个数的比
ccfeats(15) = iso_nodes / N;
% 定义第十六维特征为只含有两个细胞的集群个数
ccfeats(16) = sum(Kn == 2);
% 定义第十七维特征为两细胞集群中的细胞个数占总细胞个数的比
ccfeats(17) = sum(Kn == 2) / N;
% 定义第十八维特征为非孤立细胞到邻居们距离最大值中有几个细胞取到最小值
ccfeats(18) = sum(distances_max == min(distances_max));
% 定义第十九维特征为非孤立细胞到邻居们距离最大值中取到最小值的细胞个数占总细胞个数的比
ccfeats(19) = sum(distances_max == min(distances_max)) / N;
% 定义第二十维特征为非孤立细胞之间距离的平均值
% 定义第二十一维特征为非孤立细胞之间距离的方差
% 定义第二十二维特征为非孤立细胞之间距离的偏度
% 定义第二十三维特征为非孤立细胞之间距离的峰度
% 变量edge_lengths中存放的是非孤立细胞两两之间的距离
if N == iso_nodes
    ccfeats(20) = 0;
    ccfeats(21) = 0;
    ccfeats(22) = 0;
    ccfeats(23) = 0;
else
    edge_lengths = D .* cluster_all;
    edge_lengths(edge_lengths == 0) = [];
    ccfeats(20) = sum(edge_lengths) / length(edge_lengths);
    ccfeats(21) = std(edge_lengths);
    ccfeats(22) = skewness(edge_lengths);
    ccfeats(23) = kurtosis(edge_lengths);
end