Chartrand et al. (2023): morphology claims from the paper's feature table

Paper. Chartrand et al. (2023), Science 382, eadf0805, doi:10.1126/science.adf0805; see reproduce_chartrand_2023.ipynb for the text read and the table of all claims.

Data. These claims rest on the reconstructed morphologies, which are not in DANDI:000630 (the dandiset holds the electrophysiology only). The paper deposited the reconstructions in the Brain Image Library (doi:10.35077/g.606) and the morphology features computed from them in its code archive (Zenodo doi:10.5281/zenodo.8310226; read here from github.com/AllenInstitute/patchseq_human_L1 at commit 7777857c1b89b19735e45fd9c9bd7493cc1d805c, file data/RawFeatureWide_human+derivatives.csv). The claims are data: elsewhere and are counted apart from the claims on the dandiset. The analysis uses the paper's per-cell morphology features as deposited; it does not recompute them from the reconstructions. The dandiset is read only to list which cells it holds.

Claim	Status
fig2d-subclass-morphology	tested
fig4b-mc4r-axon	tested
fig5b-tnfaip8l3-dendrites	tested
figs6a-tumor-dendrite-branching	tested, then blocked: the medical condition is missing for most cells

import numpy as np
import pandas as pd

import l1ephys as L

pd.set_option("display.width", 200)
labels = L.load_labels()
morph = L.read_table("RawFeatureWide_human+derivatives.csv", index_col="specimen_id").drop(columns=["Unnamed: 0"])
morph["axon_emd_with_basal_dendrite"] *= 5   # l1_load.py: units to microns
MORPH = list(morph.columns)
h = labels[labels.in_strict_set].join(morph)
h["has_morph_table"] = h.index.isin(morph.index)
h["homology_type"] = pd.Categorical(h["t-type"].map(L.HOMOLOGY), categories=L.SUBCLASSES + ["other"])
h["l1_type"] = h.homology_type.notna()
h["core_l1_type"] = h.homology_type.isin(L.SUBCLASSES)
dand = L.list_assets()
dand_ids = pd.to_numeric(dand.path.str.extract(r"ses-(\d+)")[0])
h["in_dandiset"] = h.ephys_roi_id.isin(dand_ids)
print(f"{len(morph)} cells in the morphology table, {len(MORPH)} features; strict L1 set {len(h)} cells, "
      f"{int(h.has_morph_table.sum())} with morphology")
print("Cells with morphology by site:", h[h.has_morph_table].collaborator.value_counts().to_dict(),
      "| of which in dandiset 000630:", int((h.has_morph_table & h.in_dandiset).sum()))
print("By subclass:", h[h.has_morph_table].homology_type.value_counts(sort=False).to_dict())
print("Feature ranges (5 examples):")
print(h.loc[h.has_morph_table, ["axon_extent_x", "axon_max_path_distance", "basal_dendrite_num_branches",
                                "basal_dendrite_extent_x", "basal_dendrite_soma_percentile_y"]].describe().round(2).to_string())

/tmp/ipykernel_3127/497383901.py:12: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  h["has_morph_table"] = h.index.isin(morph.index)

86 cells in the morphology table, 51 features; strict L1 set 272 cells, 78 with morphology
Cells with morphology by site: {'AIBS': 53, 'Gabor': 16, 'Huib': 9} | of which in dandiset 000630: 50
By subclass: {'LAMP5': 29, 'MC4R': 12, 'PAX6': 19, 'L1 VIP': 4, 'other': 7}
Feature ranges (5 examples):
       axon_extent_x  axon_max_path_distance  basal_dendrite_num_branches  basal_dendrite_extent_x  basal_dendrite_soma_percentile_y
count          78.00                   78.00                        78.00                    78.00                             78.00
mean          623.23                 1006.66                        32.68                   183.23                              0.57
std           361.66                  413.95                        13.97                   136.58                              0.26
min           168.93                   26.52                         9.00                    32.66                              0.01
25%           415.98                  753.43                        25.00                   102.43                              0.40
50%           611.78                  967.99                        32.00                   140.94                              0.58
75%           767.62                 1210.98                        38.00                   192.04                              0.80
max          2752.96                 2257.62                        89.00                   818.76                              0.99

