Proj mat calc fix

morph_utils/ccf.py

@@ -102,6 +102,30 @@ def load_structure_graph():
     df = df.set_index('acronym')
     return df
 
+def de_layer(st):
+    """de-layer cortical projection targets
+
+    Args:
+        st (str): e.g. ipsi_VISal2/3
+
+    Returns:
+        str: e.g. ipsi_VISal
+    """
+    CTX_STRUCTS = STRUCTURE_DESCENDANTS_ACRONYM['CTX']
+    sub_st = st.replace("ipsi_","").replace("contra_","")
+    if sub_st in CTX_STRUCTS:
+
+        for l in ["1","2/3","4","5","6a","6b"]:
+            st = st.replace(l,"")
+
+        if "ENT" in st:
+            for l in ["2", "3", "5/6", "6"]:
+                st = st.replace(l,"")
+
+        return st
+    else:
+        return st
+
 
 def process_pin_jblob( slide_specimen_id, jblob, annotation, structures, prints=False) :
     """
@@ -303,6 +327,7 @@ def get_ccf_structure(voxel, name_map=None, annotation=None, coordinate_to_voxel
     return name_map[structure_id]
 
 def projection_matrix_for_swc(input_swc_file, mask_method = "tip_and_branch", 
+                              apply_mask_at_cortical_parent_level=False,
                               count_method = "node", annotation=None, 
                               annotation_path = None, volume_shape=(1320, 800, 1140),
                               resolution=10, node_type_list=[2],
@@ -316,6 +341,9 @@ def projection_matrix_for_swc(input_swc_file, mask_method = "tip_and_branch",
         'tip_and_branch' will return a projection matrix masking only structures with tip and branch nodes. If 'tip'
         will only look at structures with tip nodes. And last, if 'branch' will only look at structures with 
         branch nodes.
+        apply_mask_at_cortical_parent_level (bool): If True, the `mask_method` will be applied to aggregated cortical
+        regions. E.g. if `mask_method`='tip_and_branch' and apply_mask_at_cortical_parent_level = True, then 
+        the tip-and-branch mask will be enforced at the (e.g.) VISp level, instead of in VISp1, VISp2/3 etc. independantly
         count_method (str): ['node','tip','branch']. When 'node', will measure axon length directly.
         Otherwise will return the count of tip or branch nodes in each structure
         annotation (array, optional): 3 dimensional ccf annotation array. Defaults to None.
@@ -402,11 +430,21 @@ def node_ider(morph,i):
 
     # determine ipsi/contra projections
     morph_df["ccf_structure_sided"] = morph_df.apply(lambda row: "ipsi_{}".format(row.ccf_structure) if row.z<z_midline else "contra_{}".format(row.ccf_structure), axis=1)
-    # mask the morphology dataframe accordinagly
-    if mask_method!=None:
-
+    if apply_mask_at_cortical_parent_level:   
+        morph_df['ccf_structure_rollup'] =  morph_df['ccf_structure'].map(de_layer) 
+        morph_df["ccf_structure_sided_rollup"] = morph_df.apply(lambda row: "ipsi_{}".format(row.ccf_structure_rollup) if row.z<z_midline else "contra_{}".format(row.ccf_structure_rollup), axis=1)
+
+    # mask the morphology dataframe accordingly 
+    if mask_method is not None:
         keep_structs = []
-        for struct, struct_df in morph_df.groupby("ccf_structure_sided"):
+        for struct in morph_df['ccf_structure_sided'].unique():
+
+            if apply_mask_at_cortical_parent_level:
+                sided_parent_struct = de_layer(struct)
+                struct_df = morph_df[morph_df['ccf_structure_sided_rollup']==sided_parent_struct]
+            else:
+                struct_df = morph_df[morph_df['ccf_structure_sided']==struct]
+
             node_types_in_struct = struct_df.node_type.unique().tolist()
             if (mask_method == 'tip') and ("tip" in node_types_in_struct):
                 keep_structs.append(struct)

morph_utils/executable_scripts/distance_between_nodes_for_directory.py

@@ -10,7 +10,7 @@
 class IO_Schema(ags.ArgSchema):
     swc_input_directory = ags.fields.InputDir(description='directory with swc files')
     output_file = ags.fields.OutputFile(descripion='output csv with distances between files')
-    compartment_types = ags.fields.List(default=[3, 4], cls_or_instance=ags.fields.Int)
+    compartment_types = ags.fields.List(default=[2, 3, 4], cls_or_instance=ags.fields.Int)
     use_multiprocessing = ags.fields.Boolean(default=True)
 
 
@@ -29,6 +29,7 @@ def main(swc_input_directory, output_file, compartment_types, use_multiprocessin
     # all_combinations = [c for c in all_combinations] # if c[0] != c[1]]
 
     print("{} Comparisons to analyze".format(len(all_combinations)))
+    print(f"Comparing nodes of type: {compartment_types}")
     reslist = []
     for combo in all_combinations:
         file_1 = os.path.join(swc_input_directory, combo[0])

morph_utils/executable_scripts/projection_matrix_for_single_cell.py

@@ -7,16 +7,17 @@
 class IO_Schema(ags.ArgSchema):
     input_swc_file = ags.fields.InputFile(description='directory with micron resolution ccf registered files')
     output_projection_csv = ags.fields.OutputFile(description="output projection csv")
+    mask_method = ags.fields.Str(description = " 'tip_and_branch', 'branch', 'tip', or 'tip_or_branch' ")
+    apply_mask_at_cortical_parent_level = ags.fields.Bool( descriptions='If True, the `mask_method` will be applied at aggregated cortical regions')
+
     projection_threshold = ags.fields.Int(default=0)
     normalize_proj_mat = ags.fields.Boolean(default=True)
-    mask_method = ags.fields.Str(default="tip_and_branch",description = " 'tip_and_branch', 'branch', 'tip', or 'tip_or_branch' ")
     count_method = ags.fields.String(default="node", description="should be a member of ['node','tip','branch']")
     annotation_path = ags.fields.Str(default="",description = "Optional. Path to annotation .nrrd file. Defaults to 10um ccf atlas")
     resolution = ags.fields.Int(default=10, description="Optional. ccf resolution (micron/pixel")
     volume_shape = ags.fields.List(ags.fields.Int, default=[1320, 800, 1140], description = "Optional. Size of input annotation")
     resample_spacing = ags.fields.Float(allow_none=True, default=None, description = 'internode spacing to resample input morphology with')
 
-
 def normalize_projection_columns_per_cell(input_df, projection_column_identifiers=['ipsi', 'contra']):
     """
     :param input_df:  input projection df
