create_sample

Author: EliHei2

scripts/create_data_fast_sample.py

@@ -2,12 +2,16 @@
 from path import Path
 from segger.data.utils import calculate_gene_celltype_abundance_embedding
 import scanpy as sc
+import pandas as pd
+import math
+import numpy as np
+from segger.data.parquet._utils import get_polygons_from_xy
 
 xenium_data_dir = Path('data_raw/breast_cancer/Xenium_FFPE_Human_Breast_Cancer_Rep1/outs/')
-segger_data_dir = Path('data_tidy/pyg_datasets/bc_fast_data_emb_major')
+segger_data_dir = Path('data_tidy/pyg_datasets/bc_rep1_emb')
 
 
-scrnaseq_file = Path('data_tidy/benchmarks/xe_rep1_bc/scRNAseq.h5ad')
+scrnaseq_file = Path('/omics/groups/OE0606/internal/tangy/tasks/schier/data/atals_filtered.h5ad')
 celltype_column = 'celltype_major'
 gene_celltype_abundance_embedding = calculate_gene_celltype_abundance_embedding(
     sc.read(scrnaseq_file),
@@ -21,16 +25,69 @@
     weights=gene_celltype_abundance_embedding, # uncomment if gene-celltype embeddings are available
 )
 
+transcripts = pd.read_parquet(
+    xenium_data_dir / 'transcripts.parquet',
+    filters=[[('overlaps_nucleus', '=', 1)]]
+)
+boundaries = pd.read_parquet(xenium_data_dir / 'nucleus_boundaries.parquet')
+
+sizes = transcripts.groupby('cell_id').size()
+polygons = get_polygons_from_xy(boundaries, 'vertex_x', 'vertex_y', 'cell_id')
+densities = polygons[sizes.index].area / sizes
+bd_width = polygons.minimum_bounding_radius().median() * 2
+
+# 1/4 median boundary diameter
+dist_tx = bd_width / 4
+# 90th percentile density of bounding circle with radius=dist_tx
+k_tx = math.ceil(np.quantile(dist_tx ** 2 * np.pi * densities, 0.9)) 
+
+print(k_tx)
+print(dist_tx)
+
+
 sample.save(
       data_dir=segger_data_dir,
       k_bd=3,
       dist_bd=15.0,
-      k_tx=20,
-      dist_tx=3,
-      tile_width=220,
-      tile_height=220,
+      k_tx=dist_tx,
+      dist_tx=k_tx,
+      tile_width=120,
+      tile_height=120,
       neg_sampling_ratio=5.0,
       frac=1.0,
       val_prob=0.1,
       test_prob=0.1,
-)
+)
+
+
+xenium_data_dir = Path('data_tidy/bc_5k')
+segger_data_dir = Path('data_tidy/pyg_datasets/bc_5k_emb')
+
+
+
+sample = STSampleParquet(
+    base_dir=xenium_data_dir,
+    n_workers=1,
+    sample_type='xenium',
+    weights=gene_celltype_abundance_embedding, # uncomment if gene-celltype embeddings are available
+)
+
+
+transcripts = pd.read_parquet(
+    xenium_data_dir / 'transcripts.parquet',
+    filters=[[('overlaps_nucleus', '=', 1)]]
+)
+boundaries = pd.read_parquet(xenium_data_dir / 'nucleus_boundaries.parquet')
+
+sizes = transcripts.groupby('cell_id').size()
+polygons = get_polygons_from_xy(boundaries, 'vertex_x', 'vertex_y', 'cell_id')
+densities = polygons[sizes.index].area / sizes
+bd_width = polygons.minimum_bounding_radius().median() * 2
+
+# 1/4 median boundary diameter
+dist_tx = bd_width / 4
+# 90th percentile density of bounding circle with radius=dist_tx
+k_tx = math.ceil(np.quantile(dist_tx ** 2 * np.pi * densities, 0.9)) 
+
+print(k_tx)
+print(dist_tx)