Skip to content

Commit 1f3edda

Browse files
JBeckwithJBeckwith
committed
Refactoring
1 parent 6253772 commit 1f3edda

19 files changed

Lines changed: 1449 additions & 168 deletions

‎Example_Database_Queries.ipynb‎

Lines changed: 553 additions & 0 deletions
Large diffs are not rendered by default.

‎oligomer_concentration_by_region.pdf‎

29 KB
Binary file not shown.

‎oligomer_concentration_by_region_and_celltype.pdf‎

33.9 KB
Binary file not shown.

‎oligomer_concentration_sex_by_region.pdf‎

89.6 KB
Binary file not shown.

‎scripts/build_database.py‎

Lines changed: 489 additions & 0 deletions
Large diffs are not rendered by default.

‎scripts/run_cell_analysis_for_db.py‎

Lines changed: 155 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,155 @@
1+
#!/usr/bin/env python3
2+
"""
3+
Re-run the cell coincidence analysis for all NewCellAnalysis directories to generate:
4+
1. Updated percentile-0 cell CSV with a cell_id_in_image column
5+
2. oligomers_in_cells_percentile_0.csv — one row per in-cell oligomer with cell_id_in_image
6+
7+
Usage:
8+
python scripts/run_cell_analysis_for_db.py --test # cingulate/microglia only
9+
python scripts/run_cell_analysis_for_db.py --full # all 38 analysis directories
10+
11+
Outputs are written alongside the existing NewCellAnalysis CSV files.
12+
Directories that already have oligomers_in_cells_percentile_0.csv are skipped.
13+
"""
14+
15+
import argparse
16+
import os
17+
import re
18+
import sys
19+
20+
import polars as pl
21+
22+
sys.path.insert(0, os.path.join(os.path.dirname(__file__), "..", "src"))
23+
from AnalysisFunctions import Analysis_Functions
24+
25+
DATA_ROOT = "/scratch/sycamore-asap/ASAP_Imaging_Data/Main_Survey/20241105_oligomers_in_cells"
26+
ANALYSIS_SUBDIR = "NewCellAnalysis_OnlyLowerThreshold_NewPhotonThreshold"
27+
28+
# Analysis parameters matching the original runs
29+
PHOTON_THRESHOLD = 150.0
30+
LOWER_CELL_SIZE = 5000
31+
UPPER_CELL_SIZE = 200000
32+
SPACING = (0.5, 0.11, 0.11)
33+
BLUR_DEGREE = 4
34+
CELL_STRING = "C0"
35+
PROTEIN_STRING = "C1"
36+
37+
P0_PATTERN = re.compile(
38+
r"single_cell_coincidence_.*_percentile_0_.*abovethreshold\.csv$"
39+
)
40+
41+
ALL_TARGETS = [
42+
("caudate", "microglia"),
43+
("caudate", "neurons"),
44+
("cingulate", "astrocytes"),
45+
("cingulate", "microglia"),
46+
("cingulate", "neurons"),
47+
("frontal", "astrocytes"),
48+
("frontal", "microglia"),
49+
("frontal", "neurons"),
50+
("parahippocampal", "microglia"),
51+
("parahippocampal", "neurons"),
52+
("parietal", "astrocytes"),
53+
("parietal", "microglia"),
54+
("parietal", "neurons"),
55+
("putamen", "astrocytes"),
56+
("putamen", "microglia"),
57+
("putamen", "neurons"),
58+
("substantiaNigra", "microglia"),
59+
("temporal", "microglia"),
60+
("temporal", "neurons"),
61+
]
62+
63+
64+
def process(af, region, cell_type, condition):
65+
analysis_dir = os.path.join(
66+
DATA_ROOT, region, cell_type, f"{condition}_analysis"
67+
)
68+
spot_csv = os.path.join(analysis_dir, "spot_analysis.csv")
69+
new_cell_dir = os.path.join(analysis_dir, ANALYSIS_SUBDIR)
70+
oligo_out = os.path.join(new_cell_dir, "oligomers_in_cells_percentile_0.csv")
71+
72+
if not os.path.isfile(spot_csv):
73+
print(f" SKIP no spot_analysis.csv: {analysis_dir}")
74+
return
75+
if not os.path.isdir(new_cell_dir):
76+
print(f" SKIP no NewCellAnalysis dir: {analysis_dir}")
77+
return
78+
if os.path.isfile(oligo_out):
79+
print(f" SKIP already complete: {oligo_out}")
80+
return
81+
82+
# Find the percentile-0 cell CSV to overwrite
83+
p0_candidates = [f for f in os.listdir(new_cell_dir) if P0_PATTERN.match(f)]
84+
if not p0_candidates:
85+
print(f" SKIP no percentile-0 cell CSV in: {new_cell_dir}")
86+
return
87+
cell_csv_out = os.path.join(new_cell_dir, p0_candidates[0])
88+
89+
print(f" Loading spots from {spot_csv} ...")
90+
spot_data = (
91+
pl.scan_csv(spot_csv)
92+
.filter(pl.col("sum_intensity_in_photons") > PHOTON_THRESHOLD)
93+
.collect()
94+
)
95+
print(f" {len(spot_data):,} spots above threshold")
96+
97+
result = af.number_of_puncta_per_segmented_cell_with_threshold(
98+
analysis_file=spot_csv,
99+
analysis_data_raw=spot_data,
100+
threshold_lower=PHOTON_THRESHOLD,
101+
lower_cell_size_threshold=LOWER_CELL_SIZE,
102+
upper_cell_size_threshold=UPPER_CELL_SIZE,
103+
blur_degree=BLUR_DEGREE,
104+
cell_string=CELL_STRING,
105+
protein_string=PROTEIN_STRING,
106+
dims=3,
107+
spacing=SPACING,
108+
collect_oligomers=True,
109+
)
110+
111+
if result is None or not isinstance(result, tuple):
112+
print(f" WARNING no output for {region}/{cell_type}/{condition}")
113+
return
114+
115+
cell_df, oligo_df = result
116+
117+
if cell_df is not None and isinstance(cell_df, pl.DataFrame):
118+
cell_df.write_csv(cell_csv_out)
119+
print(f" Saved {len(cell_df):,} cells → {cell_csv_out}")
120+
else:
121+
print(f" WARNING no cell DataFrame returned")
122+
123+
if oligo_df is not None and isinstance(oligo_df, pl.DataFrame):
124+
oligo_df.write_csv(oligo_out)
125+
print(f" Saved {len(oligo_df):,} in-cell oligomers → {oligo_out}")
126+
else:
127+
print(f" WARNING no in-cell oligomers found (check cell masks exist)")
128+
129+
130+
def main():
131+
parser = argparse.ArgumentParser(
132+
description="Re-run cell analysis to generate cell_id_in_image and oligomer linkage files."
133+
)
134+
group = parser.add_mutually_exclusive_group(required=True)
135+
group.add_argument(
136+
"--test", action="store_true", help="Run on cingulate/microglia only"
137+
)
138+
group.add_argument(
139+
"--full", action="store_true", help="Run on all 38 analysis directories"
140+
)
141+
args = parser.parse_args()
142+
143+
af = Analysis_Functions()
144+
targets = [("cingulate", "microglia")] if args.test else ALL_TARGETS
145+
146+
for region, cell_type in targets:
147+
for condition in ("HC", "PD"):
148+
print(f"\n=== {region} / {cell_type} / {condition} ===")
149+
process(af, region, cell_type, condition)
150+
151+
print("\nDone.")
152+
153+
154+
if __name__ == "__main__":
155+
main()

