Read SWC files in nanometer units

[dp4846]

Thanks for the great package! I wanted to give it a try on a neuron morphology but when I tried to load an swc file from flywire using :

import jaxley as jx
nm = '720575940590515268.swc'
cell = jx.read_swc(nm, nseg=1)

I get the error:

{
	"name": "IndexError",
	"message": "arrays used as indices must be of integer (or boolean) type",
	"stack": "---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
File /Users/dean/Desktop/code/science/modules/equi_dyn_imp/test_jaxley.py:4
      2 import jaxley as jx
      3 nm = '/Users/dean/Desktop/code/science/modules/equi_dyn_imp/sk_lod1_783_healed/720575940623940963.swc'
----> 4 cell = jx.read_swc(nm, nseg=1)

File ~/opt/anaconda3/envs/jaxley_env/lib/python3.10/site-packages/jaxley/modules/cell.py:390, in read_swc(fname, nseg, max_branch_len, min_radius, assign_groups)
    362 def read_swc(
    363     fname: str,
    364     nseg: int,
   (...)
    367     assign_groups: bool = False,
    368 ) -> Cell:
    369     \"\"\"Reads SWC file into a `jx.Cell`.
    370 
    371     Jaxley assumes cylindrical compartments and therefore defines length and radius
   (...)
    388         A `jx.Cell` object.
    389     \"\"\"
--> 390     parents, pathlengths, radius_fns, types, coords_of_branches = swc_to_jaxley(
    391         fname, max_branch_len=max_branch_len, sort=True, num_lines=None
    392     )
    393     nbranches = len(parents)
    395     comp = Compartment()

File ~/opt/anaconda3/envs/jaxley_env/lib/python3.10/site-packages/jaxley/utils/swc.py:38, in swc_to_jaxley(fname, max_branch_len, sort, num_lines)
     31 if is_single_point_soma:
     32     # Warn here, but the conversion of the length happens in `_compute_pathlengths`.
     33     warn(
     34         \"Found a soma which consists of a single traced point. `Jaxley` \"
     35         \"interprets this soma as a spherical compartment with radius \"
     36         \"specified in the SWC file, i.e. with surface area 4*pi*r*r.\"
     37     )
---> 38 sorted_branches, types = _split_into_branches_and_sort(
     39     content,
     40     max_branch_len=max_branch_len,
     41     is_single_point_soma=is_single_point_soma,
     42     sort=sort,
     43 )
     45 parents = _build_parents(sorted_branches)
     46 each_length = _compute_pathlengths(
     47     sorted_branches, content[:, 1:6], is_single_point_soma=is_single_point_soma
     48 )

File ~/opt/anaconda3/envs/jaxley_env/lib/python3.10/site-packages/jaxley/utils/swc.py:93, in _split_into_branches_and_sort(content, max_branch_len, is_single_point_soma, sort)
     86 def _split_into_branches_and_sort(
     87     content: np.ndarray,
     88     max_branch_len: float,
     89     is_single_point_soma: bool,
     90     sort: bool = True,
     91 ) -> Tuple[np.ndarray, np.ndarray]:
     92     branches, types = _split_into_branches(content, is_single_point_soma)
---> 93     branches, types = _split_long_branches(
     94         branches,
     95         types,
     96         content,
     97         max_branch_len,
     98         is_single_point_soma=is_single_point_soma,
     99     )
    101     if sort:
    102         first_val = np.asarray([b[0] for b in branches])

File ~/opt/anaconda3/envs/jaxley_env/lib/python3.10/site-packages/jaxley/utils/swc.py:131, in _split_long_branches(branches, types, content, max_branch_len, is_single_point_soma)
    129 num_subbranches += 1
    130 split_branch = _split_branch_equally(branch, num_subbranches)
--> 131 lengths_of_subbranches = _compute_pathlengths(
    132     split_branch,
    133     coords=content[:, 1:6],
    134     is_single_point_soma=is_single_point_soma,
    135 )
    136 lengths_of_subbranches = [
    137     np.sum(length_traced) for length_traced in lengths_of_subbranches
    138 ]
    139 length = max(lengths_of_subbranches)

File ~/opt/anaconda3/envs/jaxley_env/lib/python3.10/site-packages/jaxley/utils/swc.py:296, in _compute_pathlengths(all_branches, coords, is_single_point_soma)
    294 branch_pathlengths = []
    295 for b in all_branches:
--> 296     coords_in_branch = coords[np.asarray(b) - 1]
    297     if len(coords_in_branch) > 1:
    298         # If the branch starts at a different neurite (e.g. the soma) then NEURON
    299         # ignores the distance from that initial point. To reproduce, use the
   (...)
    302         # 2 2 9.00 0.0 0.0 0.5 1
    303         # 3 2 10.0 0.0 0.0 0.3 2
    304         types = coords_in_branch[:, 0]

IndexError: arrays used as indices must be of integer (or boolean) type"
}

