Specific expansion of motor cortical projections in a singing mouse.

Code ↔ Paper

7 matches between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 7 matches

[1] § Methods › MAPseq data analysis › Heatmaps and cell-type labelling ↔ MAPseq/ED_fig5_UMI_thresholds_sexes.ipynb, lines 76–132 · score 0.64 · CT cells, PT cells, HB, striatum, thalamus, HY
[2] § Methods › MAPseq data analysis › Heatmaps and cell-type labelling ↔ MAPseq/ED_fig4_infectivity_downsampling.ipynb, lines 100–156 · score 0.63 · CT cells, PT cells, striatum, thalamus, infectivity, HY
[3] § Methods › MAPseq data analysis › Degree and motif analyses ↔ MAPseq/fig4_motifs.ipynb, lines 354–406 · score 0.56 · unobserved neurons, Motif proportions, population, probabilities, cells, brain
[4] § Methods › MAPseq data analysis › Degree and motif analyses ↔ MAPseq/ED_fig4_infectivity_downsampling.ipynb, lines 513–565 · score 0.56 · unobserved neurons, Motif proportions, population, probabilities, cells, brain
[5] § Testing models of brain-wide cortical projections ↔ MAPseq/ED_fig4_infectivity_downsampling.ipynb, lines 100–156 · score 0.54 · pyramidal tract, intratelencephalic, striatum, thalamus, cortical, PT
[6] § Testing models of brain-wide cortical projections ↔ MAPseq/ED_fig5_UMI_thresholds_sexes.ipynb, lines 76–132 · score 0.54 · pyramidal tract, intratelencephalic, striatum, thalamus, cortical, PT
[7] § Methods › MAPseq tissue processing ↔ MAPseq/ED_fig4_infectivity_downsampling.ipynb, lines 59–96 · score 0.53 · projecting neurons, Injection sites, Allen, infections, dissection

Paper

Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The paper is loaded by your browser, with JavaScript: read it at doi.org or on Europe PMC.

The authors' code

Jupyter notebook · 998 lines · 35 KB · no license · 4 matches