/tmp/ipykernel_3127/497383901.py:18: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  h["in_dandiset"] = h.ephys_roi_id.isin(dand_ids)

The morphology table holds the reconstructed cells of all three sites, linked to the cell table by specimen id. Of the 78 cells with morphology in the strict L1 set, 53 are from Seattle and 50 of these are cells of the dandiset, but the morphologies themselves are not in it.

fig2d-subclass-morphology

Human L1 subclasses differ morphologically in vertical orientation of axons and dendrites, axon extent and shape, and dendrite branching (Kruskal-Wallis, FDR p < 1e-3).

As l1_subclass_distinctness.ipynb cells 19 to 21: strict L1 cells with morphology in the four core subclasses; Kruskal-Wallis per morphology feature, Benjamini-Hochberg; the features Fig. 2D plots (cell 21) are the fraction of dendrite below the soma, the axon horizontal extent and the number of dendrite branches. The claims file gives no criterion; the claim is judged on these three having FDR p < 1e-3, with the other features named in the statement reported.

sub = h[h.has_morph_table & h.core_l1_type].copy()
sub["homology_type"] = sub.homology_type.cat.remove_unused_categories()
res, pairs = L.run_kw_dunn(sub, MORPH, "homology_type", posthoc="dunn")
print(f"{len(sub)} cells", sub.homology_type.value_counts(sort=False).to_dict(), f"| {len(res)} features, "
      f"{int((res.pval_fdr < 1e-3).sum())} with FDR p < 1e-3")
fig2d = ["basal_dendrite_soma_percentile_y", "axon_extent_x", "basal_dendrite_num_branches"]
print(res.loc[fig2d, ["pval", "pval_fdr", "epsilon2", "n"]].to_string())
others = [f for f in ["axon_soma_percentile_y", "axon_bias_y", "axon_extent_y", "axon_extent_x_over_y",
                      "basal_dendrite_bias_y", "axon_max_path_distance"] if f in res.index]
print("Other orientation, extent and shape features:"); print(res.loc[others, ["pval_fdr", "epsilon2"]].to_string())
print("Medians:"); print(sub.groupby("homology_type", observed=True)[fig2d].median().round(2).to_string())

64 cells {'LAMP5': 29, 'MC4R': 12, 'PAX6': 19, 'L1 VIP': 4} | 51 features, 4 with FDR p < 1e-3
                                      pval  pval_fdr  epsilon2   n
feature                                                           
basal_dendrite_soma_percentile_y  0.000008  0.000203  0.418712  64
axon_extent_x                     0.000024  0.000405  0.382557  64
basal_dendrite_num_branches       0.000605  0.003430  0.275033  64
Other orientation, extent and shape features:
                        pval_fdr  epsilon2
feature                                   
axon_soma_percentile_y  0.008276  0.236277
axon_bias_y             0.091546  0.121413
axon_extent_y           0.047302  0.155453
axon_extent_x_over_y    0.008276  0.235750
basal_dendrite_bias_y   0.002760  0.290706
axon_max_path_distance  0.003197  0.281326
Medians:
               basal_dendrite_soma_percentile_y  axon_extent_x  basal_dendrite_num_branches
homology_type                                                                              
LAMP5                                      0.73         587.00                         40.0
MC4R                                       0.20         260.26                         30.5
PAX6                                       0.71         683.42                         27.0
L1 VIP                                     0.38         632.40                         24.5

Result: partially reproduced. Across the four core subclasses (64 reconstructed cells) the fraction of dendrite below the soma (FDR p = 2.0e-4) and the axon horizontal extent (4.1e-4) differ below 1e-3. The number of dendrite branches differs with FDR p = 0.0034, below 0.05 but not below the claim's 1e-3, and the axon's vertical orientation features have FDR p 0.008 to 0.09.

fig4b-mc4r-axon

