Processing API

Processing submodule.

smooth_cmd(path, out_dir, sigma=2)

Generate a smoothed DeltaF/F recording using a Laplace of Gaussian filter.

Source code in src/mesoscopy/process/__init__.py

@process_cmd.command("smooth")
@click.argument(
    "path",
    type=click.Path(exists=True),
)
@click.option(
    "-o",
    "--out_dir",
    type=click.Path(dir_okay=True),
    default="./",
    help="Output directory for smoothed recording.",
)
@click.option(
    "-s",
    "--sigma",
    type=int,
    default=2,
    help="Standard deviation of the Gaussian kernel, in pixels.",
)
def smooth_cmd(path: str, out_dir: str, sigma: int = 2) -> None:
    """Generate a smoothed DeltaF/F recording using a Laplace of Gaussian filter."""
    import mesoscopy.process.smooth as psm

    if not Path(out_dir).exists():
        click.echo(f"Creating output directory {out_dir}...")
        Path(out_dir).mkdir(parents=True)

    click.echo(f"Loading preprocessed recording from {path}...")
    # Determine whether we're working with an NWB file
    nwb = bool(path.endswith(".nwb"))
    session_id, deltaf_series, timestamps = io.load_deltaf(path, nwb=nwb)
    aligned = io.read_timestamps_aligned(path)

    outpath = out_dir + os.sep + session_id + "_smoothed.h5"

    with timer.Timer(message="Smoothing with LoG"):
        smoothed_deltaf = psm.laplace_gaussian(deltaf_series, sigma=sigma)

        outpath = io.write_h5(
            path=outpath,
            data={
                "/F": smoothed_deltaf,
                "/timestamps": timestamps,
            },
        )
        if aligned is not None:
            io.write_timestamps_aligned(outpath, *aligned)
    click.echo(f"Saved smoothed recording at {outpath}")

zscore_cmd(path, out_dir)

Pixel-wise z-score ∆F/F signal.

Source code in src/mesoscopy/process/__init__.py

@process_cmd.command("zscore")
@click.argument(
    "path",
    type=click.Path(exists=True),
)
@click.option(
    "-o",
    "--out_dir",
    type=click.Path(dir_okay=True),
    default="./",
    help="Output directory for z-scored recording.",
)
def zscore_cmd(path: str, out_dir: str) -> None:
    """Pixel-wise z-score ∆F/F signal."""
    from pynwb.image import ImageSeries

    import mesoscopy.process.zscore as pzs

    if not Path(out_dir).exists():
        click.echo(f"Creating output directory {out_dir}...")
        Path(out_dir).mkdir(parents=True)

    click.echo(f"Loading preprocessed recording from {path}...")
    # Determine whether we're working with an NWB file
    nwb = bool(path.endswith(".nwb"))
    session_id, deltaf_series, timestamps = io.load_deltaf(path, nwb=nwb)
    aligned = io.read_timestamps_aligned(path)

    h5_outpath = out_dir + os.sep + session_id + "_zscored.h5"
    with timer.Timer(message="Z-scoring DeltaF/F"):
        zscored = pzs.zscore_deltaf(deltaf_series)
        h5_outpath = io.write_h5(
            path=h5_outpath,
            data={
                "/F": zscored,
                "/timestamps": timestamps,
            },
        )
        if aligned is not None:
            io.write_timestamps_aligned(h5_outpath, *aligned)

    click.echo(f"Saved z-scored recording at {h5_outpath}")

    # Append to NWB file
    if nwb:
        click.echo("Appending to NWB file...")
        nwbfile, nwbio = io.read_nwb(path, return_io=True)
        f = io.read_h5(h5_outpath)
        try:
            ophys_module = nwbfile.create_processing_module(
                name="ophys", description="optical physiology processed data"
            )
        except ValueError:
            click.echo("Processing module already exists...")
            ophys_module = nwbfile.processing["ophys"]

        zscored_series = ImageSeries(
            name="zScoredDeltaF",
            data=f["/F"],
            timestamps=f["/timestamps"],
            unit="df/f",
            description="z-scored dF/F widefield cortical imaging series",
        )

        ophys_module.add(zscored_series)

        io.write_nwb(path, nwbfile, io=nwbio)

regions_cmd(path, out_dir, mask_paths, include_aba)

Extract ∆F signal averages from ABA-defined regions, custom region masks, or both.

Source code in src/mesoscopy/process/__init__.py

