细胞邻域分析

import argparse
import scanpy as sc
import numpy as np
import pandas as pd
from sklearn.cluster import AgglomerativeClustering
from sklearn.metrics.pairwise import cosine_similarity
import matplotlib.pyplot as plt
import os

def get_cells_within_radius(center, radius, positions):
    """
    获取某个中心点附近指定半径范围内的细胞
    :param center: 中心点坐标 (x, y)
    :param radius: 范围半径（单位：微米）
    :param positions: 细胞空间坐标
    :return: 半径范围内的细胞索引
    """
    distances = np.linalg.norm(positions - center, axis=1)
    return np.where(distances <= radius)[0]

def hierarchical_clustering(input_file, radius, output_dir):
    """
    使用层次聚类根据邻域细胞类型进行聚类
    :param input_file: 输入的AnnData文件路径（h5ad格式）
    :param radius: 邻域半径（单位：微米）
    :param output_dir: 输出目录，用于保存结果
    """
    # 加载空间转录组数据
    adata = sc.read(input_file)

    # 假设空间坐标在adata.obsm['spatial']中，细胞类型在adata.obs['cell_type']中
    positions = adata.obsm['spatial']  # 空间坐标
    cell_types = adata.obs['cell_type']  # 细胞类型信息

    # 获取肿瘤细胞索引
    tumor_cells_indices = adata.obs.query("cell_type == 'Tumor'").index

    # 提取肿瘤细胞邻域的所有细胞类型
    neighbor_cell_types = []
    for tumor_cell_idx in tumor_cells_indices:
        center = positions[tumor_cell_idx]
        nearby_cells = get_cells_within_radius(center, radius, positions)
        neighbor_cell_types.extend(cell_types[nearby_cells].values)

    # 计算邻域细胞类型的频率分布
    neighbor_cell_types = pd.Series(neighbor_cell_types)
    cell_type_counts = neighbor_cell_types.value_counts()

    # 创建细胞类型的邻接矩阵（协同频率矩阵）
    cell_type_list = cell_type_counts.index
    type_to_idx = {cell_type: idx for idx, cell_type in enumerate(cell_type_list)}
    adj_matrix = np.zeros((len(cell_type_list), len(cell_type_list)))

    # 填充邻接矩阵
    for tumor_cell_idx in tumor_cells_indices:
        center = positions[tumor_cell_idx]
        nearby_cells = get_cells_within_radius(center, radius, positions)
        nearby_cell_types = cell_types[nearby_cells].values
        for cell_type in nearby_cell_types:
            for other_cell_type in nearby_cell_types:
                if cell_type != other_cell_type:
                    i, j = type_to_idx[cell_type], type_to_idx[other_cell_type]
                    adj_matrix[i, j] += 1
                    adj_matrix[j, i] += 1  # 对称矩阵

    # 归一化邻接矩阵
    adj_matrix /= adj_matrix.sum(axis=1, keepdims=True)

    # 使用余弦相似度计算细胞类型的相似性
    similarity_matrix = cosine_similarity(adj_matrix)

    # 层次聚类
    clustering = AgglomerativeClustering(n_clusters=None, distance_threshold=0.5, affinity='precomputed', linkage='average')
    clustering.fit(1 - similarity_matrix)  # 因为我们使用余弦相似度，所以相似度越高，距离越小，取1减去相似度

    # 将聚类结果添加到adata中
    adata.obs['cluster'] = 'Unclassified'
    adata.obs.loc[tumor_cells_indices, 'cluster'] = clustering.labels_

    # 创建输出目录（如果不存在）
    if not os.path.exists(output_dir):
        os.makedirs(output_dir)

    # 输出聚类结果到输出目录
    output_file = os.path.join(output_dir, "clustered_cells.h5ad")
    adata.write(output_file)
    print(f"聚类结果已保存至 {output_file}")

    # 可视化聚类结果
    sc.pl.spatial(adata, color='cluster', title="Hierarchical Clustering Based on Cell Type", save=f"{output_dir}/spatial_cluster.png")

    # 可选：显示聚类树状图
    plt.figure(figsize=(10, 7))
    plt.title("Hierarchical Clustering Dendrogram")
    plt.xlabel("Cell Types")
    plt.ylabel("Distance")
    from scipy.cluster.hierarchy import dendrogram
    dendrogram(clustering.children_)
    plt.savefig(f"{output_dir}/dendrogram.png")
    plt.close()

def main():
    # 使用argparse处理命令行参数
    parser = argparse.ArgumentParser(description="根据细胞类型进行邻域层次聚类")
    parser.add_argument('input_file', type=str, help="输入AnnData文件（h5ad格式）")
    parser.add_argument('--radius', type=float, default=200, help="邻域半径（单位：微米），默认为200μm")
    parser.add_argument('--output_dir', type=str, required=True, help="输出目录，用于保存结果")

    # 解析命令行参数
    args = parser.parse_args()

    # 调用聚类函数
    hierarchical_clustering(args.input_file, args.radius, args.output_dir)

if __name__ == "__main__":
    main()