MC4R cells are quantitatively distinct from other L1 types in maximum axonal path distance and axon branch frequency.

As l1_rosehip.ipynb cells 4, 11 and 17: strict L1 cells, MC4R against the rest; branch frequency is the number of axon branches per 100 um of axon; exact Mann-Whitney over all morphology features and the branch frequency; Benjamini-Hochberg. Judged on FDR p < 1e-4 (Fig. 4B legend) for both.

d = h[h.has_morph_table].copy()
d["axon_branch_freq"] = d["axon_num_branches"] / d["axon_total_length"] * 100
d["in_type"] = np.where(d.homology_type == "MC4R", "MC4R", "rest")
res = L.run_twosamp(d, MORPH + ["axon_branch_freq"], "in_type")
med = d.groupby("in_type")[["axon_max_path_distance", "axon_branch_freq"]].median()
print(f"MC4R {int((d.in_type == 'MC4R').sum())} cells against {int((d.in_type == 'rest').sum())}")
print(res.loc[["axon_max_path_distance", "axon_branch_freq"], ["pval_mw", "pval_mw_fdr", "mw_r", "nobs"]].to_string())
print("Medians:"); print(med.round(2).to_string())

MC4R 12 cells against 66
                             pval_mw  pval_mw_fdr      mw_r  nobs
feature                                                          
axon_max_path_distance  3.313772e-04     0.002872  0.626263    78
axon_branch_freq        2.839605e-08     0.000001  0.876263    78
Medians:
         axon_max_path_distance  axon_branch_freq
in_type                                          
MC4R                     662.30              3.01
rest                     991.93              1.36

Result: partially reproduced. MC4R cells (12 against 66) have a higher axon branch frequency (3.01 against 1.36 per 100 um, FDR p = 1.2e-6, below 1e-4) and a shorter maximum axonal path (662 against 992 um), but the path difference has FDR p = 0.0029, not below 1e-4.

fig5b-tnfaip8l3-dendrites

PAX6 TNFAIP8L3 cells have more dendritic branches and a horizontal dendritic extent over 550 um wide.

As l1_tnfaip.ipynb cells 16 and 17: strict L1 cells, PAX6 TNFAIP8L3 against the rest; exact Mann-Whitney over the morphology features; Benjamini-Hochberg; judged on FDR p < 0.05 (Fig. 5B legend) with TNFAIP8L3 higher for both features, and each TNFAIP8L3 cell's horizontal extent above 550 um.

d = h[h.has_morph_table].copy()
d["in_type"] = np.where(d["t-type"] == "PAX6 TNFAIP8L3", "TNFAIP8L3", "rest")
res = L.run_twosamp(d, MORPH, "in_type")
f2 = ["basal_dendrite_num_branches", "basal_dendrite_extent_x"]
print(f"TNFAIP8L3 {int((d.in_type == 'TNFAIP8L3').sum())} cells against {int((d.in_type == 'rest').sum())}")
print(res.loc[f2, ["pval_mw", "pval_mw_fdr", "mw_r", "nobs"]].to_string())
print("TNFAIP8L3 cells:"); print(d.loc[d.in_type == "TNFAIP8L3", f2 + ["collaborator"]].round(1).to_string())
print("Rest, medians:", d.loc[d.in_type == "rest", f2].median().round(1).to_dict(), "| rest maximum extent:",
      round(d.loc[d.in_type == "rest", "basal_dendrite_extent_x"].max(), 1))

TNFAIP8L3 2 cells against 76
                              pval_mw  pval_mw_fdr      mw_r  nobs
feature                                                           
basal_dendrite_num_branches  0.013320     0.084915  0.907895    78
basal_dendrite_extent_x      0.000666     0.016983  1.000000    78
TNFAIP8L3 cells:
           basal_dendrite_num_branches  basal_dendrite_extent_x collaborator
spec_id                                                                     
840142358                         50.0                    583.2         AIBS
828758585                         68.0                    818.8        Gabor

