计算某个区域中细胞的复杂性(2)

上一篇帖子提到了猴脑文章中某个区域中细胞复杂性的计算，后面又看到一张图，才是发现自己的计算少了一些东西，这里补充一下。

在原文的图四中还有一个图片，讲述了较为详细的计算过程，如下图所示：

首先是复杂性的获得，我的上一篇相关的博客已经讲了，但是那篇博客的问题也是在于此，我是直接获得了复杂性，把其他的信息都丢掉了，KD树的缺陷也在于此，于是不得已在此进行重构：

# 首先是包的引入，这里可能有很多包都是不需要的，但是我懒得删了，就这样吧
import numpy as np
from scipy.spatial.distance import cdist
import anndata as ad
import pandas as pd
import matplotlib.pyplot as plt
from collections import Counter
import seaborn as sns
from matplotlib.ticker import FuncFormatter
from matplotlib.colors import ListedColormap
from operator import itemgetter

然后是计算函数的定义：

# Function to calculate the neighborhood complexity
def calculate_neighborhood_complexity(coordinates, cell_types):
    '''
    coordinates: 细胞的坐标 list
    cell_types: 细胞的类型 list
    '''
    num_cells = len(coordinates)
    unique_cell_types = np.unique(cell_types)
    num_cell_types = len(unique_cell_types)

    cell_type_indices = {cell_type: i for i, cell_type in enumerate(unique_cell_types)}

    complexity_matrix = np.zeros((num_cells, num_cell_types), dtype=int)

    for i in range(num_cells):
        center_cell_type = cell_types[i]
        center_coord = coordinates[i]

        distances = cdist([center_coord], coordinates)[0]
        neighborhood_indices = np.where(distances <= 200)[0]  # 100 mm = 200 pixels

        neighborhood_cell_types = cell_types[neighborhood_indices]

        unique_clusters, cluster_counts = np.unique(neighborhood_cell_types, return_counts=True)
        
        for cell_type, count in zip(unique_clusters, cluster_counts):
            column_index = cell_type_indices[cell_type]
            complexity_matrix[i, column_index] += count
    complexity_matrix = pd.DataFrame(complexity_matrix, columns=unique_cell_types, )
    complexity = (complexity_matrix != 0).sum(axis=1)
    complexity_matrix['complexity'] = complexity
    complexity_matrix['cell_types'] = cell_types
    return complexity_matrix

如何使用呢？如下所示：

adata = ad.read('/data/cellbin_transAnno.h5ad')
# adata.obsm['spatial'] 中存储的是每一个细胞的空间坐标信息
# adata.obs['celltype'] 中存储的是每一个细胞的细胞类型注释信息
coordinates = adata_temp.obsm['spatial']
cell_types = np.array(adata_temp.obs['celltype'].tolist())
complexity_matrix  = calculate_neighborhood_complexity(coordinates, cell_types)
complexity_matrix.to_csv('/data/complexity_matrix.csv')

这里的 complexity_matrix 是一个 pandas 的 DataFrame，不会用 python 画图的友友可以保存为 csv 文件，然后用 R 语言画图，本人不是很会使用 R 语言，所以就不提供 R 语言的代码了，dataframe 概览如下所示，每一行为一个细胞，每一列为一个细胞类型，最后两列为每个细胞的复杂性和细胞类型，中间的值为 distances 在 200 之内的 某个细胞类型的数量。

本人绘图的代码如下所示

adata_temp = adata
coordinates = adata_temp.obsm['spatial']
cell_types = np.array(adata_temp.obs['celltype'].tolist())

complexity_matrix  = calculate_neighborhood_complexity(coordinates, cell_types)

# 选取自己感兴趣的区域/细胞类型
ex = ['X1', 'X2', 'X3']



complexity_matrix = complexity_matrix[complexity_matrix['cell_types'].isin(ex)]

fig, ax = plt.subplots(figsize=(9, 6))

for j in ex:
    temp = complexity_matrix[complexity_matrix['cell_types'] == j]
    # category = temp['cell_types'].tolist()
    if len(temp)>0:

        neighborhood_complexities = temp['complexity'].tolist()


        sns.kdeplot(neighborhood_complexities, cumulative=True, fill=False, ax=ax, label=j, color = celltypes_colors[j])
    else:
        continue

ax.set_xlabel('Neighborhood Complexity')
ax.set_ylabel('Probability (%)')
ax.set_title('ALL Neighborhood Complexity Distribution')
# Position the legend on the outside right
ax.legend(loc='upper left', bbox_to_anchor=(1.02, 1))

ax.yaxis.set_major_formatter(FuncFormatter(lambda x, _: f'{x*100:.0f}'))
ax.xaxis.set_major_locator(plt.MaxNLocator(integer=True))
plt.tight_layout()
plt.savefig('Neighborhood_Complexity_Distribution.pdf')
plt.close()

绘图结果如下所示：

以上