@process_cmd.command("regions")
@click.argument(
    "path",
    type=click.Path(exists=True),
)
@click.option(
    "-o",
    "--out_dir",
    type=click.Path(dir_okay=True),
    default="./",
    help="Output directory for region activity file.",
)
@click.option(
    "-m",
    "--mask",
    "mask_paths",
    type=click.Path(exists=True),
    multiple=True,
    help=(
        "Path to a custom region mask file (NPY, NPZ or TIFF), drawn in registered frame coordinates. A boolean mask"
        " is extracted as a single region, named after the file (or, for NPZ files, after each key); an"
        " integer-labelled mask is extracted as one region per label, named '<name>_<label>'. Masks are used as"
        " drawn and are not mirrored across hemispheres. May be passed multiple times to combine masks from several"
        " files. Supplying a mask replaces the ABA region extraction unless --include-aba is given."
    ),
)
@click.option(
    "--include-aba",
    is_flag=True,
    default=False,
    help="Extract ABA-defined regions alongside any custom masks. Has no effect if no mask is supplied.",
)
# A click command's docstring doubles as its --help text, so it deliberately carries no Raises: section.
def regions_cmd(path: str, out_dir: str, mask_paths: tuple[str, ...], include_aba: bool) -> None:
    """Extract ∆F signal averages from ABA-defined regions, custom region masks, or both."""  # noqa: DOC501
    import pandas as pd

    import mesoscopy.process.region as pr

    if not Path(out_dir).exists():
        click.echo(f"Creating output directory {out_dir}...")
        Path(out_dir).mkdir(parents=True)

    masks: dict[str, np.ndarray] = {}
    for mask_path in mask_paths:
        click.echo(f"Loading region masks from {mask_path}...")
        try:
            file_masks = io.read_mask(mask_path)
        except ValueError as err:
            msg = f"Could not read region masks from {mask_path}: {err}"
            raise click.ClickException(msg) from err

        duplicates = [name for name in file_masks if name in masks]
        if duplicates:
            msg = (
                f"Region name(s) {', '.join(duplicates)} from {mask_path} are already defined by an earlier mask"
                " file. Rename the mask file, or its NPZ keys, so that every region name is unique."
            )
            raise click.ClickException(msg)

        masks.update(file_masks)

    if masks:
        click.echo(f"Loaded {len(masks)} region mask(s): {', '.join(masks)}")

    click.echo(f"Loading preprocessed recording from {path}...")
    # Determine whether we're working with an NWB file
    nwb = bool(path.endswith(".nwb"))
    session_id, deltaf_series, timestamps = io.load_deltaf(path, nwb=nwb)
    aligned = io.read_timestamps_aligned(path)

    outpath = out_dir + os.sep + session_id + "_regions.csv"

    with timer.Timer(message="Extracting region activity"):
        activity = []
        # Custom masks stand in for the ABA regions unless the ABA regions are explicitly asked for as well.
        if include_aba or not masks:
            activity.append(pd.DataFrame(pr.extract_all_regions(deltaf_series, as_dataframe=True)))
        if masks:
            try:
                activity.append(pd.DataFrame(pr.extract_all_masks(deltaf_series, masks, as_dataframe=True)))
            except ValueError as err:
                msg = str(err)
                raise click.ClickException(msg) from err

        region_activity = pd.concat(activity, ignore_index=True)
        time_idx = region_activity["time_idx"].to_numpy()
        region_activity["time_idx"] = [str(timestamp, encoding="utf-8") for timestamp in timestamps[time_idx]]
        region_activity.rename(columns={"time_idx": "timestamp"}, inplace=True)
        if aligned is not None:
            region_activity.insert(
                region_activity.columns.get_loc("timestamp") + 1, "time_aligned", aligned[0][time_idx]
            )
        region_activity.to_csv(outpath, index=False)

    click.echo(f"Saved region activity at {outpath}")

regression_cmd(recording_path, regressor_path, out_dir, alpha, nuisance_regressor_paths, fast, file_format)

Perform pixel-wise ridge regression on a preprocessed ∆F/F recording.

Source code in src/mesoscopy/process/__init__.py