‎src/AnalysisFunctions.py‎

Lines changed: 84 additions & 59 deletions
Original file line numberDiff line numberDiff line change
@@ -49,65 +49,31 @@ def count_spots(self, database, threshold=None):
4949
Returns:
5050
polars.DataFrame: DataFrame with number of spots per z-plane.
5151
"""
52-
data = None
53-
results = None
52+
if database is None or database.is_empty():
53+
return None
5454
if threshold is None:
55-
columns = ["n_spots", "z", "image_filename"]
55+
data = (
56+
database.group_by(["image_filename", "z"])
57+
.agg(pl.len().alias("n_spots"))
58+
.sort(["image_filename", "z"])
59+
.select(["n_spots", "z", "image_filename"])
60+
)
61+
return H_F.clean_database(data, ["n_spots", "z", "image_filename"])
5662
else:
57-
columns = ["n_spots_above", "n_spots_below", "z", "image_filename"]
58-
for image_filename in np.sort(np.unique(database["image_filename"])):
59-
dataslice = database.filter(pl.col("image_filename") == image_filename)
60-
z_planes = np.sort(np.unique(dataslice["z"].to_numpy()))
61-
if threshold is None:
62-
spots_per_plane = [
63-
len(dataslice.filter(pl.col("z") == z)["sum_intensity_in_photons"])
64-
for z in z_planes
65-
]
66-
else:
67-
spots_above = [
68-
np.sum(
69-
dataslice.filter(pl.col("z") == z)[
70-
"sum_intensity_in_photons"
71-
].to_numpy()
72-
> threshold
73-
)
74-
for z in z_planes
75-
]
76-
spots_below = [
77-
np.sum(
78-
dataslice.filter(pl.col("z") == z)[
79-
"sum_intensity_in_photons"
80-
].to_numpy()
81-
<= threshold
82-
)
83-
for z in z_planes
84-
]
85-
spots_per_plane = np.vstack([spots_above, spots_below])
86-
if results is not None:
87-
results = np.hstack(
88-
[
89-
results,
90-
np.vstack(
91-
[
92-
spots_per_plane,
93-
z_planes,
94-
np.full(len(z_planes), image_filename),
95-
]
96-
),
97-
]
98-
)
99-
else:
100-
results = np.vstack(
101-
[spots_per_plane, z_planes, np.full(len(z_planes), image_filename)]
102-
)
103-
if results is not None:
104-
spot_save = {}
105-
for i, column in enumerate(columns):
106-
spot_save[column] = results[i, :]
107-
data = H_F.clean_database(pl.DataFrame(data=spot_save), columns)
108-
else:
109-
data = None
110-
return data
63+
data = (
64+
database.group_by(["image_filename", "z"])
65+
.agg([
66+
(pl.col("sum_intensity_in_photons") > threshold)
67+
.sum().alias("n_spots_above"),
68+
(pl.col("sum_intensity_in_photons") <= threshold)
69+
.sum().alias("n_spots_below"),
70+
])
71+
.sort(["image_filename", "z"])
72+
.select(["n_spots_above", "n_spots_below", "z", "image_filename"])
73+
)
74+
return H_F.clean_database(
75+
data, ["n_spots_above", "n_spots_below", "z", "image_filename"]
76+
)
11177

11278
def generate_indices(self, data, image_size, is_mask=False, is_lo=False):
11379
"""
@@ -893,6 +859,7 @@ def number_of_puncta_per_segmented_cell_with_threshold(
893859
sigma1=2.0,
894860
sigma2=40.0,
895861
spacing=(0.5, 0.11, 0.11),
862+
collect_oligomers=False,
896863
):
897864
"""
898865
Does analysis of number of oligomers in a mask area per "segmented" cell area.
@@ -943,6 +910,7 @@ def number_of_puncta_per_segmented_cell_with_threshold(
943910

944911
files = np.unique(analysis_data["image_filename"].to_numpy())
945912
cell_punctum_analysis = None
913+
oligomers_df = None
946914
start = time.time()
947915
analysis_directory = os.path.split(analysis_file)[0]
948916

@@ -1080,6 +1048,19 @@ def number_of_puncta_per_segmented_cell_with_threshold(
10801048
)
10811049
filename_tosave = np.full_like(centroids[:, 0], file, dtype="object")
10821050

1051+
# Collect per-spot arrays (bounds-filtered) for oligomers-per-cell output
1052+
if collect_oligomers and dims == 3:
1053+
all_spot_flat = np.ravel_multi_index(
1054+
np.vstack([z.astype(int), x.astype(int), y.astype(int)]),
1055+
image_size,
1056+
order="F",
1057+
)
1058+
sum_photons_b = subset["sum_intensity_in_photons"].to_numpy()[bounds]
1059+
bg_punctum_b = subset["bg_per_punctum"].to_numpy()[bounds]
1060+
bg_pixel_b = subset["bg_per_pixel"].to_numpy()[bounds]
1061+
zi_b = subset["zi"].to_numpy()[bounds]
1062+
zf_b = subset["zf"].to_numpy()[bounds]
1063+
10831064
def areaanalysis(coords):
10841065
if coords.shape[-1] == 2:
10851066
xm, ym = coords[:, 0], coords[:, 1]
@@ -1115,6 +1096,47 @@ def areaanalysis(coords):
11151096
n_cell_ratios = np.array([r[0] for r in results])
11161097
n_spots_in_object = np.array([r[1] for r in results])
11171098

1099+
# Build cell→oligomer mapping via a label map over the image
1100+
if collect_oligomers and dims == 3 and len(pil_mask) > 0:
1101+
n_pixels_total = int(np.prod(image_size))
1102+
cell_label_map = np.zeros(n_pixels_total, dtype=np.int32)
1103+
for j_cell, coords in enumerate(pil_mask):
1104+
c = coords.astype(int)
1105+
if (
1106+
np.any(c[:, 0] >= z_lim)
1107+
or np.any(c[:, 1] >= x_lim)
1108+
or np.any(c[:, 2] >= y_lim)
1109+
):
1110+
continue
1111+
flat_px = np.ravel_multi_index(c.T, image_size, order="F")
1112+
cell_label_map[flat_px] = j_cell + 1 # 1-based label
1113+
1114+
spot_cell_ids = cell_label_map[all_spot_flat]
1115+
del cell_label_map
1116+
in_any_cell = spot_cell_ids > 0
1117+
1118+
if in_any_cell.any():
1119+
n_in = int(in_any_cell.sum())
1120+
img_oligo_df = pl.DataFrame(
1121+
{
1122+
"z": z[in_any_cell].astype("float64"),
1123+
"x": x[in_any_cell].astype("float64"),
1124+
"y": y[in_any_cell].astype("float64"),
1125+
"sum_intensity_in_photons": sum_photons_b[in_any_cell],
1126+
"bg_per_punctum": bg_punctum_b[in_any_cell],
1127+
"bg_per_pixel": bg_pixel_b[in_any_cell],
1128+
"zi": zi_b[in_any_cell],
1129+
"zf": zf_b[in_any_cell],
1130+
"image_filename": [file] * n_in,
1131+
"cell_id_in_image": spot_cell_ids[in_any_cell].tolist(),
1132+
}
1133+
)
1134+
oligomers_df = (
1135+
pl.concat([oligomers_df, img_oligo_df])
1136+
if oligomers_df is not None
1137+
else img_oligo_df
1138+
)
1139+
11181140
if len(areas) > 0:
11191141
if dims == 2:
11201142
data = {
@@ -1127,6 +1149,7 @@ def areaanalysis(coords):
11271149
else "cell_protein_load"
11281150
): n_cell_ratios,
11291151
"n_puncta_in_cell": n_spots_in_object,
1152+
"cell_id_in_image": np.arange(1, len(areas) + 1, dtype=np.int32),
11301153
"image_filename": filename_tosave,
11311154
}
11321155
else:
@@ -1141,6 +1164,7 @@ def areaanalysis(coords):
11411164
else "cell_protein_load"
11421165
): n_cell_ratios,
11431166
"n_puncta_in_cell": n_spots_in_object,
1167+
"cell_id_in_image": np.arange(1, len(areas) + 1, dtype=np.int32),
11441168
"image_filename": filename_tosave,
11451169
}
11461170

@@ -1167,8 +1191,9 @@ def areaanalysis(coords):
11671191

11681192
if cell_punctum_analysis is None:
11691193
print("WARNING: No cells with valid areas found in any files.")
1170-
return None
1171-
return cell_punctum_analysis.rechunk()
1194+
return (None, None) if collect_oligomers else None
1195+
result = cell_punctum_analysis.rechunk()
1196+
return (result, oligomers_df) if collect_oligomers else result
11721197

11731198
def colocalise_spots_with_threshold(
11741199
self,

‎src/HelperFunctions.py‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -126,7 +126,7 @@ def make_datarray_largeobjects(
126126
else:
127127
dataarray = np.squeeze(stack.T)
128128
if dataarray is not None:
129-
dataarray = np.asarray(np.matrix(np.squeeze(dataarray)).transpose())
129+
dataarray = np.squeeze(dataarray).T
130130
if len(dataarray.shape) > 1:
131131
return pl.DataFrame(
132132
data=dataarray,
@@ -201,7 +201,7 @@ def make_datarray_spot(
201201
else:
202202
dataarray = da
203203
if dataarray is not None:
204-
dataarray = np.asarray(np.matrix(np.squeeze(dataarray)).transpose())
204+
dataarray = np.squeeze(dataarray).T
205205
if len(dataarray.shape) > 1:
206206
return pl.DataFrame(
207207
data=dataarray,

0 commit comments

Comments
 (0)