Rest, medians: {'basal_dendrite_num_branches': 32.0, 'basal_dendrite_extent_x': 139.9} | rest maximum extent: 571.0

Result: partially reproduced. The two PAX6 TNFAIP8L3 cells have horizontal dendritic extents of 583 and 819 um, both over 550 um, against a median of 140 um in the other 76 cells (FDR p = 0.017). They have more dendrite branches (50 and 68 against a median of 32), but with FDR p = 0.085, not below 0.05.

figs6a-tumor-dendrite-branching

Cells from tumor patients have higher dendritic branching than cells from epilepsy patients, not explained by brain area or subclass.

As l1_metadata_effects.ipynb cells 65, 66, 72 and 73: strict L1 cells with an L1 t-type and morphology, epileptic against non-epileptic donors (code's epileptic(), on the public table's medical_conditions in place of the code's LIMS records); exact Mann-Whitney over the morphology features; Benjamini-Hochberg. Judged on the number of dendrite branches having FDR p < 0.05 and being higher in non-epileptic (tumor) donors, overall and within MTG LAMP5 cells (cell 72, the check that area and subclass do not explain it).

def epileptic(x):
    if pd.isna(x):
        return x
    return "epileptic" if "epilep" in x.lower() else "non-epileptic"

d = h[h.has_morph_table & h.l1_type].copy()
d["epileptic"] = d.medical_conditions.map(epileptic)
print("Medical condition in the public table, cells with morphology:", d.medical_conditions.value_counts(dropna=False).to_dict())
for name, data in [("all", d), ("MTG LAMP5", d[(d.structure == "MTG") & (d.homology_type == "LAMP5")])]:
    data = data.dropna(subset=["epileptic"])
    if data.epileptic.nunique() < 2:
        print(f"{name}: only one group with a known condition ({data.epileptic.value_counts().to_dict()}); not testable")
        continue
    res = L.run_twosamp(data, MORPH, "epileptic")
    med = data.groupby("epileptic").basal_dendrite_num_branches.median()
    print(f"{name}: {data.epileptic.value_counts().to_dict()}; dendrite branches median {med.round(1).to_dict()}, "
          f"p {res.loc['basal_dendrite_num_branches', 'pval_mw']:.3g}, FDR p {res.loc['basal_dendrite_num_branches', 'pval_mw_fdr']:.3g}; "
          f"smallest FDR p {res.pval_mw_fdr.min():.3g} ({res.pval_mw_fdr.idxmin()})")

Medical condition in the public table, cells with morphology: {nan: 51, 'Epilepsy': 10, 'Tumor status - primary': 4, 'Tumor': 3, 'Epileptiform Events': 3}
all: {'epileptic': 13, 'non-epileptic': 7}; dendrite branches median {'epileptic': 26.0, 'non-epileptic': 37.0}, p 0.351, FDR p 0.942; smallest FDR p 0.568 (axon_exit_theta)

MTG LAMP5: only one group with a known condition ({'epileptic': 7}); not testable

Result: blocked, metadata insufficient. The public cell table gives a medical condition for only 20 of the 71 reconstructed cells with an L1 t-type (13 epileptic, 7 not; 51 have none), where the paper's code read the conditions for all donors from the Allen Institute's donor database. On these 20 cells the number of dendrite branches is higher in non-epileptic donors (median 37 against 26) without significance (FDR p = 0.94), and the check within MTG LAMP5 cells cannot be run (only epileptic donors are labelled). The test cannot be repeated on the paper's cells.

Summary

Claim	Outcome
fig2d-subclass-morphology	partially reproduced
fig4b-mc4r-axon	partially reproduced
fig5b-tnfaip8l3-dendrites	partially reproduced
figs6a-tumor-dendrite-branching	blocked, metadata insufficient

A failure to reproduce here may be a failure of this attempt, not of the paper: see the reasoning and the differences from the paper's method.

The analysis uses the paper's deposited morphology features, not the reconstructions. The strict L1 set holds 78 reconstructed cells in the feature table, against the 71 of Table S2; which cells the paper counted is not settled here.