adding simple spatial capability

Author: scyrusm

panopticon/spatial.py

@@ -0,0 +1,48 @@
+import pandas as pd
+import numpy as np
+import matplotlib.pyplot as plt
+
+def get_labels_and_polys(tif_file):
+    from panopticon.utilities import import_check
+    exit_code = import_check(
+        "tifffile", 'pip install tifffile')
+    if exit_code != 0:
+        return
+    exit_code = import_check(
+        "stardist", 'pip install stardist')
+    if exit_code != 0:
+        return
+    from tifffile import imread, imwrite, imshow
+    from stardist.models import StarDist2D
+    img = imread(tif_file)
+    img = img.swapaxes(0,2) #Assumes to be CXY, we switch to YXC
+    model = StarDist2D.from_pretrained('2D_versatile_he')
+    
+    img = img/255
+    labels, polys = 
+    labels, polys = model.predict_instances_big(img, axes='YXC', block_size=4096, prob_thresh=0.01,nms_thresh=0.001, min_overlap=128, context=128, normalizer=None, n_tiles=(4,4,1))
+
+def get_gdf(labels, polys):
+    exit_code = import_check(
+        "geopandas", 'pip install geopandas')
+    if exit_code != 0:
+        return
+    import geopandas as gpd
+
+    # Creating a list to store Polygon geometries
+    geometries = []
+    
+    # Iterating through each nuclei in the 'polys' DataFrame
+    for nuclei in range(len(polys['coord'])):
+    
+        # Extracting coordinates for the current nuclei and converting them to (y, x) format
+        coords = [(y, x) for x, y in zip(polys['coord'][nuclei][0], polys['coord'][nuclei][1])]
+    
+        # Creating a Polygon geometry from the coordinates
+        geometries.append(Polygon(coords))
+    
+    # Creating a GeoDataFrame using the Polygon geometries
+    gdf = gpd.GeoDataFrame(geometry=geometries)
+    gdf['id'] = [f"ID_{i+1}" for i, _ in enumerate(gdf.index)]
+    return gdf
+    

scripts/panopticon

@@ -45,19 +45,26 @@ def combine_looms(input_looms, output_loom, output_type='loom'):
     'Whitelist a gene to include (will blacklist all other genes).  Intended for getting antibody counts from empty droplets, where other types of counts are not needed.  Usage --whitelist_gene Hash123',
     multiple=True,
     default=None)
+@click.option(
+    '--exclude_feature_type',
+    help=
+    'Exclude a feature type from conversion (e.g. "Gene Expression", "Peaks", etc.).  Usage --exclude_feature_type Peaks',
+    multiple=True,
+    default=None)
 @click.option(
     '--write_chunked',
     help='flag to write output loom in chunked fashion (to reduce memory burden)',
     is_flag=True)
-def convert_10x_h5(h5file, output_loom, label, convert_to_ca, whitelist_gene,write_chunked):
+def convert_10x_h5(h5file, output_loom, label, convert_to_ca, whitelist_gene,write_chunked, exclude_feature_type):
     from panopticon.utilities import convert_10x_h5
     convert_10x_h5(h5file,
                    output_loom,
                    label[0],
                    label[1],
                    genes_as_ca=convert_to_ca,
                    gene_whitelist=whitelist_gene,
-                   write_chunked=write_chunked)
+                   write_chunked=write_chunked,
+                   exclude_feature_type=exclude_feature_type)
 
 
 @cli.command(
@@ -136,6 +143,25 @@ def reaptec(fastq_dir, cellranger_output_dir, star_reference_dir,
     reaptec_main(fastq_dir, cellranger_output_dir, star_reference_dir,
                  reference_url, genome_url)
 
+@cli.command(
+        help="Generating a geopandas dataframe with location of called nuclei")
+@click.argument('tif_file', type=click.Path(exists=True),nargs=1)
+@click.argument('output_file_name',nargs=1)
+def generate_segmentation_gdf(tif_file, output_file_name):
+    from panopticon.spatial import get_labels_and_polys, get_gdf
+    if output_file_name.split('.')[-1] not in ['pkl','csv','parquet']:
+        raise Exception("output_file_name should end in one of ['csv','pkl','parquet']")
+    labels, polys = get_labels_and_polys(tif_file)
+    gdf = get_gdf(labels, polys)
+    if output_file_name.endswith('.pkl'):
+        gdf.to_pickle(output_file_name)
+    elif output_file_name.endswith('.csv'):
+        gdf.to_csv(output_file_name)
+    elif output_file_name.endswith('.parquet'):
+        gdf.to_parquet(output_file_name)
+
+
+
 
 if __name__ == '__main__':
     cli()