@@ -42,12 +43,15 @@ def main(input_swc_file,
          annotation_path,
          volume_shape,
          resample_spacing,
+         apply_mask_at_cortical_parent_level,
          **kwargs):
 
     if annotation_path == "":
         annotation_path = None
-
-    if mask_method not in [None,'tip_and_branch', 'branch', 'tip', 'tip_or_branch']:
+
+    if mask_method is None:
+        mask_method = "None"
+    if mask_method not in [None, 'None', 'tip_and_branch', 'branch', 'tip', 'tip_or_branch']:
         raise ValueError(f"Invalid mask_method provided {mask_method}")  
 
     results = []
@@ -58,7 +62,8 @@ def main(input_swc_file,
                                     annotation_path = annotation_path, 
                                     volume_shape=volume_shape,
                                     resolution=resolution,
-                                    resample_spacing=resample_spacing)
+                                    resample_spacing=resample_spacing,
+                                    apply_mask_at_cortical_parent_level=apply_mask_at_cortical_parent_level)
     results = [res]
 
     output_projection_csv = output_projection_csv.replace(".csv", f"_{mask_method}.csv")

morph_utils/executable_scripts/projection_matrix_from_swc_directory.py

@@ -5,16 +5,17 @@
 import time 
 import subprocess
 from morph_utils.ccf import projection_matrix_for_swc