# %%
###### load packages
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import math
from itertools import combinations, chain
from scipy import stats
from upsetplot import from_memberships
from scipy.stats import chi2, binomtest
from statsmodels.distributions.empirical_distribution import ECDF # for generating empirical cdfs
import matplotlib.lines as mlines # needed for custom legend
import os
# needed for editable text in svg format
plt.rcParams['font.family'] = ['sans-serif']
plt.rcParams['font.sans-serif'] = ['Arial']
plt.rcParams['text.usetex'] = False
plt.rcParams['svg.fonttype'] = 'none'
# load metadata
metadata = pd.read_csv("metadata.csv")
# import custom colormaps
from colormaps import blue_cmp, orange_cmp
# specify blue/orange colors for qualitative intervals/data
blue_qual = [blue_cmp.colors[50], blue_cmp.colors[100], blue_cmp.colors[150], blue_cmp.colors[200], blue_cmp.colors[250]]
orange_qual = [orange_cmp.colors[36], orange_cmp.colors[72], orange_cmp.colors[108], orange_cmp.colors[144], orange_cmp.colors[180], orange_cmp.colors[216], orange_cmp.colors[252]]
%matplotlib inline
# %% [markdown]
# # Load data
# %%
# set-up out paths
out_path = '/your/path/here/'
# %%
# load data
in_path = 'your/path/here'
# load binarized data
omc_bin = []
for idx, row in metadata.iterrows():
df = pd.read_csv(os.path.join(in_path, row["mice"]+"_"+row["dataset"]+"_bin.csv"), index_col=0)
omc_bin.append(df)
omc_countN = []
for idx, row in metadata.iterrows():
df = pd.read_csv(os.path.join(in_path, row["mice"]+"_"+row["dataset"]+"_countN.csv"), index_col=0)
omc_countN.append(df)
# %% [markdown]
# # Functions
# %% [markdown]
# ## Processing functions
# %%
def clean_up_data(df_dirty, to_drop = ['OB', 'ACAi', 'ACAc', 'HIP'], inj_site="OMCi"):
"""Clean up datasets so all matrices are in the same format. Function
(1) drops unwanted columns, e.g. negative controls or dissections of other injection sites.
(2) renames RN (allen acronym) as BS (my brainstem acronym)
(3) drops neurons w/ 0 projections after dropped columns
Args:
df_dirty (DataFrame): Pandas dataframe (neurons x area) that needs to be processed
to_drop (list, optional): columns to drop. Defaults to ['OB', 'ACAi', 'ACAc', 'HIP'].
inj_site (str, optional): Injection site. Defaults to "OMCi".
Returns:
DataFrame: Cleaned up data in dataframe format.
"""
# 1. drop unused areas
dropped = df_dirty.drop(to_drop, axis=1)
# 2. change RN to bs
replaced = dropped.rename(columns={'RN':'BS'})
# 3. drop neurons w/ 0 projections after removing negative regions
if type(inj_site)==str:
nodes = replaced.drop([inj_site], axis=1).sum(axis=1)
else:
nodes = replaced.drop(inj_site, axis=1).sum(axis=1)
n_idx = nodes > 0 # non-zero projecting neurons
clean = replaced[n_idx]
return clean
# %%
def sort_by_celltype(proj, it_areas=["OMCc", "AUD", "STR"], ct_areas=["TH"],
pt_areas=["AMY","HY","SNr","SCm","PG","PAG","BS"],
sort=True):
"""
Function takes in projection matrix and outputs matrix sorted by the 3 major celltypes:
- IT = intratelencephalic (projects to cortical and/or Striatum), type = 10
- CT = corticalthalamic (projects to thalamus w/o projection to brainstem), type = 100
- PT = pyramidal tract (projects to brainstem += other areas), type = 1000
Returns single dataframe with cells sorted and labelled by 3 cell types (IT/CT/PT)
Args:
proj (DataFrame): pd.DataFrame of BC x area. Entries can be normalized BC or binary.
it_areas (list, optional): Areas to determine IT cells. Defaults to ["OMCc", "AUD", "STR"].
ct_areas (list, optional): Areas to determine CT cells. Defaults to ["TH"]. Don't actually use this...
pt_areas (list, optional): Areas to determine PT cells. Defaults to ["AMY","HY","SNr","SCm","PG","PAG","BS"].
sort (bool, optional): Whether to sort by cell type or return w/ original index. Defaults to True.
Returns:
df_out (DataFrame): Returns dataframe with extra column (type) labelled and sorted by cell type
"""
ds=proj.copy()
# Isolate PT cells
pt_counts = ds[pt_areas].sum(axis=1)
pt_idx = ds[pt_counts>0].index
ds_pt = ds.loc[pt_idx,:]
ds_pt['type'] = "PT"
# Isolate remaining non-PT cells
ds_npt = ds.drop(pt_idx)
# Identify CT cells by thalamus projection
th_idx = ds_npt['TH'] > 0
ds_th = ds_npt[th_idx]
if sort:
ds_th = ds_th.sort_values('TH', ascending=False)
ds_th['type'] = "CT"
# Identify IT cells by the remaining cells (non-PT, non-CT)
ds_nth = ds_npt[~th_idx]
if sort:
ds_nth = ds_nth.sort_values(it_areas,ascending=False)
ds_nth['type'] = "IT"
# combine IT and CT cells
ds_npt = pd.concat([ds_nth, ds_th])
# combine IT/CT and PT cells
if sort:
sorted = pd.concat([ds_npt,ds_pt],ignore_index=True)
df_out=sorted.reset_index(drop=True)
else:
df_out = pd.concat([ds_npt,ds_pt]).sort_index()
return(df_out)
# %%
def dfs_preprocess_counts(df_list, drop=["type"],
norm_by="inj_median", inj_site="OMCi", medians=None):
"""Take dataframe and process it for downstream analysis. Can be binarized or
non-binarized data.
Args:
df_list (list): List of dataframes,
drop (list, optional): List of columns to drop before returning. Defaults to ["OMCi", "type"].
norm_by (str, optional): How to normalize counts. Can be "inj_median" or 'all_median' or "given". Defaults to "inj_median".
inj_site (str, optional): What column to use for injection norm. Defaults to "OMCi".
medians (list, optional): if norm_by=="given", uses list of medians (same length as df_list) to normalize dataframes. Defaults to None.
Returns:
out_list (list): List of dataframes same order as input list.
"""
out_list = []
norms_out = []
for i in range(len(df_list)):
df = df_list[i].drop(drop, axis=1)
if norm_by=="inj_row":
# normalize on neuron by neuron basis to injection site value
out_df = df.div(df[inj_site], axis=0)
elif norm_by=="given":
out_df = df/medians[i]
elif norm_by=="inj_mean":
vals = df[inj_site].values.flatten()
mean = np.mean(vals)
norms_out.append(mean)
out_df = df/mean
else:
# normalize by non-zero median of injection site
if norm_by=="inj_median":
vals = df[inj_site].values.flatten() # added injection site median
# normalize by median of all non-zero values median across whole dataset
elif norm_by=="all_median":
vals = df.values.flatten()
idx = vals.nonzero() # only use non-zero ncounts for determining median
plot = vals[idx]
median = np.median(plot)
norms_out.append(median)
out_df = df/median
for j in range(len(drop)):
out_df[drop[j]] = df_list[i][drop[j]] # add dropped columns back in, can preserve labels
out_list.append(out_df)
if len(norms_out)>0:
return(norms_out, out_list)
else:
return(out_list)
# %%
def dfs_to_cdf(df_list, plot_areas, resolution=1000, metadata=metadata):
"""Takes in list of DFs of count(N) data and returns dataframe w/ cdf data that can be plotted.
Returned Dataframe includes metadata
Args:
df_list (list): list of DataFrames of count(N) data.
plot_areas (list): List of strings of areas to calculate cdfs.
resolution (int, optional): Used to determine resolution of cdf line. Defaults to 1000.
medatadata (df, optional): Metadata where row corresponds to df_list indices. Defaults to metadata.
Returns:
cdf_df (DataFrame): Dataframe of cdf summary values
all_ecdfs (list): Empirical cdf
"""
# combine all DFs into one df labelled w/ metadata
all_bc = pd.DataFrame(columns=list(df_list[0].columns)+["mice", "species", "dataset"])
for i in range(metadata.shape[0]):
df = df_list[i].copy(deep=True)
df['mice'] = metadata.loc[i, 'mice']
df['species'] = metadata.loc[i, "species"]
df['dataset'] = metadata.loc[i, "dataset"]
all_bc = pd.concat([all_bc, df])
all_bc = all_bc.reset_index(drop=True)
cdf_df = pd.DataFrame(columns=["x", "cdf", "mice", "species", "dataset", "area"])
all_ecdfs = {}
# calculate cdf by area, then add by mouse
for area in plot_areas:
# just use nonzero BC
area_idx = all_bc[area] > 0
area_bc = all_bc.loc[area_idx, [area, "mice", "species", "dataset"]]
# get min/max for each area to set cdf bounds
area_min = area_bc[area].min()
area_max = area_bc[area].max()
for i in range(metadata.shape[0]):
micei = metadata.loc[i, 'mice']
mice_bc = area_bc[area_bc['mice']==micei]
if mice_bc[area].sum()==0:
print("NO BARCODES, cannot compute ECDF for", area, metadata.loc[i,'mice'])
else:
# print(area, metadata.loc[i,"mice"])
ecdf = ECDF(mice_bc[area])
x = np.logspace(np.log10(area_min), np.log10(area_max), num=resolution)
y = ecdf(x)
int = pd.DataFrame({"x":x, "cdf":y, "mice":metadata.loc[i,"mice"], "species":metadata.loc[i,"species"],
"dataset":metadata.loc[i,"dataset"], "area":area})
cdf_df = pd.concat([cdf_df, int])
all_ecdfs[micei+"_"+area] = ecdf
return(cdf_df, all_ecdfs)
# %%
def dfs_to_proportions(df_list, drop=["OMCi", "type"], keep=None, cell_type=None, meta=metadata,
inj_site="OMCi", aud_rename=False):
"""Output dataframe of proportions in format that can be plotted with seaborn
Args:
df_list (list):
- List of dataframes of neurons/BC by areas
drop (list, optional):
- Defaults to ["OMCi", "type"]
- list of areas/columns to drop before calculating proportions
cell_type (string, optional):
- Specify cell types in df, either IT, CT or PT
- Defaults to None
Returns:
plot_df (pandas_dataframe):
- returns dataframe in format for seaborn plotting
- columns = areas, and other metadata
"""
plot_df = pd.DataFrame(columns=["area", "proportion", "mice", "species", "dataset"])
if cell_type == "IT":
drop.extend([inj_site, 'TH', 'HY', 'AMY', 'SNr', 'SCm', 'PG',
'PAG', 'BS'])
elif cell_type == "PT":
if aud_rename:
drop.extend([inj_site,inj_site[:-1]+"c", 'AUD/TEa', "STR"])
else:
drop.extend([inj_site,inj_site[:-1]+"c", 'AUD', "STR"])
if keep:
drop = []
mice = meta["mice"]
species = meta["species"]
dataset = meta["dataset"]
for i in range(len(df_list)):
df = df_list[i].drop(drop, axis=1)
if keep:
df = df.loc[:, keep] # just subset keep columns
bc_sum = df.sum()
proportion = bc_sum/df.shape[0]
df_add = pd.DataFrame({"area":proportion.index.values, "proportion":proportion.values,
"mice":meta.loc[i,"mice"], "species":meta.loc[i,"species"], "dataset":meta.loc[i,"dataset"]})
# add cell_type if specified
if cell_type:
df_add["type"] = cell_type
plot_df = pd.concat([plot_df, df_add])
return plot_df
# %%
def downsample_neurons(data, meta=metadata, random_state=10, species="MMus", sample_ns=None):
"""Given list of dataframe, sample from combined neurons/cells (with replacement), in equivalent numbers
to singing mouse (or sample_ns) cells/brain. Returns list where each element is dataframe of neurons w/ numbers equivalent to ns.
Args:
data (list): List of pandas dataframes where each row is a different cell/neuron
metadata (DataFrame, optional): Metadata used to determine which df in data is lab/singing mouse.
Defaults to metadata.
random_state (int, optional): Set random state to use for repeatable sampling.
Defaults to 10.
species (str, optional): Species to sample. Defaults to "MMus".
sample_ns (list, optional): List of int use as sample size per 'brain', if none defaults to singing mouse brain size.
Defaults to None.
"""
all = [data[i] for i in range(metadata.shape[0]) if metadata.loc[i,'species']==species]
all = pd.concat(all).reset_index(drop=True)
# print("mm_all.shape[0]", mm_all.shape[0])
pool = all.copy()
samp = []
ns = []
if sample_ns:
ns = sample_ns
else:
# get steg sample neuron sizes
for i in range(metadata.shape[0]):
if metadata.loc[i,'species']=="STeg":
ns.append(data[i].shape[0])
for i in range(len(ns)):
n = ns[i]
int = pool.sample(n, random_state=random_state+i) # can't have same random_state for every round or will sample the same neurons
samp.append(int.reset_index(drop=True))
return(samp)
# %%
def statistic_testing(df, sp1="MMus", sp2="STeg", to_plot='proportion',
groupby="area", test="ttest"):
"""output dataframe based on comparison of species individuals
output dataframe can be used for making volcano plot
Args:
df (DataFrame): DataFrame of repeated values across species and inviduals.
sp1 (str, optional): Group1 to use as comparison. Defaults to "MMus".
sp2 (str, optional): Group1 to use as comparison. Defaults to "STeg".
to_plot (str, optional): Column to calculate comparison. Defaults to 'proportion'.
groupby (str, optional): Categories within to do comparison. Defaults to "area".
test (str, optional): What test to use ["ttest", "mannwhitneyu"]. Defaults to "ttest".
Returns:
plot (DataFrame): pd.DataFrame w/ pvalues, means across species, and log. Can be used for volcano plot.
"""
groups = sorted(df[groupby].unique())
# sp1
sp1_df = df[df["species"]==sp1]
sp1_array = sp1_df.pivot_table(columns='mice', values=to_plot, index=groupby).values
sp2_df = df[df["species"]==sp2]
sp2_array = sp2_df.pivot_table(columns='mice', values=to_plot, index=groupby).values
if test=="ttest":
results = stats.ttest_ind(sp1_array, sp2_array, axis=1)
elif test=="mannwhitneyu":
results = stats.mannwhitneyu(sp1_array, sp2_array, axis=1)
p_vals = results[1]
plot = pd.DataFrame({groupby:groups, "p-value":p_vals})
plot[sp1+"_mean"] = sp1_array.mean(axis=1)
plot[sp1+"_sem"] = stats.sem(sp1_array, axis=1)
plot[sp2+"_mean"] = sp2_array.mean(axis=1)
plot[sp2+"_sem"] = stats.sem(sp2_array, axis=1)
# plot["effect_size"] = (plot["st_mean"]-plot["mm_mean"]) / (plot["st_mean"] + plot["mm_mean"]) # modulation index
plot["fold_change"] = plot[sp2+"_mean"]/(plot[sp1+"_mean"])
plot["log2_fc"] = np.log2(plot["fold_change"])
plot["nlog10_p"] = -np.log10(plot["p-value"])
plot["p<0.05"] = plot["p-value"]<0.05
return(plot)
# %%
def estimate_n_total(df, plot_areas):
"""Calculate estimated N_total based on observed number of neuron per areas.
Works with IT and PT areas with as many areas as wanted. See Han et al., 2017
for more detail/formula.
Args:
df (DataFrame): Binary matrix of BC x areas
areas (list): List of areas using to calculat motifs. Should be list of strings.
Returns:
n_total (int): Estimated n total for input data.
"""
n_obs = df.shape[0]
n_areas = df.sum()
all_terms = []
for k in range(1, len(plot_areas)+1):
combos = list(combinations(plot_areas, k))
term = 0
for i in range(len(combos)):
product = 1
for j in range(len(combos[i])):
n_area = n_areas[combos[i][j]]
product = product*n_area
term = term + product
all_terms.append(term)
# need to subtract first term from n_obs
all_terms[0] = n_obs - all_terms[0]
# multiply every other by -1 starting w/ 3rd term
for l in range(len(all_terms)):
if l>=2 and l%2==0:
all_terms[l] = -1*all_terms[l]
# find roots of polynomial
roots = np.roots(all_terms)
# convert to real numbers if numbers complex
if isinstance(roots[0],complex):
reals = []
for num in roots:
if num.imag==0:
reals.append(num.real)
else:
reals = roots.copy()
# pick root that is more than n_obs
# if can't find root more than n_obs, return n_obs as n_total
n_total = n_obs
for num in reals:
if num > n_obs:
n_total = round(num)
return(n_total)
# %%
def df_to_motif_proportion(df, areas, proportion=True, subset=None):
"""Make series to feed into upset plot based on given data and area list
Args:
df (pd.DataFrame): df containing bc x area
areas (list): List of areas (str) to be plotted
proportion (bool, optional): Whether to output counts or proportions. Defaults to True.
subset (str, optional): Whether to subset on certain brain area (e.g. 'PAG'). Defaults to None.
Returns:
plot_s (DataFrame): dataframe of proprotion/count of each motif for input data.
"""
# generate all combinations of areas in true/false list
area_comb = []
for i in range(len(areas)):
n = i+1
area_comb.append(list(combinations(areas, n)))
area_comb_list = list(chain.from_iterable(area_comb)) # flatten list
memberships = from_memberships(area_comb_list) # generate true/false
area_comb_TF = memberships.index.values # extract array w/ true/false values
area_comb_names = memberships.index.names # get order of areas
# calculate number of neurons of each motif
comb_count = []
for tf in area_comb_TF:
neurons = df
for i in range(len(area_comb_names)):
# subset dataset on presence/absence of proj
neurons = neurons[neurons[area_comb_names[i]]==tf[i]]
if proportion: # return proportion
comb_count.append(neurons.shape[0]/df.shape[0])
else: # return count
comb_count.append(neurons.shape[0])
plot_s = from_memberships(area_comb_list, data=comb_count)
if subset:
motif_areas = plot_s.index.names
subset_idx = motif_areas.index(subset)
idx = [i for i, x in enumerate(plot_s.index) if x[subset_idx]]
plot_s = plot_s[idx]
return(plot_s)
# %%
def df_to_motif_estimated_proportions(data, combinations, adjust_total=False):
"""Given dataframe of cells and index of combinations (generated from df_to_motif_prportions),
return series similar to output for df_to_motif_proportions. Proportions are estimated
by multiplying the bulk/population proporitons/probabilities.
Args:
data (DataFrame): DataFrame of BC x areas, binary
combinations (MultiIndex): Index from output of df-to_motif_proportions
adjust_total (bool): whether to adjust the matrix w/ unobserved neurons. Defaults to False.
Returns:
"""
# adjust n_total (add in 0 projectors), if not done
if adjust_total:
# calculate n_total
n_obs = data.shape[0]
n_total = estimate_n_total(data, combinations.names)
n_unobs = np.array(n_total )- np.array(n_obs)
unobs_df = pd.DataFrame(0, index=np.arange(n_unobs), columns=data.columns)
df = pd.concat([data, unobs_df]).reset_index(drop=True)
else:
df = data.copy()
# get brain areas specified in combinations
names = list(combinations.names)
# subset df to just columns used in motifs
df_subset = df.loc[:,names]
# get bulk proportions across those areas
bulk_prop = df_subset.sum(axis=0)/df.shape[0]
# calculate expected motif proportion based on product of bulk proportions
expected_proportions = []
names = combinations.names
for i in range(combinations.shape[0]):
motif = combinations[i]
product = 1
for j in range(len(names)):
if motif[j]: # if area in motif, multiply by bulk probability
product = product * bulk_prop[names[j]]
else: # if area not in motif, multiply by 1-bulk probability (chance of not projecting to area)
product = product * (1-bulk_prop[names[j]])
expected_proportions.append(product)
motif_expected_prop = pd.Series(expected_proportions, index=combinations)
return(motif_expected_prop)
# %%
def TF_to_motifs(index):
"""Generate motif strings from True/False MultiIndex genreated by df_to_motif_proportion()
Args:
index (MultiIndex): True/False for area combinations
Returns:
array of strings for area combinations
"""
motifs_strings = []
for r in index:
motif = ""
for i in range(len(index.names)):
if r[i]:
motif = motif+index.names[i]+"_"
motifs_strings.append(motif)
return(motifs_strings)
# %% [markdown]
# ## Plotting functions
# %%
def dot_plot(data, subset=None, title=None, err="se", add_legend=False,
to_plot="proportion", ylim=(0), fig_size=(3.5,3.5),
jitter=False, alpha=0.3):
"""Plot closed circle of value per area for data.
Args:
data (DataFrame): Dataframe of proportions, included downsampled data.
subset (list, optional): List of what to subset on to plot. Defaults to None.
title (str, optional): Title for plot. Defaults to None.
err (str, optional): Error to plot, can be "ci", "pi", "se", or "sd". Defaults to "se".
add_legend (bool, optional): Whether to add legend labeling mean/err. Defaults to False.
to_plot (str, optional): Column to plot. Defaults to "proportion".
ylim (tuple, optional): lower bound for yaxis. Defaults to (0).
fig_size (tuple, optional): Dimensions (in inches) of plot. Defaults to (3.5,3.5).
jitter (bool, optional): Whether to have overlapping points. Defaults to False.
alpha (float, optional): How much transparency for points. Defaults to 0.3.
"""
if subset:
df = data[data[subset[0]]==subset[1]]
else:
df = data.copy()
fig, ax = plt.subplots()
sns.stripplot(x='species', y=to_plot, hue="species",
dodge=False, data=df, jitter=jitter,
size=10, alpha=alpha)
sns.pointplot(df, x="species", y=to_plot, hue="species",
dodge=False, marker="+",
ls="", zorder=10, errorbar="se")
ax.set_xlabel("")
plt.title(title, size=20)
plt.ylim(ylim) # make sure y axis starts at 0
ax.spines['right'].set_visible(False)
ax.spines['top'].set_visible(False)
# set figure size
fig = plt.gcf()
fig.set_size_inches(fig_size[0],fig_size[1])
if add_legend:
legend = mlines.Line2D([], [], color="black", marker="[]", linewidth=0, label="mean, "+err)
plt.legend(handles=[legend], loc="lower right")
else:
plt.legend([],[], frameon=False)
return(fig)
# %%
def plot_cdf(data, plot_areas, log=True, title="", color_by="species", colors=[blue_cmp.colors[255], orange_cmp.colors[255]],
individual=True, meta=metadata, legend=True, fig_size=(3,3), calc_cdf=True):
"""Takes in countN data and returns cdf plots
Args:
data (list): list of DataFrames of count(N), where each element in animal
plot_areas (list): list of strings of areas to include in final output
log (bool, optional): Whether to use log on axis scale. Defaults to True.
title (str, optional): figure title. Defaults to "".
color_by (str, optional): Can be "mice", "species", or "dataset", what to label as metadata. Defaults to "species".
colors (list, optional): colors used to label cdfs. Defaults to [blue_cmp.colors[255], orange_cmp.colors[255]].
individual (bool, optional): _description_. Defaults to True.
meta (_type_, optional): _description_. Defaults to metadata.
calc_cdf (bool, optional): Whether to calcualte cdf or not. Defaults to True.
"""
# calculate ecdf per animal and put into dataframe
if calc_cdf:
cdf_df, foo = dfs_to_cdf(data, plot_areas=plot_areas, metadata=meta)
else:
cdf_df = data.copy()
# calculate number of axes needed
n = math.ceil(len(plot_areas)/5) # round up divide by 4 = axs rows
if len(plot_areas)==1:
fig, ax = plt.subplots(1,1, figsize=(5,5))
ax_list = [ax]
else:
fig, axs = plt.subplots(n, 5, figsize=(20, 5*n))
ax_list = axs.flat
i = 0
for ax in ax_list:
if i < len(plot_areas):
area = plot_areas[i]
plot = cdf_df[cdf_df['area']==area]
groups = plot[color_by].unique()
plot_1 = plot[plot[color_by] ==groups[0]]
plot_2 = plot[plot[color_by] ==groups[1]]
if individual:
sns.lineplot(plot_1, x="x", y="cdf", estimator=None, units="mice", color=colors[0], ax=ax) # plots individual mice
sns.lineplot(plot_2, x="x", y="cdf", estimator=None, units="mice", color=colors[1], ax=ax) # plots individual mice
else:
sns.lineplot(plot_1, x="x", y="cdf", ax=ax) # plots mean ci95
sns.lineplot(plot_2, x="x", y="cdf", ax=ax) # plots mean ci95
# get rid of top and right axis
ax.spines['right'].set_visible(False)
ax.spines['top'].set_visible(False)
if log:
ax.set_xscale("log")
ax.set_xlabel("Normalized Counts")
ax.set_title(area)
i+=1
else:
ax.axis('off')
# create cutom legend
if legend:
colors = [colors[0], colors[1]]
lines = [mlines.Line2D([0], [0], color=c, linewidth=3) for c in colors]
labels = [groups[0], groups[1]]
fig.legend(lines,labels, bbox_to_anchor=(0.75, 0.935))
# increase text size
plt.rcParams.update({'font.size': 12})
if title!="":
plt.suptitle(title, y=0.93, size=20)
elif title=="":
plt.suptitle("By "+color_by, y=0.93, size=20)
# erase minor ticks
plt.minorticks_off()
# set figure size
fig = plt.gcf()
fig.set_size_inches(fig_size[0],fig_size[1])
return(fig)
# %% [markdown]
# # Data pre-processing
# %%
# initial processing
# binary processing
omc_clean = [clean_up_data(df) for df in omc_bin]
omc_type = [sort_by_celltype(df) for df in omc_clean]
# normalized count (countN) processing
omc_cleanN = [clean_up_data(df) for df in omc_countN]
omc_typeN = [sort_by_celltype(df) for df in omc_cleanN]
means, omc_preprocessN = dfs_preprocess_counts(omc_typeN, norm_by="inj_mean") # normalize by injection mean
medians, omc_preprocessN_med = dfs_preprocess_counts(omc_typeN, norm_by="inj_median")
# %% [markdown]
# # Infectivity
# - Plots for Extended Data fig. 4
# %%
# Plot unique barcodes by species
infect_df = pd.DataFrame(columns=["Unique Barcodes", "mice", "species", "dataset"])
for i in range(metadata.shape[0]):
infect_df.loc[i,"Unique Barcodes"] = omc_type[i].shape[0]
infect_df.loc[i, "mice"] = metadata.loc[i,"mice"]
infect_df.loc[i,"species"] = metadata.loc[i, "species"]
infect_df.loc[i, "dataset"] = metadata.loc[i, "dataset"]
print("Mean")
display(infect_df.groupby("species")["Unique Barcodes"].mean())
print("SEM")
infect_df.groupby("species")["Unique Barcodes"].sem()
# %%
dot_plot(infect_df, to_plot="Unique Barcodes", fig_size=(3,3))
plt.yscale("log")
plt.ylim(100, 30000)
plt.minorticks_off()
# plt.savefig(out_path+"infectivity_uBC_dotplot.svg", dpi=300, bbox_inches="tight")
plt.show()
# %%
# plot barcode count per neuron for injection site (OMCi)
# NO WITHIN SAMPLE NORM - i.e. only normalized to spike-ins
plot_cdf(omc_typeN, plot_areas=["OMCi"], color_by="species",
title=None, legend=False, individual=False, fig_size=(2,2))
plt.title("OMCi")
# plt.savefig(out_path+"omci_rnaspike_norm_cdf.svg", dpi=300, bbox_inches="tight")
plt.show()
# %%
# take medians for each animal
omc_typeN[0]["OMCi"].median()
medians_omci = []
for i in range(metadata.shape[0]):
medians_omci.append(omc_typeN[i]["OMCi"].median())
omci_med = metadata.copy()
omci_med["inj_med"] = medians_omci
omci_med["area"] = "OMCi"
# calcualte mean, sem for each species
print("Mean")
print(omci_med.groupby("species")["inj_med"].mean())
print("\nSEM")
print(omci_med.groupby("species")["inj_med"].sem())
# Mann whitney u test
mm_med = omci_med[omci_med["species"]=="MMus"]
st_med = omci_med[omci_med["species"]=="STeg"]
print(stats.mannwhitneyu(mm_med["inj_med"], st_med["inj_med"]))
dot_plot(omci_med, to_plot="inj_med")
plt.show()
# %% [markdown]
# # Downsampling
# %%
# seperate it cells
omc_it = [df[df['type']=="IT"] for df in omc_type]
# %% [markdown]
# # Downsampling proportions
# %%
# seperate binarized counts into celltypes
it_bin = [df[df["type"]=="IT"].reset_index(drop=True) for df in omc_type]
pt_bin = [df[df["type"]=="PT"].reset_index(drop=True) for df in omc_type]
it_areas = ["OMCc", "AUD", "STR"]
pt_areas = ["TH", "HY", "AMY", "PAG", "SNr", "SCm", "PG", "BS"]
# calculate area proprotions
it_prop = dfs_to_proportions(it_bin, cell_type="IT", drop=["OMCi", "type"])
pt_prop = dfs_to_proportions(pt_bin, cell_type="PT", drop=["OMCi", "type"])
# %%
# downsample mmus neurons
mmus_down_meta = pd.DataFrame({"mice":["mm_down1", "mm_down2", "mm_down3", "mm_down4",
"mm_down5", "mm_down6", "mm_down7"],
"species":["MMus\ndownsample"]*7,
"dataset":["aggregate"]*7})
it_mmus_down = downsample_neurons(it_bin)
pt_mmus_down = downsample_neurons(pt_bin)
# calculate proportions
it_mm_down_prop = dfs_to_proportions(it_mmus_down, cell_type="IT", drop=["OMCi", "type"],
meta=mmus_down_meta)
pt_mm_down_prop = dfs_to_proportions(pt_mmus_down, cell_type="PT", drop=["OMCi", "type"],
meta=mmus_down_meta)
# add to proportion dataframes
it_prop_all = pd.concat([it_prop, it_mm_down_prop])
pt_prop_all = pd.concat([pt_prop, pt_mm_down_prop])
# %%
# plot downsampled dot plots
# plot it areas
for area in it_areas:
dot_plot(it_prop_all, subset=["area",area], title=area, add_legend=False,
jitter=True)
plt.ylabel("Propotion")
# plt.savefig(out_path+"OMC_"+area+"_proportion_dot_downsample.svg", dpi=300, bbox_inches="tight")
plt.show()
# plot pt areas
for area in pt_areas:
dot_plot(pt_prop_all, subset=["area",area], title=area, add_legend=False,
jitter=True)
plt.ylabel("Propotion")
# plt.savefig(out_path+"OMC_"+area+"_proportion_dot_downsample.svg", dpi=300, bbox_inches="tight")
plt.show()
# %%
# significance testing - IT
print("\nMANNWHITNEYU - mmus v. steg")
mannwhitney = statistic_testing(it_prop_all, test="mannwhitneyu", sp1="MMus", sp2="STeg")
display(mannwhitney)
print("\nMANNWHITNEYU - mmus downsample v. steg")
mannwhitney = statistic_testing(it_prop_all, test="mannwhitneyu", sp1="MMus\ndownsample", sp2="STeg")
display(mannwhitney)
# %%
# significance testing - PT
print("\nMANNWHITNEYU - mmus v. steg")
mannwhitney = statistic_testing(pt_prop_all, test="mannwhitneyu", sp1="MMus", sp2="STeg")
display(mannwhitney)
print("\nMANNWHITNEYU - mmus downsample v. steg")
mannwhitney = statistic_testing(pt_prop_all, test="mannwhitneyu", sp1="MMus\ndownsample", sp2="STeg")
display(mannwhitney)
# %% [markdown]
# # Downsampling motifs
# %%
# calculate estimates for motifs
plot_areas = ["OMCc", "AUD", "STR"]
# Estimate n-totals
n_totals = [estimate_n_total(omc_it[i], plot_areas) for i in range(len(omc_it))]
# Count obs motifs
n_obs_motifs = [df_to_motif_proportion(df, areas=plot_areas, proportion=False) for df in omc_it]
motifs = n_obs_motifs[0].index
# convert to proportions using adjusted n_totals
p_obs_motifs = [n_obs_motifs[i]/n_totals[i] for i in range(len(n_totals))]
# calculate expected proportions based on independent bulk probabilities adjusted for n_total
p_expected_motifs = [df_to_motif_estimated_proportions(df, motifs, adjust_total=True) for df in omc_it]
motif_strings = TF_to_motifs(motifs)
# put into dataframe - note p_obs is adjusted for n_total
it_motifs_df = []
for i in range(len(n_totals)):
df = pd.DataFrame({"motifs":motif_strings, "n_shape":omc_it[i].shape[0], "n_total":n_totals[i],
"n_obs":n_obs_motifs[i], "p_obs":p_obs_motifs[i], "p_exp":p_expected_motifs[i],
"mice":metadata.loc[i,"mice"], "species":metadata.loc[i,"species"]})
it_motifs_df.append(df)
it_motifs_df = pd.concat(it_motifs_df).reset_index(drop=True)
it_motifs_df
# %%
# downsample MMus and calculate motif proportions
# downsample
mmus_down_meta = pd.DataFrame({"mice":["mm_down1", "mm_down2", "mm_down3", "mm_down4",
"mm_down5", "mm_down6", "mm_down7"],
"species":["MMus\ndownsample"]*7,
"dataset":["aggregate"]*7})
it_mmus_down = downsample_neurons(omc_it)
# %%
# calculate estimates for motifs
plot_areas = ["OMCc", "AUD", "STR"]
# Estimate n-totals
n_totals = [estimate_n_total(it_mmus_down[i], plot_areas) for i in range(len(it_mmus_down))]
# Count obs motifs
n_obs_motifs = [df_to_motif_proportion(df, areas=plot_areas, proportion=False) for df in it_mmus_down]
motifs = n_obs_motifs[0].index
# convert to proportions using adjusted n_totals
p_obs_motifs = [n_obs_motifs[i]/n_totals[i] for i in range(len(n_totals))]
# calculate expected proportions based on independent bulk probabilities adjusted for n_total
p_expected_motifs = [df_to_motif_estimated_proportions(df, motifs, adjust_total=True) for df in it_mmus_down]
motif_strings = TF_to_motifs(motifs)
# put into dataframe - note p_obs is adjusted for n_total
it_motifs_df_mmd = []
for i in range(len(n_totals)):
df = pd.DataFrame({"motifs":motif_strings, "n_shape":it_mmus_down[i].shape[0], "n_total":n_totals[i],
"n_obs":n_obs_motifs[i], "p_obs":p_obs_motifs[i], "p_exp":p_expected_motifs[i],
"mice":mmus_down_meta.loc[i,"mice"], "species":mmus_down_meta.loc[i,"species"]})
it_motifs_df_mmd.append(df)
it_motifs_df_mmd = pd.concat(it_motifs_df_mmd).reset_index(drop=True)
# it_motifs_df_mmd
# %%
# plot
plot_df = pd.concat([it_motifs_df, it_motifs_df_mmd])
for motif in motif_strings:
dot_plot(plot_df, subset=["motifs",motif], to_plot="p_obs", title=motif, add_legend=False,
jitter=True)
plt.ylabel("Propotion")
# plt.savefig(out_path+"OMC_"+motif+"_proportion_dot_downsample.svg", dpi=300, bbox_inches="tight")
plt.show()
# %%
# Significance testing
print("\nMANNWHITNEYU - mmus v. steg")
mannwhitney = statistic_testing(plot_df, to_plot="p_obs", groupby="motifs", test="mannwhitneyu", sp1="MMus", sp2="STeg")
display(mannwhitney)
print("\nMANNWHITNEYU - mmus downsample v. steg")
mannwhitney = statistic_testing(plot_df, to_plot="p_obs", groupby="motifs", test="mannwhitneyu", sp1="MMus\ndownsample", sp2="STeg")
display(mannwhitney)