@process_cmd.command("regression")
@click.argument(
    "recording_path",
    type=click.Path(exists=True),
)
@click.argument(
    "regressor_path",
    type=click.Path(exists=True),
)
@click.option(
    "-o",
    "--out_dir",
    type=click.Path(dir_okay=True),
    default="./",
    help="Output directory for regression results.",
)
@click.option(
    "-a",
    "--alpha",
    type=float,
    default=1.0,
    help="Ridge regularisation strength. Defaults to 1.0.",
)
@click.option(
    "-n",
    "--nuisance-regressors",
    "nuisance_regressor_paths",
    type=click.Path(exists=True),
    multiple=True,
    help=(
        "Path to an external nuisance regressor file (NPZ or HDF5), e.g. behavioural motion energy. Every"
        " array/dataset in the file other than 'timestamps' is treated as one nuisance regressor and interpolated"
        " onto the recording's own timestamps before being z-scored and appended to the regressor matrix. May be"
        " passed multiple times to add nuisance regressors from several files."
    ),
)
@click.option(
    "-f",
    "--fast",
    is_flag=True,
    default=False,
    help=(
        "Use fast vectorised implementation of ridge regression. This is an experimental feature and may not work for"
        " all datasets. Use with caution. Defaults to False."
    ),
)
@click.option(
    "--npz",
    "file_format",
    flag_value="npz",
    default="npz",
    help="Save regression results as a compressed NumPy .npz file. Defaults to True.",
)
@click.option(
    "--h5",
    "file_format",
    flag_value="h5",
    help="Save regression results as an HDF5 file. Defaults to False.",
)
def regression_cmd(
    recording_path: str,
    regressor_path: str,
    out_dir: str,
    alpha: float,
    nuisance_regressor_paths: tuple[str, ...],
    fast: bool,
    file_format: str,
) -> None:
    """Perform pixel-wise ridge regression on a preprocessed ∆F/F recording."""  # noqa: DOC501
    import mesoscopy.process.regression as regr

    if not Path(out_dir).exists():
        click.echo(f"Creating output directory {out_dir}...")
        Path(out_dir).mkdir(parents=True)

    click.echo(f"Loading preprocessed recording from {recording_path}...")
    # Determine whether we're working with an NWB file
    nwb = bool(recording_path.endswith(".nwb"))
    session_id, deltaf_series, timestamps = io.load_deltaf(recording_path, nwb=nwb)

    click.echo(f"Loading regressors from {regressor_path}...")
    regressors, labels, trial_idx, regressor_alignment = io.read_regressors(regressor_path)
    labels = list(labels)

    recording_alignment = io.read_timestamps_aligned(recording_path)
    try:
        for warning in regr.check_alignment(recording_alignment, regressor_alignment):
            click.echo(f"WARNING: {warning}")
    except ValueError as err:
        raise click.ClickException(str(err)) from err

    alignment_attrs: dict[str, str] = {}
    if recording_alignment is not None and regressor_alignment is not None:
        alignment_attrs = {
            "session_start_time": recording_alignment[1]["session_start_time"],
            "behaviour_session": recording_alignment[1]["behaviour_session"],
        }

    if nuisance_regressor_paths:
        # Nuisance regressors are recorded on their own clock (e.g. a behavioural camera), so both series are
        # anchored to elapsed seconds since their own first sample before interpolating one onto the other -- see
        # `regr.elapsed_seconds`.
        target_timestamps = regr.elapsed_seconds(timestamps)

    for nuisance_path in nuisance_regressor_paths:
        click.echo(f"Loading nuisance regressors from {nuisance_path}...")
        nuisance, nuisance_labels, nuisance_timestamps = io.read_nuisance_regressors(nuisance_path)
        nuisance = regr.interpolate_regressors(nuisance, regr.elapsed_seconds(nuisance_timestamps), target_timestamps)
        if trial_idx is not None:
            nuisance = nuisance[trial_idx]
        regressors, labels = regr.append_nuisance_regressors(regressors, labels, nuisance, nuisance_labels)
        click.echo(f"Added nuisance regressors: {nuisance_labels}")

    outpath = out_dir + os.sep + session_id + f"_regression.{file_format}"

    trial_idx_used = False
    if trial_idx is not None:
        deltaf_series = deltaf_series[trial_idx]
        trial_idx_used = True

    with timer.Timer(message="Running regression"):
        if fast:
            coefs, r2, mse = regr.ridge_regression_fast(deltaf_series, regressors, alpha=alpha)
        else:
            coefs, r2, mse = regr.ridge_regression(deltaf_series, regressors)

        if file_format == "npz":
            outpath = io.write_npz(
                path=outpath,
                data={
                    "coefficients": coefs,
                    "r2": r2,
                    "mse": mse,
                    "labels": labels,
                    "trial_idx": trial_idx if trial_idx_used else [],
                    **alignment_attrs,
                },
            )
        elif file_format == "h5":
            outpath = io.write_h5(
                path=outpath,
                data={
                    "/coefficients": coefs,
                    "/r2": r2,
                    "/mse": mse,
                    "/labels": np.array(labels).astype("S"),
                    "/trial_idx": trial_idx if trial_idx_used else [],
                },
                attributes=alignment_attrs,
            )
    click.echo(f"Saved regression results at {outpath}")

perievent_cmd(recording_path, trials_path, out_dir, event, pre, post, fs, method, baseline, with_metrics, response, onset, onset_sd, onset_fraction, onset_min_samples, extrapolate_range, smoothing, decay_fraction, min_trials, reliability, bootstrap)

Extract per-trial windows around a behavioural event from a behaviour-aligned recording.

RECORDING_PATH is an HDF5 recording with /F and /timestamps_aligned, or a _regions.csv with a time_aligned column. TRIALS_PATH is the *_trials.csv written by visiomode-analysis session. HDF5 input gives HDF5 output; CSV input gives long-format CSV output with the TRIALS_PATH columns joined onto every row. --with-metrics also writes the tables of process metrics from the CSV output, using --baseline as the metrics baseline window.

Source code in src/mesoscopy/process/__init__.py