The swc file:

# SWC format file
# based on specifications at http://www.neuronland.org/NLMorphologyConverter/MorphologyFormats/SWC/Spec.html
# Created on 2023-11-11 using navis (https://github.com/navis-org/navis)
# Meta: {"id": "720575940590515268", "name": "None", "units": "1 nanometer"}
# PointNo Label X Y Z Radius Parent
# Labels:
# 0 = undefined, 1 = soma, 5 = fork point, 6 = end point
1 1 191532.69 421374.75 115987.305 2238 -1
2 0 192377.12 420970.44 116211.14 2386 10
3 0 191525.58 421437.7 116000.59 2186 1
4 0 191571.89 421589.03 116106.99 2025 3
5 0 192243.7 421278.66 116512.14 111 6
6 0 192182.8 421356.8 116571.05 94 7
7 0 192063.5 421239.16 116727.695 126 34
8 0 192477.55 420327.28 116631.4 433 13
9 0 192269.16 421333.12 116449.88 142 5
10 0 192187.52 421176.38 116098.0 2307 4
11 5 192535.8 420809.06 116381.16 2455 2
12 0 192172.28 421048.7 116722.36 96 8
13 0 192213.83 420853.34 116777.36 119 33
14 0 192273.95 421553.12 116403.25 130 9
15 0 192263.62 421607.34 116336.4 113 14
16 6 192223.39 421769.06 116293.71 60 15
17 0 192501.23 420332.78 117299.76 625 23
18 0 192476.67 420347.47 117279.695 530 17
19 0 192393.4 420064.84 117279.54 1064 21
20 0 192370.98 420051.38 117301.19 942 19
21 0 192314.83 420101.5 117318.914 1215 22
22 0 192309.8 420100.78 117293.79 1320 29
23 0 192287.52 420052.16 117304.1 764 30
24 0 192547.8 420394.03 116812.12 1860 40
25 0 192482.69 420684.25 116369.88 2388 11
26 0 195497.5 421955.12 117498.63 143 11
27 0 195659.4 421068.56 116916.2 120 11
28 0 192350.75 420181.84 117196.43 1570 32
29 0 192320.97 420082.34 117256.39 1462 28
30 0 192286.86 420073.97 117325.92 865 20
31 0 192473.73 420301.16 116946.234 1833 24
32 0 192387.45 420209.53 117169.72 1686 39
33 0 192180.38 420748.78 116888.04 72 36
34 0 192081.55 421135.2 116738.04 124 12
35 0 192650.03 420560.47 116547.17 2279 38
36 0 192238.52 419984.88 117029.68 415 18
37 0 192373.19 420423.0 116943.73 2128 41
38 0 192661.17 420631.06 116422.64 2332 25
39 0 192381.53 420234.7 117085.51 1792 31
40 0 192431.44 420277.97 117007.8 2058 37
41 0 192576.64 420380.72 116928.86 2156 42
42 0 192506.31 420279.9 116684.8 2100 43
43 0 192603.6 420411.75 116670.77 2177 35
44 6 195495.88 421991.2 117509.625 91 26
45 0 195751.19 421059.1 116869.59 122 27
46 6 195878.67 421126.38 116857.23 93 45

Does the swc need to be reformatted?

Best I could tell the error was happening when in `_compute_pathlengths' :

    for b in all_branches:
        coords_in_branch = coords[np.asarray(b) - 1]

b was an empty list.

[jnsbck]

Hey, thanks for reaching out.
From your .swc file it looks like you are using $nm$ as the unit of length. This trips up the max_branch_len default, which is 300.0($\mu m$). Jaxley expects units to be given in $\mu m$, similar to NEURON, see here. I suggest changing the units in your swc file.
I hope this fixes your problem. Please lemme know if this is not the case, then I will look into it further.

Best,

Jonas

[michaeldeistler]

I think we could also directly parse this in Jaxley. We could just extract the meta:

# Meta: {"id": "720575940590515268", "name": "None", "units": "1 nanometer"}

And then use it to directly read nanometer units.