ED_fig4_infectivity_downsampling.ipynb at commit 9b5fa18, no license · at the source

Overview

Authors: Emily C. Isko^1,2, Clifford E. Harpole¹, Xiaoyue Mike Zheng^1,2, Huiqing Zhan¹, Martin B. Davis¹, Anthony M. Zador^1,2, Arkarup Banerjee^1,2

ORCID iDs: Emily C. Isko, Clifford E. Harpole, Xiaoyue Mike Zheng, Huiqing Zhan, Martin B. Davis, Anthony M. Zador, Arkarup Banerjee

Cold Spring Harbor Laboratory,Cold Spring Harbor, NY USA
School of Biological Sciences, Cold Spring Harbor Laboratory,Cold Spring Harbor, NY USA

Institutions: Cold Spring Harbor Laboratory (United States)

Journal: Nature, volume 655, issue 8122, pages 438-446

Dates: received 27 November 2024; accepted 26 March 2026; published online 6 May 2026; in print 2026

Type: Research article · Language: English

License: CC BY-NC-ND

Identifiers: DOI 10.1038/s41586-026-10458-y · PMID 42092144 · PMCID PMC13345907 · OpenAlex W7160437351

Open access: hybrid, a free copy (OpenAlex)

Status: code verified

Categories: mouse (organism), systems (subfield)