@process_cmd.command("peri-event")
@click.argument(
    "recording_path",
    type=click.Path(exists=True, dir_okay=False),
)
@click.argument(
    "trials_path",
    type=click.Path(exists=True, dir_okay=False),
)
@click.option(
    "-o",
    "--out_dir",
    type=click.Path(dir_okay=True),
    default="./",
    help="Output directory for peri-event windows.",
)
@click.option(
    "--event",
    type=click.Choice(pev.EVENTS),
    default="cue_onset",
    show_default=True,
    help="Behavioural event to align windows to.",
)
@click.option(
    "--pre",
    type=float,
    default=1.0,
    show_default=True,
    help="Seconds before the event.",
)
@click.option(
    "--post",
    type=float,
    default=3.0,
    show_default=True,
    help="Seconds after the event.",
)
@click.option(
    "--fs",
    type=float,
    default=25.0,
    show_default=True,
    help="Sampling rate of the output windows, in Hz.",
)
@click.option(
    "--method",
    type=click.Choice(["interp", "nearest"]),
    default="interp",
    show_default=True,
    help="Linear interpolation, or the nearest recorded sample.",
)
@click.option(
    "--baseline",
    type=(float, float),
    default=None,
    help="Baseline window START END, seconds relative to the event; its mean is subtracted per trial.",
)
@click.option(
    "--with-metrics",
    is_flag=True,
    default=False,
    help="Also write the response metrics of `process metrics`, joined with TRIALS_PATH. CSV input only.",
)
@_metrics_options
@_reliability_options
@_bootstrap_options
def perievent_cmd(
    recording_path: str,
    trials_path: str,
    out_dir: str,
    event: str,
    pre: float,
    post: float,
    fs: float,
    method: str,
    baseline: tuple[float, float] | None,
    with_metrics: bool,
    response: tuple[float, float] | None,
    onset: str,
    onset_sd: float,
    onset_fraction: float,
    onset_min_samples: int,
    extrapolate_range: tuple[float, float],
    smoothing: int,
    decay_fraction: float,
    min_trials: int,
    reliability: pm.ReliabilityOptions,
    bootstrap: pm.BootstrapOptions,
) -> None:
    """Extract per-trial windows around a behavioural event from a behaviour-aligned recording.

    RECORDING_PATH is an HDF5 recording with /F and /timestamps_aligned, or a _regions.csv with a time_aligned
    column. TRIALS_PATH is the *_trials.csv written by visiomode-analysis session. HDF5 input gives HDF5 output;
    CSV input gives long-format CSV output with the TRIALS_PATH columns joined onto every row. --with-metrics also
    writes the tables of `process metrics` from the CSV output, using --baseline as the metrics baseline window.
    """  # noqa: DOC501
    import pandas as pd

    is_csv = recording_path.endswith(".csv")
    if with_metrics:
        if not is_csv:
            msg = "needs a _regions.csv recording; metrics are not computed for HDF5 windows."
            raise click.BadParameter(msg, param_hint="--with-metrics")
        _validate_metrics_options(smoothing, extrapolate_range)

    if not Path(out_dir).exists():
        click.echo(f"Creating output directory {out_dir}...")
        Path(out_dir).mkdir(parents=True)

    click.echo(f"Loading trials from {trials_path}...")
    trials = pd.read_csv(trials_path)
    events, trial_index = pev.event_times(trials, event)
    grid = pev.window_grid(pre, post, fs)

    click.echo(f"Loading recording from {recording_path}...")
    session_id = Path(recording_path).stem
    outpath = out_dir + os.sep + session_id + f"_event-{pev.event_name(event)}_perievent." + ("csv" if is_csv else "h5")

    with timer.Timer(message="Extracting peri-event windows"):
        if is_csv:
            kept, windows = _perievent_csv(recording_path, outpath, events, trial_index, grid, method, baseline, trials)
        else:
            attrs = {
                "event": event,
                "pre": pre,
                "post": post,
                "fs": fs,
                "method": method,
                "baseline": np.array(baseline if baseline is not None else [], dtype=np.float64),
                "trials_path": str(trials_path),
            }
            kept = _perievent_h5(recording_path, outpath, events, trial_index, grid, method, baseline, attrs)

    n_kept = int(kept.sum())
    click.echo(f"Kept {n_kept} trials, dropped {len(trials) - n_kept}.")
    click.echo(f"Saved peri-event windows at {outpath}")

    if with_metrics and n_kept == 0:
        click.echo("No trials kept, skipping metrics.")
    elif with_metrics:
        with timer.Timer(message="Extracting metrics"), warnings.catch_warnings(record=True) as caught:
            warnings.simplefilter("always", UserWarning)
            tables = pm.metrics_tables(
                windows,
                baseline=baseline,
                response=response,
                onset=onset,
                onset_sd=onset_sd,
                onset_fraction=onset_fraction,
                onset_min_samples=onset_min_samples,
                extrapolate_range=extrapolate_range,
                decay_fraction=decay_fraction,
                smoothing=smoothing,
                event=event,
                reliability=reliability,
                min_trials=min_trials,
                bootstrap=bootstrap,
            )
        _echo_warnings(caught)
        _write_metrics(tables, out_dir, Path(outpath).stem.removesuffix("_perievent"))

metrics_cmd(path, out_dir, baseline, trials_path, response, onset, onset_sd, onset_fraction, onset_min_samples, extrapolate_range, smoothing, decay_fraction, min_trials, reliability, bootstrap)

Extract per-trial response metrics and per-session variability from peri-event traces.

PATH is the long-format *_perievent.csv written by process peri-event. Writes _metrics.csv with onset time, peak time, amplitude, area under the curve, decay time, offset time and duration per trial per region, and _metrics-session.csv with the mean, SD and CV of each across trials, plus the mean pairwise trial-trace correlation, per region and trial group. Go/no-go sessions are grouped by sdt_type, go/no-go stimulus and lever push, as well as all trials. For go/no-go sessions aligned to the cue, trial start or lever push, the session table also has trial-to-trial reliability metrics over the cue-to-lever epoch.

Also writes _metrics-boot.csv with the same metrics taken from the trial-mean trace of each region and trial group, with bootstrap percentile intervals from resampling the trials, and _traces-boot.csv with the mean traces and their intervals. --bootstrap 0 skips both.

Source code in src/mesoscopy/process/__init__.py

