Generate Cell Bin (including batch processing)

installation

conda create --name scanpy
conda activate scanpy
conda install pip
pip install scanpy
pip install igraph
pip install leidenalg

pip install ipykernel

python -m ipykernel install --user --name=scanpy
pip install jupyter

Use

jupyter notebook

pip install shapely
pip install rtree

Here is the example code to cut the cell bin based on identified nuclei center and Voronoi segmentation.

import numpy as np
import scanpy as sc
import json
import pandas as pd
from shapely.geometry import Point, Polygon
from scipy.sparse import vstack
from rtree import index

def load_h5ad(filename):
    return sc.read(filename)

def load_json_shapes(json_file_path):
    with open(json_file_path, 'r') as file:
        json_data = json.load(file)
    return json_data.get('shapes', [])

def merge_cells_based_on_polygons(adata, shapes):
    new_data = []
    new_obs = []
    
    # Resetting x and y coordinates
    x_coords = adata.obs['x'].values
    y_coords = adata.obs['y'].values

    min_x = int(x_coords.min())
    min_y = int(y_coords.min())

    x_coords -= min_x
    y_coords -= min_y

    # Create an R-tree index for cells
    idx = index.Index()
    for i, (x, y) in enumerate(zip(x_coords, y_coords)):
        idx.insert(i, (x, y, x, y))
    
    total_shapes = len(shapes)
    for shape_num, shape in enumerate(shapes, 1):
        print(f"Processing shape {shape_num}/{total_shapes}...")
        
        polygon_coords = shape['points']
        polygon = Polygon(polygon_coords)
        
        merged_cell_data = None
        for i in idx.intersection(polygon.bounds):
            x, y = x_coords[i], y_coords[i]
            point = Point(x, y)
            if polygon.contains(point):
                cell_data = adata.X[i]
                if merged_cell_data is None:
                    merged_cell_data = cell_data
                else:
                    merged_cell_data += cell_data
        
        if merged_cell_data is not None:
            new_data.append(merged_cell_data)
            new_obs.append({
                'x': np.mean([p[0] for p in polygon_coords]),
                'y': np.mean([p[1] for p in polygon_coords]),
                'polygon': json.dumps(polygon_coords)  # Store polygon coordinates as a JSON string
            })

    if not new_data:
        print("No data to merge. Exiting.")
        return None

    new_X = vstack(new_data)
    new_adata = sc.AnnData(new_X)
    new_adata.var_names = adata.var_names  # Retain original gene names
    new_adata.obs = pd.DataFrame(new_obs)
    return new_adata

filename = r'G:\My Drive\Project\[1] Spinal Cord Spatial Transcriptome\h5ad\Bin1_visulization/L1_ur_B01204F1.cellbin.h5ad'
adata = load_h5ad(filename)

json_file_path = r'G:\My Drive\Project\[1] Spinal Cord Spatial Transcriptome\h5ad\Bin1_visulization/L1_ur_B01204F1.cellbin.h5ad_gene_counts.json'
shapes = load_json_shapes(json_file_path)

new_adata = merge_cells_based_on_polygons(adata, shapes)

if new_adata is not None:
    # Save the new adata object
    new_filename = filename.replace('.h5ad', '_merged_voronoi_v3.h5ad')
    new_adata.write(new_filename)
    print(f"New AnnData object saved to {new_filename}")
else:
    print("No AnnData object created.")

batch processing script

import numpy as np
import scanpy as sc
import json
import pandas as pd
from shapely.geometry import Point, Polygon
from scipy.sparse import vstack
from rtree import index
import os

def load_h5ad(filename):
    return sc.read(filename)

def load_json_shapes(json_file_path):
    with open(json_file_path, 'r') as file:
        json_data = json.load(file)
    return json_data.get('shapes', [])

def merge_cells_based_on_polygons(adata, shapes):
    new_data = []
    new_obs = []

    # Resetting x and y coordinates
    x_coords = adata.obs['x'].values
    y_coords = adata.obs['y'].values

    min_x = int(x_coords.min())
    min_y = int(y_coords.min())
    
    #reset the coordiante of h5ad files
    x_coords -= min_x
    y_coords -= min_y
    
    # Create an R-tree index for cells
    idx = index.Index()
    for i, (x, y) in enumerate(zip(x_coords, y_coords)):
        idx.insert(i, (x, y, x, y))
    
    total_shapes = len(shapes)
    for shape_num, shape in enumerate(shapes, 1):
        
        # Skip shapes with empty points
        if not shape['points']:
            print(f"Skipping shape {shape_num}/{total_shapes} due to empty points...")
            continue

        print(f"Processing shape {shape_num}/{total_shapes}...")
        
        polygon_coords = shape['points']
        polygon = Polygon(polygon_coords)
        
        # Ensure bounding box coordinates are in the right order
        minx, miny, maxx, maxy = polygon.bounds
        bounds = (min(minx, maxx), min(miny, maxy), max(minx, maxx), max(miny, maxy))
        
        merged_cell_data = None
        for i in idx.intersection(bounds):
            x, y = x_coords[i], y_coords[i]
            point = Point(x, y)
            if polygon.contains(point):
                cell_data = adata.X[i]
                if merged_cell_data is None:
                    merged_cell_data = cell_data
                else:
                    merged_cell_data += cell_data
        
        if merged_cell_data is not None:
            new_data.append(merged_cell_data)
            new_obs.append({
                'x': np.mean([p[0] for p in polygon_coords]),
                'y': np.mean([p[1] for p in polygon_coords]),
                'polygon': json.dumps(polygon_coords)  # Store polygon coordinates as a JSON string
            })

    if not new_data:
        print("No data to merge. Exiting.")
        return None

    new_X = vstack(new_data)
    new_adata = sc.AnnData(new_X)
    new_adata.var_names = adata.var_names  # Retain original gene names
    new_adata.obs = pd.DataFrame(new_obs)
    return new_adata



def get_prefix_from_filename(filename):
    # Return the prefix before .cellbin_v2
    return filename.split('.cellbin_v2')[0]

def process_files_in_directory(directory):
    # List all .h5ad files in the directory
    h5ad_files = [f for f in os.listdir(directory) if f.endswith('.h5ad')]

    for h5ad_file in h5ad_files:
        prefix = get_prefix_from_filename(h5ad_file)
        # Search for a matching JSON file
        matching_json_files = [f for f in os.listdir(directory) if f.startswith(prefix) and f.endswith('.json')]
        
        if matching_json_files:
            json_file = matching_json_files[0]
            print(f"Processing {h5ad_file} and {json_file}...")
            
            adata = load_h5ad(os.path.join(directory, h5ad_file))
            shapes = load_json_shapes(os.path.join(directory, json_file))
            
            new_adata = merge_cells_based_on_polygons(adata, shapes)
            
            if new_adata is not None:
                # Save the new adata object
                new_filename = h5ad_file.replace('.h5ad', '_merged_voronoi.h5ad')
                new_adata.write(os.path.join(directory, new_filename))
                print(f"New AnnData object saved to {new_filename}")
            else:
                print("No AnnData object created for", h5ad_file)
        else:
            print(f"No matching JSON found for {h5ad_file}. Skipping.")

directory = r'G:\My Drive\Project\[1] Spinal Cord Spatial Transcriptome\h5ad\section_h5ad_files_v2'

process_files_in_directory(directory)