Methods: Connectivity, Statistics, Spectral & time-frequency, Evoked potentials, Graphs

Keywords: Neuroscience, Speciation

MeSH: Motor Cortex*, Rodentia*, Vocalization, Animal*, Animals, Auditory Cortex, Female, Male, Mice, Neural Pathways, Species Specificity (* major topic)

Topic: Animal Vocal Communication and Behavior (Developmental Biology, Biochemistry, Genetics and Molecular Biology), according to OpenAlex

Funding: NINDS NIH HHS (RF1 NS132046)

Citations: cited by 2 papers (Europe PMC); 81 references in the paper

Abstract

The abstract is not reproduced here: the paper's license (CC BY-NC-ND) does not allow it. Read it in the paper, at the publisher or on Europe PMC.

Repository

Its files are read in the Code ↔ Paper reader above, with 7 matches between paragraphs and lines of code.

singingmicelab/Isko-2026-MAPseq-code

License: none: the authors keep all their rights

State: the link answers, verified on 28 September 2026

Evidence: files inventoried

Commit: 9b5fa186b8cb4fee3e0f806829ff0eacb97c62ff, 5 February 2026

Languages: Jupyter (4), Python (1)

Size: 10 files, 5 scripts

Software Heritage: not archived

Found in: “Code availability”