@process_cmd.command("metrics")
@click.argument(
    "path",
    type=click.Path(exists=True, dir_okay=False),
)
@click.option(
    "-o",
    "--out_dir",
    type=click.Path(dir_okay=True),
    default="./",
    help="Output directory for metrics tables.",
)
@click.option(
    "--baseline",
    type=(float, float),
    default=None,
    help="Baseline window START END, seconds relative to the event, end exclusive. Defaults to all pre-event samples.",
)
@click.option(
    "-t",
    "--trials",
    "trials_path",
    type=click.Path(exists=True, dir_okay=False),
    default=None,
    help="Trials CSV to join onto the per-trial table by trial_index, for peri-event files without trials columns.",
)
@_metrics_options
@_reliability_options
@_bootstrap_options
def metrics_cmd(
    path: str,
    out_dir: str,
    baseline: tuple[float, float] | None,
    trials_path: str | None,
    response: tuple[float, float] | None,
    onset: str,
    onset_sd: float,
    onset_fraction: float,
    onset_min_samples: int,
    extrapolate_range: tuple[float, float],
    smoothing: int,
    decay_fraction: float,
    min_trials: int,
    reliability: pm.ReliabilityOptions,
    bootstrap: pm.BootstrapOptions,
) -> None:
    """Extract per-trial response metrics and per-session variability from peri-event traces.

    PATH is the long-format *_perievent.csv written by `process peri-event`. Writes <stem>_metrics.csv with onset
    time, peak time, amplitude, area under the curve, decay time, offset time and duration per trial per region,
    and <stem>_metrics-session.csv with the mean, SD and CV of each across trials, plus the mean pairwise
    trial-trace correlation, per region and trial group. Go/no-go sessions are grouped by sdt_type, go/no-go
    stimulus and lever push, as well as all trials. For go/no-go sessions aligned to the cue, trial start or lever
    push, the session table also has trial-to-trial reliability metrics over the cue-to-lever epoch.

    Also writes <stem>_metrics-boot.csv with the same metrics taken from the trial-mean trace of each region and
    trial group, with bootstrap percentile intervals from resampling the trials, and <stem>_traces-boot.csv with
    the mean traces and their intervals. --bootstrap 0 skips both.
    """  # noqa: DOC501
    import pandas as pd

    _validate_metrics_options(smoothing, extrapolate_range)

    if not Path(out_dir).exists():
        click.echo(f"Creating output directory {out_dir}...")
        Path(out_dir).mkdir(parents=True)

    click.echo(f"Loading peri-event traces from {path}...")
    perievent = pd.read_csv(path)
    if perievent.empty:
        msg = f"{path} has no rows; `process peri-event` kept no trials."
        raise click.ClickException(msg)
    trials = None
    if trials_path is not None:
        click.echo(f"Loading trials from {trials_path}...")
        trials = pd.read_csv(trials_path)

    with timer.Timer(message="Extracting metrics"), warnings.catch_warnings(record=True) as caught:
        warnings.simplefilter("always", UserWarning)
        try:
            tables = pm.metrics_tables(
                perievent,
                baseline=baseline,
                response=response,
                trials=trials,
                onset=onset,
                onset_sd=onset_sd,
                onset_fraction=onset_fraction,
                onset_min_samples=onset_min_samples,
                extrapolate_range=extrapolate_range,
                decay_fraction=decay_fraction,
                smoothing=smoothing,
                reliability=reliability,
                min_trials=min_trials,
                bootstrap=bootstrap,
            )
        except ValueError as error:
            msg = f"{path}: {error} Expected the columns written by `process peri-event`."
            raise click.ClickException(msg) from error
    _echo_warnings(caught)

    _write_metrics(tables, out_dir, Path(path).stem.removesuffix("_perievent"))

connectivity_cmd(path, out_dir, response, trials_path, min_trials, mask_response, min_rt, response_pad, max_lag, mi_neighbours, mutual_info, te_history, te_surrogates, transfer_entropy, rolling, rolling_window, rolling_step, seed)

Measure pairwise connectivity between regions, over the cue-to-response epoch or over whole traces.

PATH is the long-format _perievent.csv written by process peri-event, or the _regions.csv written by process regions. Writes _connectivity.csv with one row per region pair and trial group: the Pearson correlation of the pooled epoch samples, the same after removing the group's mean response, the mean of the per-trial correlations, the mutual information of the pooled samples and of the residuals in bits, the partial correlation given every other region, the lag and value of the peak cross-correlation, and with --with-te the transfer entropy of the residuals in each direction with its z-score against surrogates. Epochs run from the cue to each trial's response, as for the reliability metrics of process metrics, extended past the response by --response-pad; without the cue_onset and response_time columns the whole response window is used. Go/no-go sessions are grouped by sdt_type, go/no-go stimulus and lever push, as well as all trials. With --with-rolling it also writes _connectivity-rolling.csv, with the correlations and mutual information in windows sliding over the whole peri-event window, one row per region pair, trial group and window.

A regions CSV gives the metrics over the whole recording instead, as the all-trials rows only, without the residual and per-trial columns. The trial options and --with-rolling do not apply.

Source code in src/mesoscopy/process/__init__.py

