Visium HD使用Cellpose进行细胞分割获取单细胞精度空间数据

# 导入必要的库import numpy as npimport pandas as pdimport scanpy as scimport tifffileimport geopandas as gpdfrom shapely.geometry import Pointfrom cellpose import models, core
# 定义读取图像的函数，并可选择下采样因子def read_dapi_image(img_path: str, downscale_factor: int = 2) -> np.ndarray:    img = io.imread(img_path)  # 读取图像文件    print(img.shape)  # 打印图像的形状    return downscale_image(img, downscale_factor=downscale_factor)  # 返回下采样后的图像
# 定义下采样图像的函数def downscale_image(img: np.ndarray, downscale_factor: int = 2) -> np.ndarray:    # 计算每个轴上所需的填充量    pad_height = (downscale_factor - img.shape[0] % downscale_factor) % downscale_factor    pad_width = (downscale_factor - img.shape[1] % downscale_factor) % downscale_factor    new_shape = (img.shape[0] + pad_height, img.shape[1] + pad_width, img.shape[2])  # 新图像形状    new_image = np.zeros(new_shape, dtype=img.dtype)  # 创建新的零填充图像    for channel in range(img.shape[2]):        new_image[:, :, channel] = np.pad(            img[:, :, channel], ((0, pad_height), (0, pad_width)), mode="constant"        )  # 对每个通道进行零填充    return new_image  # 返回填充后的图像
# 定义运行Cellpose进行细胞分割的函数def run_cellpose(img: np.ndarray, model_path: str):    use_GPU = core.use_gpu()  # 检查是否使用GPU    model = models.CellposeModel(gpu=use_GPU, pretrained_model=model_path)  # 加载Cellpose模型    channels = [0, 0]  # 设置通道（这里假设为单通道）    masks, flows, styles = model.eval(        [img],        channels=channels,        diameter=model.diam_labels,        flow_threshold=1,        cellprob_threshold=0,        batch_size=4,    )  # 执行细胞分割并返回掩膜、流和样式信息    return (masks, flows, styles)  # 返回分割结果
# 设置Visium HD数据的文件路径并读取图像img_path = "./data/visiumHD_colorectal_cancer/P1_CRC/Visium_HD_Human_Colon_Cancer_P1_tissue_image.btf"maxed_visium = read_dapi_image(img_path, downscale_factor=1)  # 读取并下采样图像
# 设置分块参数和模型路径chunk_per_axis = 2  # 每个轴上分块数目mp = 'cyto'  # 使用细胞质模型进行分割
# 对图像进行分块并调用Cellpose进行细胞分割masks_top_left, flows, styles = run_cellpose(    maxed_visium[        : np.shape(maxed_visium)[0] // chunk_per_axis,        : np.shape(maxed_visium)[1] // chunk_per_axis,        :,    ],    model_path=mp,)masks_top_right, flows, styles = run_cellpose(    maxed_visium[        : np.shape(maxed_visium)[0] // chunk_per_axis,        np.shape(maxed_visium)[1] // chunk_per_axis :,        :,    ],    model_path=mp,)masks_bottom_left, flows, styles = run_cellpose(    maxed_visium[        np.shape(maxed_visium)[0] // chunk_per_axis :,        : np.shape(maxed_visium)[1] // chunk_per_axis,        :,    ],    model_path=mp,)masks_bottom_right, flows, styles = run_cellpose(    maxed_visium[        np.shape(maxed_visium)[0] // chunk_per_axis :,        np.shape(maxed_visium)[1] // chunk_per_axis :,        :,    ],    model_path=mp,)
# 初始化全掩膜数组，准备合并分块结果constant = 1000000  # 用于区分不同块的常量值full_mask = np.zeros_like(maxed_visium[:, :, 0], dtype=np.uint32)  # 创建与原图相同大小的掩膜数组
# 将各个块的掩膜合并到全掩膜中，并加上常量以区分不同区域full_mask[    : np.shape(maxed_visium)[0] // chunk_per_axis,    : np.shape(maxed_visium)[1] // chunk_per_axis,] = masks_top_left[0]full_mask[    : np.shape(maxed_visium)[0] // chunk_per_axis,    np.shape(maxed_visium)[1] // chunk_per_axis :,] = masks_top_right[0] + (constant)full_mask[    np.shape(maxed_visium)[0] // chunk_per_axis :,    : np.shape(maxed_visium)[1] // chunk_per_axis,] = masks_bottom_left[0] + (constant * 2)full_mask[    np.shape(maxed_visium)[0] // chunk_per_axis :,    np.shape(maxed_visium)[1] // chunk_per_axis :,] = masks_bottom_right[0] + (constant * 3)
# 将全掩膜中值为常量的部分设置为0，表示没有细胞full_mask = np.where(full_mask % constant == 0, 0, full_mask)
# 将全掩膜保存为PNG格式文件tifffile.imsave("./visium_hd/segmentation.png", full_mask)# 设置Visium HD数据目录路径以加载空间转录组数据dir_base = "./data/visiumHD_colorectal_cancer/P1_CRC/outs/binned_outputs/square_002um/"
# 加载Visium HD数据中的基因表达矩阵（HDF5格式）raw_h5_file = dir_base + "filtered_feature_bc_matrix.h5"adata = sc.read_10x_h5(raw_h5_file)  # 使用Scanpy读取数据
# 加载组织位置坐标文件（Parquet格式）tissue_position_file = dir_base + "spatial/tissue_positions.parquet"df_tissue_positions = pd.read_parquet(tissue_position_file)
# 将数据框的索引设置为条形码，以便后续合并操作使用df_tissue_positions = df_tissue_positions.set_index("barcode")
# 在数据框中创建索引以检查合并情况df_tissue_positions["index"] = df_tissue_positions.index
# 将组织位置添加到AnnData对象的元数据中（obs）adata.obs = pd.merge(adata.obs, df_tissue_positions, left_index=True, right_index=True)
# 从坐标数据框创建GeoDataFrame，以便进行空间分析geometry = [    Point(xy)    for xy in zip(        df_tissue_positions["pxl_col_in_fullres"],        df_tissue_positions["pxl_row_in_fullres"],    )]gdf_coordinates = gpd.GeoDataFrame(df_tissue_positions, geometry=geometry)
## 将捕获区域分配给单个细胞# 检查点掩膜是否在全掩膜范围内，以确保坐标有效性point_mask = (    gdf_coordinates["pxl_row_in_fullres"].values.astype(int) < np.shape(full_mask)[0]) & (gdf_coordinates["pxl_col_in_fullres"].values.astype(int) < np.shape(full_mask)[1])
# 根据有效坐标从全掩膜中提取细胞ID信息cells = full_mask[    gdf_coordinates["pxl_row_in_fullres"].values.astype(int)[point_mask],    gdf_coordinates["pxl_col_in_fullres"].values.astype(int)[point_mask],]gdf_coordinates = gdf_coordinates[point_mask]gdf_coordinates["cells"] = cells  # 将细胞ID添加到GeoDataFrame中assigned_regions = gdf_coordinates[gdf_coordinates["cells"] > 0]
# 将细胞ID信息合并到AnnData对象中，以便后续分析使用adata.obs = adata.obs.merge(    assigned_regions[["cells"]], left_index=True, right_index=True, how="left")adata = adata[~pd.isna(adata.obs["cells"])]  # 移除没有细胞ID的数据行adata.obs["cells"] = adata.obs["cells"].astype(int)  # 将细胞ID转换为整数类型
## 在捕获区域内汇总转录本计数from scipy import sparse
# 按唯一细胞ID对数据进行分组，以便汇总转录本计数groupby_object = adata.obs.groupby(["cells"], observed=True)counts = adata.X  # 提取基因表达计数矩阵spatial_coords = adata.obs[["pxl_row_in_fullres", "pxl_col_in_fullres"]].valuesN_groups = groupby_object.ngroups  # 获取唯一细胞ID数量（组数）N_genes = counts.shape[1]  # 获取基因数量
# 初始化稀疏矩阵以存储每个细胞的基因计数汇总结果summed_counts = sparse.lil_matrix((N_groups, N_genes))polygon_id = []  # 存储多边形ID列表cell_coords = []  # 存储细胞坐标列表row = 0
# 遍历每个唯一多边形，计算基因计数的总和。for polygons, idx_ in groupby_object.indices.items():    summed_counts[row] = counts[idx_].sum(0)  # 汇总该多边形内所有细胞的基因计数    cell_coords.append(np.mean(spatial_coords[idx_], axis=0))  # 获取该多边形内细胞坐标均值作为代表坐标    row += 1    polygon_id.append(polygons)  # 添加多边形ID到列表中cell_coords = np.array(cell_coords)  summed_counts = summed_counts.tocsr()  # 转换为压缩稀疏行矩阵格式
# 创建一个新的AnnData对象，用于存储汇总后的基因计数矩阵及其元数据。grouped_adata = anndata.AnnData(    X=summed_counts,    obs=pd.DataFrame(polygon_id, columns=["cells"], index=polygon_id),    var=adata.var,)grouped_adata.write('adata_cellpose.h5ad')