Holds: README, 4 notebooks

Not found: license file, CITATION.cff, environment file, tests, continuous integration, documentation

Tools: Matplotlib (5 files), NumPy (5 files), pandas (4 files), SciPy (4 files), seaborn (4 files), statsmodels (3 files)

Availability: 1 check, the latest on 28 September 2026: the link answers

28 September 2026: the link answers

6 files

MAPseq/ED_fig4_infectivity_downsampling.ipynb, Jupyter, 998 lines, 4 matches
MAPseq/ED_fig5_UMI_thresholds_sexes.ipynb, Jupyter, 926 lines, 2 matches
MAPseq/colormaps.py, Python, 47 lines
MAPseq/fig3_MAPseq.ipynb, Jupyter, 1,016 lines
MAPseq/fig4_motifs.ipynb, Jupyter, 1,099 lines, 1 match
README.md, Text, 8 lines

At the source: github.com/singingmicelab/Isko-2026-MAPseq-code.git

Code availability statement

The paper has a code availability statement. Its license (CC BY-NC-ND) does not allow reproducing it here; in short, from what the harvester recognized in it:

it points to the authors' code: singingmicelab/Isko-2026-MAPseq-code

Read it in the paper: doi.org/10.1038/s41586-026-10458-y.

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

1 repository of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
5 scripts, each with its path and the digest of its content;
7 matches between paragraphs of the paper and lines of the code (method lexical-v1);
neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