@process_cmd.command("connectivity")
@click.argument(
    "path",
    type=click.Path(exists=True, dir_okay=False),
)
@click.option(
    "-o",
    "--out_dir",
    type=click.Path(dir_okay=True),
    default="./",
    help="Output directory for the connectivity table.",
)
@click.option(
    "--response",
    type=(float, float),
    default=None,
    help="Response window START END, seconds relative to the event, end inclusive. Defaults to all post-event samples."
    " Peri-event input only.",
)
@click.option(
    "-t",
    "--trials",
    "trials_path",
    type=click.Path(exists=True, dir_okay=False),
    default=None,
    help="Trials CSV to join by trial_index, for peri-event files without trials columns. Peri-event input only.",
)
@click.option(
    "--min-trials",
    type=click.IntRange(min=1),
    default=10,
    show_default=True,
    help="Trial groups with fewer trials get empty metrics. The all-trials metrics are always taken. Peri-event"
    " input only.",
)
@click.option(
    "--mask-response/--no-mask-response",
    default=True,
    show_default=True,
    help="End cue-aligned epochs at each trial's response; trials without one end at the median response time."
    " Peri-event input only.",
)
@click.option(
    "--min-rt",
    type=click.FloatRange(min=0),
    default=0.2,
    show_default=True,
    help="Drop trials with a response time below this, in seconds. Peri-event input only.",
)
@click.option(
    "--response-pad",
    type=click.FloatRange(min=0),
    default=0.0,
    show_default=True,
    help="Extend each epoch past the response by this, in seconds. Cue-aligned epochs still end within the response"
    " window. Peri-event input only.",
)
@click.option(
    "--max-lag",
    type=click.FloatRange(min=0),
    default=0.5,
    show_default=True,
    help="Largest lag searched for the peak cross-correlation, in seconds.",
)
@click.option(
    "--mi-neighbours",
    type=click.IntRange(min=1),
    default=3,
    show_default=True,
    help="Neighbours for the mutual information estimator.",
)
@click.option(
    "--with-mi/--no-mi",
    "mutual_info",
    default=True,
    show_default=True,
    help="Estimate the mutual information of each pair, the slowest part.",
)
@click.option(
    "--te-history",
    type=click.IntRange(min=1),
    default=1,
    show_default=True,
    help="Past samples of the target the transfer entropy conditions on.",
)
@click.option(
    "--te-surrogates",
    type=click.IntRange(min=0),
    default=200,
    show_default=True,
    help="Surrogates for the transfer entropy z-scores: trial shuffles on peri-event input, circular shifts on a"
    " regions CSV. 0 leaves the z-scores out.",
)
@click.option(
    "--with-te/--no-te",
    "transfer_entropy",
    default=False,
    show_default=True,
    help="Estimate the transfer entropy of each pair in both directions.",
)
@click.option(
    "--with-rolling/--no-rolling",
    "rolling",
    default=False,
    show_default=True,
    help="Also take r, r_residual, mi and mi_residual in windows sliding over the whole peri-event window, written"
    " to <stem>_connectivity-rolling.csv. Peri-event input only.",
)
@click.option(
    "--rolling-window",
    type=click.FloatRange(min=0, min_open=True),
    default=0.5,
    show_default=True,
    help="Length of the rolling windows, in seconds, rounded to whole samples.",
)
@click.option(
    "--rolling-step",
    type=click.FloatRange(min=0, min_open=True),
    default=0.1,
    show_default=True,
    help="Time between the starts of consecutive rolling windows, in seconds, rounded to whole samples.",
)
@click.option(
    "--seed",
    type=int,
    default=42,
    show_default=True,
    help="Random seed for the surrogates.",
)
def connectivity_cmd(
    path: str,
    out_dir: str,
    response: tuple[float, float] | None,
    trials_path: str | None,
    min_trials: int,
    mask_response: bool,
    min_rt: float,
    response_pad: float,
    max_lag: float,
    mi_neighbours: int,
    mutual_info: bool,
    te_history: int,
    te_surrogates: int,
    transfer_entropy: bool,
    rolling: bool,
    rolling_window: float,
    rolling_step: float,
    seed: int,
) -> None:
    """Measure pairwise connectivity between regions, over the cue-to-response epoch or over whole traces.

    PATH is the long-format *_perievent.csv written by `process peri-event`, or the *_regions.csv written by
    `process regions`. Writes <stem>_connectivity.csv with one row per region pair and trial group: the Pearson
    correlation of the pooled epoch samples, the same after removing the group's mean response, the mean of the
    per-trial correlations, the mutual information of the pooled samples and of the residuals in bits, the
    partial correlation given every other region, the lag and value of the peak cross-correlation, and with --with-te
    the transfer entropy of the residuals in each direction with its z-score against surrogates. Epochs run from
    the cue to each trial's response, as for the reliability metrics of `process metrics`, extended past the
    response by --response-pad; without the cue_onset and response_time columns the whole response window is used.
    Go/no-go sessions are grouped by sdt_type, go/no-go stimulus and lever push, as well as all trials. With
    --with-rolling it also writes <stem>_connectivity-rolling.csv, with the correlations and mutual information in
    windows sliding over the whole peri-event window, one row per region pair, trial group and window.

    A regions CSV gives the metrics over the whole recording instead, as the all-trials rows only, without the
    residual and per-trial columns. The trial options and --with-rolling do not apply.
    """  # noqa: DOC501
    import pandas as pd

    if not Path(out_dir).exists():
        click.echo(f"Creating output directory {out_dir}...")
        Path(out_dir).mkdir(parents=True)

    whole_trace = Path(path).stem.endswith("_regions")
    if rolling and whole_trace:
        click.echo("Warning: --with-rolling takes a peri-event file; the rolling table is skipped.")
        rolling = False
    click.echo(f"Loading {'region' if whole_trace else 'peri-event'} traces from {path}...")
    table_in = pd.read_csv(path)
    if table_in.empty:
        msg = f"{path} has no rows; `process {'regions' if whole_trace else 'peri-event'}` kept no frames."
        raise click.ClickException(msg)
    trials = None
    if trials_path is not None and not whole_trace:
        click.echo(f"Loading trials from {trials_path}...")
        trials = pd.read_csv(trials_path)

    with timer.Timer(message="Extracting connectivity"), warnings.catch_warnings(record=True) as caught:
        warnings.simplefilter("always", UserWarning)
        try:
            options: dict[str, typing.Any] = {
                "max_lag": max_lag,
                "mi_neighbours": mi_neighbours,
                "mutual_info": mutual_info,
                "te_history": te_history,
                "te_surrogates": te_surrogates,
                "transfer_entropy": transfer_entropy,
                "seed": seed,
            }
            if whole_trace:
                table = pc.trace_table(table_in, **options)
            else:
                table = pc.connectivity_table(
                    table_in,
                    response=response,
                    trials=trials,
                    mask_response=mask_response,
                    min_rt=min_rt,
                    response_pad=response_pad,
                    min_trials=min_trials,
                    **options,
                )
        except ValueError as error:
            stage = "regions" if whole_trace else "peri-event"
            msg = f"{path}: {error} Expected the columns written by `process {stage}`."
            raise click.ClickException(msg) from error
    _echo_warnings(caught)

    stem = out_dir + os.sep + Path(path).stem.removesuffix("_perievent")
    output = stem + "_connectivity.csv"
    table.to_csv(output, index=False)
    click.echo(f"Saved connectivity metrics at {output}")
    if not rolling:
        return

    echoed = {str(warning.message) for warning in caught}
    with timer.Timer(message="Extracting rolling connectivity"), warnings.catch_warnings(record=True) as caught:
        warnings.simplefilter("always", UserWarning)
        try:
            rolling_table = pc.rolling_table(
                table_in,
                window=rolling_window,
                step=rolling_step,
                trials=trials,
                min_rt=min_rt,
                min_trials=min_trials,
                mi_neighbours=mi_neighbours,
            )
        except ValueError as error:
            msg = f"{path}: {error} Check --rolling-window and --rolling-step."
            raise click.ClickException(msg) from error
    _echo_warnings([warning for warning in caught if str(warning.message) not in echoed])

    output = stem + "_connectivity-rolling.csv"
    rolling_table.to_csv(output, index=False)
    click.echo(f"Saved rolling connectivity metrics at {output}")