-from morph_utils.proj_mat_utils import roll_up_proj_mat
+from morph_utils.proj_mat_utils import roll_up_proj_mat, normalize_projection_columns_per_cell
 
 
 class IO_Schema(ags.ArgSchema):
     ccf_swc_directory = ags.fields.InputDir(description='directory with micron resolution ccf registered files')
     output_directory = ags.fields.OutputDir(description="output directory")
+    mask_method = ags.fields.Str(description = " 'tip_and_branch', 'branch', 'tip', or 'tip_or_branch' ")
+    apply_mask_at_cortical_parent_level = ags.fields.Bool( descriptions='If True, the `mask_method` will be applied at aggregated cortical regions')
+    count_method = ags.fields.String(default="node", description="should be a member of ['node','tip','branch']")
     projection_threshold = ags.fields.Int(default=0)
     normalize_proj_mat = ags.fields.Boolean(default=True)
-    mask_method = ags.fields.Str(default="tip_and_branch",description = " 'tip_and_branch', 'branch', 'tip', or 'tip_or_branch' ")
-    count_method = ags.fields.String(default="node", description="should be a member of ['node','tip','branch']")
     annotation_path = ags.fields.Str(default="",description = "Optional. Path to annotation .nrrd file. Defaults to 10um ccf atlas")
     resolution = ags.fields.Int(default=10, description="Optional. ccf resolution (micron/pixel")
     volume_shape = ags.fields.List(ags.fields.Int, default=[1320, 800, 1140], description = "Optional. Size of input annotation")
@@ -23,20 +24,6 @@ class IO_Schema(ags.ArgSchema):
     virtual_env_name = ags.fields.Str(default='skeleton_keys_4',description='Name of virtual conda env to activate on hpc. not needed if running local')
     output_projection_csv = ags.fields.OutputFile(description="output projection csv, when running local only")
 