Datasets cited

doi:10.5061/dryad.8kprr4z2p, at Dryad; found in “Data availability”

Data availability statement

The paper has a data availability statement. Its license (CC BY-NC-ND) does not allow reproducing it here; in short, from what the harvester recognized in it:

it points to a dataset: Dryad 10.5061/dryad.8kprr4z2p

Read it in the paper: doi.org/10.1038/s41586-026-10458-y.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 2, 28 September 2026

Publisher: n/a → Nature Portfolio

Version 1, 28 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 7 authors, 2 keywords, 10 MeSH terms, 1 funder, 75 references.

Cite

This paper

Isko, E. C., Harpole, C. E., Zheng, X. M., Zhan, H., Davis, M. B., Zador, A. M., & Banerjee, A. (2026). Specific expansion of motor cortical projections in a singing mouse. Nature, 655(8122), 438-446. https://doi.org/10.1038/s41586-026-10458-y

BibTeX

@article{isko2026specific, author = {Isko, Emily C. and Harpole, Clifford E. and Zheng, Xiaoyue Mike and Zhan, Huiqing and Davis, Martin B. and Zador, Anthony M. and Banerjee, Arkarup}, title = {{Specific expansion of motor cortical projections in a singing mouse}}, journal = {Nature}, year = {2026}, month = may, volume = {655}, number = {8122}, pages = {438--446}, publisher = {Nature Portfolio}, issn = {0028-0836}, doi = {10.1038/s41586-026-10458-y}, url = {https://doi.org/10.1038/s41586-026-10458-y}, pmid = {42092144}, pmcid = {PMC13345907} }

RIS

TY - JOUR AU - Isko, Emily C. AU - Harpole, Clifford E. AU - Zheng, Xiaoyue Mike AU - Zhan, Huiqing AU - Davis, Martin B. AU - Zador, Anthony M. AU - Banerjee, Arkarup TI - Specific expansion of motor cortical projections in a singing mouse T2 - Nature J2 - Nature PY - 2026 DA - 2026/05/06 VL - 655 IS - 8122 SP - 438 EP - 446 SN - 0028-0836 PB - Nature Portfolio DO - 10.1038/s41586-026-10458-y UR - https://doi.org/10.1038/s41586-026-10458-y LA - en ER -

CSL-JSON

{ "id": "10.1038/s41586-026-10458-y", "type": "article-journal", "title": "Specific expansion of motor cortical projections in a singing mouse", "container-title": "Nature", "author": [ { "family": "Isko", "given": "Emily C." }, { "family": "Harpole", "given": "Clifford E." }, { "family": "Zheng", "given": "Xiaoyue Mike" }, { "family": "Zhan", "given": "Huiqing" }, { "family": "Davis", "given": "Martin B." }, { "family": "Zador", "given": "Anthony M." }, { "family": "Banerjee", "given": "Arkarup" } ], "container-title-short": "Nature", "volume": "655", "issue": "8122", "page": "438-446", "DOI": "10.1038/s41586-026-10458-y", "PMID": "42092144", "PMCID": "PMC13345907", "ISSN": "0028-0836", "publisher": "Nature Portfolio", "URL": "https://doi.org/10.1038/s41586-026-10458-y", "language": "en", "issued": { "date-parts": [ [ 2026, 5, 6 ] ] } }

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.