decode_cmd(path, out_dir, label, response, trials_path, min_trials, mask_response, per_trial_end, min_rt, response_pad, folds, repeats, shuffles, rolling, rolling_window, rolling_step, baseline, seed)

Decode the stimulus or the lever push of each trial from region activity.

PATH is the long-format *_perievent.csv written by process peri-event for a go/no-go session. A logistic regression and a linear discriminant decoder are fitted to each region's mean activity between the cue and the median response time, and to all regions together, under repeated stratified cross-validation, over all trials and within the groups that hold the other variable constant. Writes _label-_decoding.csv with the scores and their chance level per decoder, region and trial group, and _label-_decoding-weights.csv with the population weights. With --with-rolling the scores are also taken in windows sliding over the peri-event window and written to _label-_decoding-rolling.csv.

Source code in src/mesoscopy/process/__init__.py

@process_cmd.command("decode")
@click.argument(
    "path",
    type=click.Path(exists=True, dir_okay=False),
)
@click.option(
    "-o",
    "--out_dir",
    type=click.Path(dir_okay=True),
    default="./",
    help="Output directory for the decoding tables.",
)
@click.option(
    "-l",
    "--label",
    type=click.Choice(list(pdc.LABELS), case_sensitive=False),
    default="stim",
    show_default=True,
    help="What to decode. stim is the go against the no-go stimulus, response the lever push against no push.",
)
@click.option(
    "--response",
    type=(float, float),
    default=None,
    help="Response window START END, seconds relative to the event, end inclusive. Defaults to all post-event samples.",
)
@click.option(
    "-t",
    "--trials",
    "trials_path",
    type=click.Path(exists=True, dir_okay=False),
    default=None,
    help="Trials CSV to join by trial_index, for peri-event files without trials columns.",
)
@click.option(
    "--min-trials",
    type=click.IntRange(min=1),
    default=10,
    show_default=True,
    help="Trials each class needs in a trial group. Groups with fewer get empty scores.",
)
@click.option(
    "--mask-response/--no-mask-response",
    default=True,
    show_default=True,
    help="End the epoch at the response rather than at the response window end.",
)
@click.option(
    "--per-trial-end",
    is_flag=True,
    default=False,
    help="End each trial's epoch at its own response instead of the median response time.",
)
@click.option(
    "--min-rt",
    type=click.FloatRange(min=0),
    default=0.2,
    show_default=True,
    help="Drop trials with a response time below this, in seconds.",
)
@click.option(
    "--response-pad",
    type=click.FloatRange(min=0),
    default=0.0,
    show_default=True,
    help="Extend the epoch past the response by this, in seconds, within the response window.",
)
@click.option(
    "--folds",
    type=click.IntRange(min=2),
    default=5,
    show_default=True,
    help="Folds of the stratified cross-validation.",
)
@click.option(
    "--repeats",
    type=click.IntRange(min=1),
    default=10,
    show_default=True,
    help="Repeats of the cross-validation, each with a fresh split.",
)
@click.option(
    "--shuffles",
    type=click.IntRange(min=0),
    default=200,
    show_default=True,
    help="Label shuffles for the chance level and p-values. 0 leaves the columns out.",
)
@click.option(
    "--with-rolling/--no-rolling",
    "rolling",
    default=False,
    show_default=True,
    help="Also decode in windows sliding over the whole peri-event window, written to"
    " <stem>_label-<label>_decoding-rolling.csv.",
)
@click.option(
    "--rolling-window",
    type=click.FloatRange(min=0, min_open=True),
    default=0.5,
    show_default=True,
    help="Length of the rolling windows, in seconds, rounded to whole samples.",
)
@click.option(
    "--rolling-step",
    type=click.FloatRange(min=0, min_open=True),
    default=0.1,
    show_default=True,
    help="Time between the starts of consecutive rolling windows, in seconds, rounded to whole samples.",
)
@click.option(
    "--baseline",
    type=(float, float),
    default=None,
    help="Window START END, seconds relative to the event, end exclusive, whose mean gives the chance level of the"
    " rolling table. Defaults to every sample before the event.",
)
@click.option(
    "--seed",
    type=int,
    default=42,
    show_default=True,
    help="Random seed for the splits and the shuffles.",
)
def decode_cmd(
    path: str,
    out_dir: str,
    label: str,
    response: tuple[float, float] | None,
    trials_path: str | None,
    min_trials: int,
    mask_response: bool,
    per_trial_end: bool,
    min_rt: float,
    response_pad: float,
    folds: int,
    repeats: int,
    shuffles: int,
    rolling: bool,
    rolling_window: float,
    rolling_step: float,
    baseline: tuple[float, float] | None,
    seed: int,
) -> None:
    """Decode the stimulus or the lever push of each trial from region activity.

    PATH is the long-format *_perievent.csv written by `process peri-event` for a go/no-go session. A logistic
    regression and a linear discriminant decoder are fitted to each region's mean activity between the cue and the
    median response time, and to all regions together, under repeated stratified cross-validation, over all trials
    and within the groups that hold the other variable constant. Writes <stem>_label-<label>_decoding.csv with the
    scores and their chance level per decoder, region and trial group, and <stem>_label-<label>_decoding-weights.csv
    with the population weights. With --with-rolling the scores are also taken in windows sliding over the
    peri-event window and written to <stem>_label-<label>_decoding-rolling.csv.
    """  # noqa: DOC501
    import pandas as pd

    label = label.lower()
    if not Path(out_dir).exists():
        click.echo(f"Creating output directory {out_dir}...")
        Path(out_dir).mkdir(parents=True)

    click.echo(f"Loading peri-event traces from {path}...")
    perievent = pd.read_csv(path)
    if perievent.empty:
        msg = f"{path} has no rows; `process peri-event` kept no trials."
        raise click.ClickException(msg)
    trials = None
    if trials_path is not None:
        click.echo(f"Loading trials from {trials_path}...")
        trials = pd.read_csv(trials_path)

    with timer.Timer(message="Decoding"), warnings.catch_warnings(record=True) as caught:
        warnings.simplefilter("always", UserWarning)
        try:
            tables = pdc.decoding_table(
                perievent,
                label=label,
                response=response,
                trials=trials,
                mask_response=mask_response,
                per_trial_end=per_trial_end,
                min_rt=min_rt,
                response_pad=response_pad,
                min_trials=min_trials,
                folds=folds,
                repeats=repeats,
                shuffles=shuffles,
                seed=seed,
            )
        except ValueError as error:
            msg = f"{path}: {error} Expected the columns written by `process peri-event` with the trials columns."
            raise click.ClickException(msg) from error
    _echo_warnings(caught)

    stem = out_dir + os.sep + Path(path).stem.removesuffix("_perievent") + f"_label-{label}"
    output = stem + "_decoding.csv"
    tables.table.to_csv(output, index=False)
    click.echo(f"Saved decoding scores at {output}")
    output = stem + "_decoding-weights.csv"
    tables.weights.to_csv(output, index=False)
    click.echo(f"Saved population weights at {output}")
    if not rolling:
        return

    echoed = {str(warning.message) for warning in caught}
    with timer.Timer(message="Decoding in rolling windows"), warnings.catch_warnings(record=True) as caught:
        warnings.simplefilter("always", UserWarning)
        try:
            rolling_table = pdc.rolling_table(
                perievent,
                label=label,
                trials=trials,
                window=rolling_window,
                step=rolling_step,
                baseline=baseline,
                min_rt=min_rt,
                min_trials=min_trials,
                folds=folds,
                repeats=repeats,
                shuffles=shuffles,
                seed=seed,
            )
        except ValueError as error:
            msg = f"{path}: {error} Check --rolling-window, --rolling-step and --baseline."
            raise click.ClickException(msg) from error
    _echo_warnings([warning for warning in caught if str(warning.message) not in echoed])

    output = stem + "_decoding-rolling.csv"
    rolling_table.to_csv(output, index=False)
    click.echo(f"Saved rolling decoding scores at {output}")