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)

Wiki Page: Spatial Transcriptomics Cell Merging

Overview

The provided script is engineered to process spatial transcriptomics data. Its main goal is to merge cells based on polygon shapes defined in a JSON file. The merging is performed for each cell that resides within these polygons.

Dependencies

• numpy: For numerical operations and data structures.
• scanpy: A library to handle single-cell genomics data.
• json: For reading and processing JSON files.
• pandas: For data manipulation and structure.
• shapely: Provides geometric objects and operations, specifically used for polygons.
• scipy: Handles sparse matrix operations.
• rtree: Provides efficient spatial indexing.
• os: Interacts with the operating system, especially for file and directory operations.

Key Functions

1. load_h5ad(filename):

   • Input: Path to the .h5ad file.
   • Output: An AnnData object.
   • Description: Loads an .h5ad file which is typically used in single-cell genomics.

2. load_json_shapes(json_file_path):

   • Input: Path to the JSON file.
   • Output: A list of shapes defined in the JSON. Each shape contains a list of points.
   • Description: Reads a JSON file and retrieves the polygon shapes.

3. merge_cells_based_on_polygons(adata, shapes):

   • Input: An AnnData object and a list of shapes.
   • Output: A new AnnData object with merged cells.
   • Description: Merges cells based on polygons. If a cell's coordinates fall within a polygon, its data is merged. The function utilizes an R-tree spatial index for efficient querying.

4. get_prefix_from_filename(filename):

   • Input: Filename string.
   • Output: Extracted prefix of the filename.
   • Description: Strips out the .cellbin_v2 from the filename to get the prefix.

5. process_files_in_directory(directory):

   • Input: Path to the directory containing .h5ad and .json files.
   • Description: Processes all .h5ad files in the specified directory. For each .h5ad file, it searches for a corresponding .json file containing polygon shapes. Cells in the .h5ad file are then merged based on these polygons, and the result is saved as a new .h5ad file.

Execution Flow

The script starts by setting the directory variable to a path containing .h5ad and .json files. The process_files_in_directory(directory) function is then called to process all the files in the set directory. For each .h5ad file, the script:

• Searches for a corresponding .json file.
• Loads the AnnData object from the .h5ad file and the shapes from the .json file.
• Merges the cells based on the polygons.
• Saves the merged data into a new .h5ad file.

Conclusion

This script is a comprehensive solution for researchers working with spatial transcriptomics data who wish to merge cells based on predefined polygon shapes. It leverages efficient spatial querying using R-trees to make the process fast and efficient.