-
-def normalize_projection_columns_per_cell(input_df, projection_column_identifiers=['ipsi', 'contra']):
-    """
-    :param input_df:  input projection df
-    :param projection_column_identifiers: list of identifiers for projection columns. i.e. strings that identify projection columns from metadata columns
-    :return: normalized projection matrix
-    """
-    proj_cols = [c for c in input_df.columns if any([ider in c for ider in projection_column_identifiers])]
-    input_df[proj_cols] = input_df[proj_cols].fillna(0)
-
-    res = input_df[proj_cols].T / input_df[proj_cols].sum(axis=1)
-    input_df[proj_cols] = res.T
-
-    return input_df
 
 
 def main(ccf_swc_directory, 
@@ -45,6 +32,7 @@ def main(ccf_swc_directory,
          projection_threshold, 
          normalize_proj_mat,
          mask_method,
+         apply_mask_at_cortical_parent_level,
          count_method,
          annotation_path,
          volume_shape,
@@ -56,7 +44,10 @@ def main(ccf_swc_directory,
 
     if run_host not in ['local','hpc']:
         raise ValueError(f"Invalid run_host parameter entered ({run_host})")
-    if mask_method not in [None,'tip_and_branch', 'branch', 'tip', 'tip_or_branch']:
+
+    if mask_method is None:
+        mask_method = "None"
+    if mask_method not in ["None",'tip_and_branch', 'branch', 'tip', 'tip_or_branch']:
         raise ValueError(f"Invalid mask_method provided {mask_method}")
 
     if annotation_path == "":
@@ -85,7 +76,9 @@ def main(ccf_swc_directory,
                                 annotation_path = annotation_path, 
                                 volume_shape=volume_shape,
                                 resolution=resolution,
-                                resample_spacing= resample_spacing)
+                                resample_spacing= resample_spacing,
+                                apply_mask_at_cortical_parent_level = apply_mask_at_cortical_parent_level,
+                                )
             results.append(res)
 
         else:
@@ -106,7 +99,7 @@ def main(ccf_swc_directory,
                 command = command+ f" --count_method {count_method}"
                 command = command+ f" --annotation_path {annotation_path}"
                 command = command+ f" --resolution {resolution}"
-                # command = command+ f" --volume_shape {volume_shape}"
+                command = command+ f" --apply_mask_at_cortical_parent_level {apply_mask_at_cortical_parent_level}"
                 if resample_spacing is not None:
                     command = command+ f" --resample_spacing {resample_spacing}"
 

morph_utils/proj_mat_utils.py

@@ -0,0 +1,65 @@
+import os
+import numpy as np
+import pandas as pd
+
+from morph_utils.ccf import de_layer 
+
+
+def roll_up_proj_mat(infile, outfile):
+
+    df = pd.read_csv(infile, index_col=0)
+    df.index = df.index.map(os.path.basename)
+
+    non_proj_cols = [f for f in df.columns if not any([i in f for i in ["ipsi","contra"]])]
+    new_df = df[non_proj_cols].copy()
+
+    proj_cols = [f for f in df.columns if any([i in f for i in ["ipsi","contra"]])]
+    de_layer_dict = {p:de_layer(p) for p in proj_cols}
+
+    parent_names = list(de_layer_dict.values())
+    unique_parent_names = np.unique(parent_names)
+    unique_parent_names = sorted(unique_parent_names, key=lambda x:parent_names.index(x))
+
+    roll_up_records = {}
+    for low_res_struct in unique_parent_names:
+        children = [k for k,v in de_layer_dict.items() if v==low_res_struct ]
+        roll_up_records[low_res_struct] = children
+
+
+
+    # for parent, child_list in roll_up_records.items():
+    #     new_df[parent] = df[child_list].sum(axis=1)
+    new_cols = {
+        parent: df[child_list].sum(axis=1)
+        for parent, child_list in roll_up_records.items()
+    }
+    new_cols_df = pd.DataFrame(new_cols)
+    new_df = pd.concat([new_df, new_cols_df], axis=1)
+
+    # sanity check
+    for n_struct,old_list in roll_up_records.items():
+        sum_old = df[old_list].sum(axis=1)
+        sum_new = new_df[n_struct]
+        assert sum(sum_old==sum_new) == len(df)
+
+
+
+    # print(outfile)
+    # print()
+    assert os.path.abspath(outfile) != os.path.abspath(infile)
+    new_df.to_csv(outfile)
+
+
+def normalize_projection_columns_per_cell(input_df, projection_column_identifiers=['ipsi', 'contra']):
+    """
+    :param input_df:  input projection df
+    :param projection_column_identifiers: list of identifiers for projection columns. i.e. strings that identify projection columns from metadata columns
+    :return: normalized projection matrix
+    """
+    proj_cols = [c for c in input_df.columns if any([ider in c for ider in projection_column_identifiers])]
+    input_df[proj_cols] = input_df[proj_cols].fillna(0)
+
+    res = input_df[proj_cols].T / input_df[proj_cols].sum(axis=1)
+    input_df[proj_cols] = res.T
+
+    return input_df

setup.cfg

@@ -29,3 +29,4 @@ console_scripts =
     morph_utils_aggregate_single_cell_projs = morph_utils.executable_scripts.aggregate_single_cell_projection_csvs:console_script
     morph_utils_move_somas_left_hemisphere = morph_utils.executable_scripts.move_somas_to_left_hemisphere_swc_directory:console_script
     morph_utils_local_crop_ccf_swcs = morph_utils.executable_scripts.local_crop_ccf_swc_directory:console_script
+    morph_utils_dsit_btwn_nodes_directory = morph_utils.executable_scripts.distance_between_nodes_for_directory:console_script