图宾根神经数据科学笔记-全-
图宾根神经数据科学笔记（全）
001：峰值检测与特征提取 🧠
在本节课中，我们将学习神经数据科学中的第一个核心问题：如何从记录的神经元电信号中检测出动作电位（即“峰值”），并从中提取关键特征，为后续的神经元分类（即“峰值排序”）做准备。
概述
神经元电活动的记录方法可以追溯到约70-80年前，由David Hubel和Torsten Wiesel等先驱开发。他们发明了钨丝电极和尖锐电极，能够非常接近大脑中的神经元，记录其膜电位的变化。如今，记录技术已大大发展，但核心目标不变：从复杂的原始信号中识别出单个神经元的电活动。
本课程将聚焦于在较短时间内（毫秒级）和较小空间尺度上记录神经元活动的技术，例如使用微电极阵列记录单个动作电位（峰值）。我们将从四通道的微丝电极阵列数据入手，学习如何从原始信号中检测峰值并提取特征。
记录技术与数据
记录大脑活动的方法多种多样，它们在时间和空间尺度上差异很大。例如，fMRI和PET成像覆盖的是相对较大尺度的活动。本课程主要讨论的是能够监测单个神经元快速活动（即动作电位或峰值）的技术。
- 历史技术：早期的钨丝微电极结构简单，由一个被绝缘体包裹的导体构成，仅在尖端暴露以测量周围电活动。
- 微丝电极阵列：这是钨丝电极的进一步发展。它由四根微丝捆绑而成，每根微丝都是一个独立的记录通道。由于四根微丝在空间上的排列，附近每个神经元的电压信号在每个通道上的特征会略有不同。我们将利用这一原理进行峰值排序。
- 硅探针与Neuropixels：现代技术如硅探针，尤其是Neuropixels探针，实现了工程学上的飞跃。它们能在单个探针上集成成百上千个记录位点，从而同时记录数百个神经元的活动。
在本课程中，我们将使用微丝电极阵列获取的数据。我们的目标是：从类似下图的四通道同步记录数据中，找出有多少个以及哪些神经元贡献了动作电位。
我们的最终目标是将原始数据转化为清晰标注的数据集，其中包含已识别出的不同神经元对应的峰值波形。
微丝电极阵列峰值排序的原理是三角测量。我们利用的是：电极附近的神经元距离某些通道更近，而距离其他通道较远。因此，不同神经元在不同通道上产生的峰值幅度特征不同，这将成为我们区分它们的依据。
峰值检测流程
峰值排序的任务流程是：从原始数据出发，检测峰值 -> 提取特征 -> 对峰值进行聚类（即排序到不同的神经元）-> 验证聚类结果的生物学合理性。
让我们从第一步——峰值检测开始。我们拥有的原始信号如下图所示，其中包含我们想要的峰值（向下的小脉冲），也包含非常缓慢的低频波动（局部场电位）以及可能的高频噪声。
为了更清晰地看到峰值，我们需要滤除这些干扰。我们的目标是得到下图这样的信号，其中峰值清晰可见。
信号滤波
我们通过滤波来去除不需要的频率成分。以下是对时间序列滤波的简要介绍。
- 低通滤波：去除高频成分，保留低频趋势。其传递函数显示，低频部分增益高，高频部分增益低。
- 高通滤波：去除低频成分，保留高频细节。其传递函数与低通滤波相反。
- 带通滤波：只允许特定频率范围内的信号通过。这对于峰值检测非常有用，因为神经元动作电位的主要能量通常集中在1000 Hz到3000 Hz之间。
我们使用巴特沃斯滤波器来实现带通滤波。这种滤波器的传递函数在通带内尽可能平坦，在截止频率处有较陡的衰减。
在代码练习中，你可以尝试调整这些滤波范围，并查阅文献了解常用的参数。
应用带通滤波后，我们得到主要包含峰值成分的信号。
设置检测阈值
接下来，我们需要设置一个阈值，当信号幅度超过（通常是低于，因为记录到的峰值常为负向）此阈值时，我们认为检测到了一个峰值。
一个常见的方法是估计背景噪声的标准差，然后将阈值设为该标准差的若干倍（例如3或4倍）。但问题在于，如果数据中包含大量峰值，会使得标准差的估计产生偏差。
因此，我们使用一个更稳健的估计量，公式如下：
稳健标准差估计值 = median(|X - median(X)|) / 0.6745
这个公式的推导基于数据服从高斯分布的假设。中位数绝对偏差（MAD）除以0.6745，可以得到对标准差稳健的估计。
设置好阈值后，我们标记所有幅度超过阈值的时刻，将其视为峰值出现的时间点。这样，我们就从一个30 kHz采样的连续记录，得到了一系列峰值出现的时间点列表。
特征提取
检测到峰值后，我们需要提取每个峰值在四个通道上的波形。对于每个峰值时间点，我们截取它周围一小段时间窗口内的信号（例如，峰值前后各0.5毫秒）。
在现代电极（如Neuropixels）记录中，数据维度极高（上百个通道）。为了可视化，常将多通道信号随时间的变化绘制成图像，从中可以清晰看到不同神经元的“空间足迹”。
如果我们直接使用所有数据（例如，4通道 × 32个时间点 = 128维）进行后续分析，可能会面临计算复杂、内存消耗大、模型参数过多等问题。因此，我们需要进行降维，提取一小部分能够稳健地区分不同神经元的特征。
好的特征应该能清晰地区分不同神经元的峰值。过去人们使用过一些描述性特征，如峰峰值、波宽、能量等，但这些特征可能相关性强、分布非高斯，有时甚至无法定义。
主成分分析
一种常用的降维方法是主成分分析。PCA试图找到数据中方差最大的方向，并以此构建一组正交（不相关）的新基。第一个主成分是方差最大的方向，第二个是与第一个正交且方差次大的方向，依此类推。
其目标可以形式化地表示为寻找投影方向，以最大化投影后数据的方差，或最小化重建误差。
通过求解数据协方差矩阵的特征值和特征向量，我们可以得到主成分。对于峰值数据，使用前几个主成分通常就能重构出大部分方差（例如，前3个主成分可解释95%的方差），更高阶的成分可能主要包含噪声。
然而，PCA的局限性在于：方差最大的方向不一定是最有利于聚类分离的方向。有时，方差较小的方向反而能更好地区分不同类别的数据点。
非线性特征提取方法
为了克服PCA的线性限制，可以使用非线性特征提取方法：
- 核主成分分析：通过非线性核函数计算数据间的相似性（核矩阵），然后对核矩阵进行特征分解。这能够处理原始空间中非线性结构的数据（例如同心圆分布）。
- 自编码器网络：这是一种神经网络结构，包含一个将输入压缩到低维表示的编码器和一个从低维表示重建输入的解码器。通过训练网络最小化重建误差，编码器输出的低维表示就是提取的特征。线性自编码器等价于PCA，而使用非线性激活函数的深度自编码器能够学习更复杂的非线性特征表示。
总结
本节课我们一起学习了神经数据科学中峰值处理流程的前两步：
- 峰值检测：我们了解了如何对原始神经电信号进行带通滤波，以突出峰值成分；并学习了如何稳健地估计噪声标准差，从而设置阈值来检测峰值。
- 特征提取：我们探讨了为何需要降维，并介绍了主成分分析这一线性降维方法的基本原理及其局限性。最后，我们简要提及了核PCA和自编码器等非线性特征提取技术，它们能捕获数据中更复杂的结构。
至此，我们已经完成了从原始数据到峰值检测和特征提取的流程。在接下来的课程中，我们将基于这些提取出的特征，进入峰值聚类阶段，学习如何将峰值归类到不同的神经元。
002：尖峰分类 🧠
在本节课中，我们将学习神经数据科学流程中的核心环节——尖峰分类。我们将探讨如何将从原始数据中检测并提取出的尖峰特征，通过聚类算法归类到不同的神经元上。
上一周，我们从原始数据出发，完成了尖峰检测和特征提取。本节中，我们将聚焦于尖峰分类或聚类过程。我们的目标是，将上周从单个神经元提取出的不同尖峰或动作电位的特征，按照它们所对应的神经元进行分组。
这意味着，我们需要将四通道记录中这些未分类的尖峰波形，转化为清晰分离的类别，如下图所示。例如，这里的紫色神经元在第3通道上有高信号，而黄色簇在所有通道上的信号都低得多。那么，我们具体该如何实现呢？
数据可视化与手动分类的局限性 📊
以下是我们在课程和练习中会反复查看的散点图。图中绘制了不同通道上第一主成分的所有配对关系，例如通道1对通道2、通道1对通道3等。由于有四个通道，因此总共有六张图。
当我们绘制这些第一主成分时，可以清晰地看到数据中的一些结构。例如，我用红色圈出的这个点簇在第3通道上振幅很高，并且与其余数据明显分离。可能还有另一个簇，在观察这些二维边缘分布时也能被较好地分离出来。
那么，我们如何找出这些簇呢？一些电生理记录设备公司提供的软件中，常包含手动聚类切割功能。这相当于我刚刚在散点图中所做的操作：简单地画一个椭圆，并声称“这是一个簇”。然而，由于特征空间通常是高维的，人类很难在低维投影上完成这项工作。
以下是一个手动分类软件的示例。通过圈出所有这些黄色的尖峰，可以识别出一个神经元。
然而，这种策略的效果如何呢？已有研究对此进行了验证，结果表明这并非一个好方法。它主观性强、错误率相当高、人们划定的边界往往不是最优的、极其耗时，并且在不同实验者之间不可重复，也非基于模型。
研究中创建了一些合成数据来验证。在五个模拟通道上，评分者平均发现的单元数量多于实际存在的数量。此外，他们还测量了假阳性和假阴性率，结果明显不理想，且个体间差异很大。
如今，考虑到硅探针或神经像素探针拥有成百上千个通道，这种方法变得完全不可行。因此，我们需要采用算法化和自动化的方法。
另一篇论文也阐明了这一点。在图中，人类操作者有着相当高的假阳性和假阴性率。而如果使用一定数量的主成分进行自动聚类，则可以做得相当好。这是一个拥有真实数据（ground truth）的例子，我们稍后会详细讨论。
自动聚类方法 🚀
相比之下，自动聚类方法可以更客观、基于模型、且可重复。只要公布分析所用的代码和随机种子，任何人都能在同一数据集上运行并获得相同结果。同时，它也便于量化，可以直截了当地计算一些质量指标。
然而，聚类本身的问题在于，无论使用何种方法，它都是一个无监督学习问题。虽然针对该问题有许多解决方案，但也存在许多悬而未决的难题。例如，如何选择聚类数量就是一个重大的未解之谜，目前尚无原则性的解决方案。
我们将使用诸如K-Means（你可能已经了解）或高斯混合模型等方法。如果你对这些基础算法完全不熟悉，我建议你回顾冬季学期的机器学习课程，观看关于高斯混合模型和K-Means聚类的讲座视频，链接已提供。
接下来，我将简要回顾高斯混合模型最重要的部分。
高斯混合模型简介 📈
高斯混合模型将数据点的密度描述为多个高斯分量的加权和。
模型的参数包括：混合比例、高斯分量的均值及其协方差矩阵。例如，下图中的数据是由二维空间中的三个高斯分布生成的，它们具有三个不同的均值和相等的混合比例，其中一个的协方差矩阵比其他的大。
我们如何找到这样一个高斯混合模型的参数呢？通常，我们尝试优化模型下数据的对数似然。将其写出来，形式如下：
L(θ) = Σ_n log( Σ_k π_k N(x_n | μ_k, Σ_k) )
其中，π_k 是混合系数，μ_k 是簇的均值，Σ_k 是簇的协方差。
由于对数内部存在求和，我们无法得到参数的封闭形式解。因此，我们需要使用一种称为期望最大化的算法。
期望最大化算法 🔄
期望最大化算法可用于寻找任何潜变量模型的最大似然估计。潜变量模型是指，并非所有指定模型所需的变量都被观测到。在我们的例子中，潜变量就是簇分配，即我们不知道每个数据点具体来自哪个高斯分量。
在EM算法中，我们有一些模型参数（混合系数、均值、协方差）和未观测到的潜变量（数据点到混合分量的分配），以及观测到的数据（尖峰特征）。
EM算法交替执行以下两步，直到算法收敛（例如似然不再变化）：
- E步：在给定当前参数估计的情况下，估计潜变量。对于尖峰分类，即基于当前簇参数估计簇分配。
- M步：在给定当前潜变量（簇分配）估计的情况下，更新参数以最大化似然。
让我们详细看看这两步。
在E步中，我们使用当前参数值评估责任值。即计算数据点 n 属于簇 k 的后验概率。公式如下：
γ_{nk} = π_k N(x_n | μ_k, Σ_k) / Σ_j π_j N(x_n | μ_j, Σ_j)
例如，在这个小例子中，我们可视化了特定位置上的蓝色和红色簇。在第一个E步中，我们计算这些责任值，并根据责任值为点上色。靠近蓝色圆圈的数据点被染成蓝色，意味着它们有更高的责任值（后验概率）被分配到蓝色簇；反之，远处的红点则更可能被分配到红色簇。
在M步中，我们更新均值、协方差矩阵和混合系数。公式实际上相当简单：
- 均值是数据点的加权平均，权重就是责任值。
- 协方差矩阵同理。
- 混合系数是每个簇的责任值之和除以总数据点数。
然后，我们带着新的参数回到E步，再次计算责任值，如此迭代进行。
下图展示了这个过程：首先初始化模型，然后在第一步的M步中更新参数，移动椭圆到对应数据点的均值位置，接着计算新的责任值，并不断迭代直到收敛。
确定聚类数量与模型选择 🎯
如前所述，如果聚类数量未知，如何确定它尚无普遍接受的解决方案。观察这个数据集，你可能会说需要在这里放一个高斯，那里放一个，下面这里放一个或两个。
一个相对有原则的方法是优化贝叶斯信息准则。为什么不能只用对数似然呢？因为如果我们不断增加簇的数量，负对数似然会一直下降。理论上，如果在每个数据点上放一个单独的高斯，我们能得到最优的负对数似然，但这完全不能揭示数据的真实底层结构。
因此，我们必须平衡模型对数据的拟合程度和模型复杂度。BIC正是通过惩罚参数数量来实现这一点。在BIC公式中：
BIC = -2 * log(L) + P * log(N)
第一项对应对数似然（随聚类数增加而下降），第二项惩罚模型参数数量 P（随聚类数增加而上升）。BIC平衡了这两种趋势，给出一个折中的最优解。当然，这并非唯一方法，在特定假设下可证明其最优性，但对于有限的具体数据，它本质上也是一种启发式方法。你可以选择最小化BIC的聚类数，或者选择BIC不再显著下降的点。最终，你需要通过生物学知识或领域知识来说服自己，聚类结果是合理的。
处理局部最优解 🧩
该算法的另一个问题是可能存在局部最优解。这很容易理解，例如，你可以称这里的簇为簇1，那里的为簇2，但交换它们同样好。除了这种平凡情况，还存在非平凡的局部最优。例如，在这个数据集上多次运行聚类，可能有一个解收敛到看起来像这样的状态，而真实解更可能是那样。
我们如何确保找到尽可能好的解呢？
一种方法是简单地用不同的初始化进行多次独立运行，然后根据BIC挑选最好的结果。这概念上非常简单：从不同的起点出发，优化似然函数可能会陷入不同的局部最优，我们只需选择最好的那个。
另一种算法称为分裂与合并，它也是一种试图跳出局部最优的启发式方法。其工作原理如下：在EM算法收敛后，你可以测试在算法的每次迭代中，是否可以通过分裂一个候选簇（用两个簇替换它，并固定其他分量，运行部分EM）或合并一对候选簇（用一个簇替换两个簇）来进一步降低BIC。通过这种人为干预模型优化的方式，有时可以将模型从局部最优中“推”出来。
处理异常值与非高斯分布 📉
我们在记录中有时遇到的另一个问题是异常值。实际上，我们的数据并非严格高斯分布，它们可能具有重尾。例如，绘制数据点到均值的马氏距离（一种考虑了数据协方差结构的距离）与其概率的关系图。如果是高斯分布，你会期望一条直线。但在许多记录中，我们发现数据点远离均值的概率高于偶然预期。
这表明数据中存在异常值，或者说数据分布具有重尾。这在下图中也有显示：绘制残差直方图或观察混合模型对尖峰数据的拟合情况，可以看到对于这些较高的电压值，其出现概率高于高斯模型的预测。
处理这个问题的一种方法是将高斯混合分量替换为t分布。t分布是高斯分布的推广，它多了一个称为自由度的参数。当自由度趋于无穷时，t分布就是高斯分布。t分布可以被视为具有相同均值的多个高斯分布的混合，这使其能够建模具有重尾或异常值的数据。有论文表明，许多尖峰数据集用t分布拟合比用高斯分布更好。如果你想实现t分布混合模型，也可以在编程练习中尝试。
处理波形漂移 🕰️
我们在尖峰分类数据集中有时遇到的另一个问题是波形漂移。观察所有这些尖峰在特征空间中的位置，可以看到两个清晰的簇。但如果我们观察这些尖峰被采集的时间顺序，会发现其中一个相当稳定，而另一个呈拉长的香蕉形。这显然是由于电极和神经元之间发生了相对移动，例如电极向某个方向漂移，导致该神经元在某一通道上的能量随时间越来越小。
我们也可以将这种漂移纳入模型。一种方法是将高斯混合模型扩展为卡尔曼滤波器混合模型。我们采用标准高斯混合模型，但让簇均值随时间漂移。在每个离散时间点，均值等于前一时间点的均值加上一个漂移项。然后，我们可以在M步中使用这个漂移的均值，并利用标准的前向-后向卡尔曼滤波递归进行更新。如果你不了解卡尔曼滤波的细节也没关系，这里只想说明，你可以通过这种方式扩展模型以解决特定问题。
这种方法效果相当好。如图所示，高斯混合模型的解决方案在这里建模了一个漂亮的分量和一个巨大的分量；而卡尔曼滤波混合模型对于这个红色神经元，几乎将所有高斯分量重叠放置，对于黑色神经元，则允许它们随时间漂移。
最后，我们可以将所有方法结合起来，开发一个漂移的t分布混合模型。如果你对该方法的细节感兴趣，可以查阅相关论文或参考我同事Alexander Ecker的GitHub页面上的实现。
处理同步尖峰与生成模型方法 ⚡
我想讨论的最后一个问题是同步尖峰。大多数情况下，尖峰是相当短暂的事件，彼此孤立发生。但是，特别是当电极附近有很多神经元时，可能会出现尖峰重叠，即多个神经元的多个尖峰在你用作特征的时间窗口内重叠。
为了说明这可能产生的影响，假设一个神经元的尖峰簇在这里，另一个的在下边。那么，在很小一部分尖峰中，你会得到这些混合信号。你甚至会观察到这些混合信号，如果它们频繁发生，甚至可能被误认为是位于中间的第三个神经元。你还可以有各种时间偏移的叠加，从而产生这些模糊的异常值。
如果你想处理这些重叠尖峰，可以看到例如这里有尖峰，那里也有尖峰，你甚至可以通过这两个神经元在特征空间中的不同抖动来预测位置。当然，如果使用多通道探针，这个问题会更加突出，因为在某些通道上发现至少部分重叠的尖峰的概率大大增加。
因此，与我们在讲座主要部分讨论的聚类方法不同，我们可以采用一种相当不同的尖峰分类方法来处理这个问题（也适用于海量多通道数据）。与其使用聚类方法或先提取特征再聚类的方法，我们可以转向使用数据的生成模型来最终找到尖峰序列。
在这种情况下，我们对随时间变化的电压信号进行建模，将其视为尖峰的叠加。公式如下：
V(t) = Σ_n Σ_i w_i * x_n(t - τ_i) + ε(t)
其中，w_i 是尖峰波形模板，x_n 是尖峰序列（二进制信号，具有伯努利先验），ε(t) 是噪声。
如果我们将其用作尖峰分类算法的生成模型，那么我们就可以处理重叠尖峰了。如何优化这样的模型呢？这实际上并不容易，原因在于尖峰是发生在离散时间点上的二进制信号。
我们通过交替优化来解决：
- 假设尖峰时间已知，估计波形（可通过线性回归完成）。
- 假设波形已知且固定，通过模板匹配确定尖峰时间。我们使用一种称为二元追踪的算法来实现。
二元追踪算法在这个生成模型下，我们有一个同时是尖峰序列和波形函数的损失函数。我们希望尽可能好地建模观测到的电压轨迹。我们可以解析地计算在给定时间仓中添加或删除一个尖峰所引起的损失函数变化。因此，可以实现一种贪婪算法：查看所有可能的更改，翻转能带来最大损失下降的尖峰，然后更新损失函数，并迭代进行，同时估计尖峰波形。
算法的初始化也很重要，我们可以采用与讲座第一部分类似的方法。然后，重复以下步骤：给定电压数据和假设的尖峰时间，通过线性回归估计波形；计算残差和每个神经元的先验概率；接着使用二元追踪贪婪算法再次估计尖峰时间。
这可以作为我们第一部分介绍的高斯混合、t分布混合或漂移t分布混合方法之上的一个增强步骤。如果这样做，结果会相当不错。下图是一个来自32通道记录的示例，使用该算法可以识别出分布在多个通道上的相当数量的神经元。放大后，可以查看发现的不同神经元的特征。例如，这个黄色神经元确实有更高的特征，但许多其他神经元也清晰地从背景噪声中显现出来。
本节课中，我们一起学习了尖峰分类的核心概念，从手动方法的局限性，到自动聚类算法（特别是高斯混合模型与EM算法），再到确定聚类数量、处理局部最优、异常值、波形漂移以及同步尖峰等挑战。我们还简要介绍了基于生成模型的先进方法。
在下节课中，我们将探讨如何确保从尖峰分类算法中得到的结果具有生物学意义，并且确实对应于真实的神经元。祝你本周愉快！
003：单神经元识别
在本节课中，我们将学习如何从尖峰数据中识别出单个神经元。我们将重点介绍如何验证自动聚类步骤的生物学意义，确保我们找到的簇确实对应着单个神经元的活动。
概述
上一节我们介绍了特征提取和聚类。本节中，我们来看看如何验证聚类结果，区分哪些簇代表单个神经元（单单元簇），哪些代表多个神经元（多单元簇）。我们将引入不应期和交叉相关图等生物学概念作为验证工具。
簇的初步观察
首先，我们回顾一下聚类后的结果。下图展示了在三个通道上，沿各自第一主成分的二维投影。可以清晰地看到，蓝色簇与其他簇分离得很好，而橙色和粉色簇彼此之间则存在重叠。
单单元与多单元簇
并非所有簇都必然代表单个神经元的尖峰。理解这一点需要回顾微电极记录的工作原理。
- 单单元簇：距离电极尖端几十微米范围内的神经元会产生非常大的信号，通常可以被分离成独立的簇。
- 多单元簇：在更大半径范围内的神经元，由于其信号衰减较多，波形相似，常常无法彼此区分，从而被归入同一个簇，这代表来自特定位置的一组神经元的活动。
利用不应期进行验证
为了评估一个簇是单单元还是多单元，我们需要引入尖峰序列的时间结构知识，主要是寻找不应期。如果某个簇中的所有尖峰都来自同一个神经元，那么每个尖峰之后必须有一段特定的静默期，即不应期。
动作电位发生后，钠离子通道会进入失活状态，无法立即再次开放，这构成了不应期的生理基础。
自相关图与交叉相关图
检测不应期的一个工具是自相关图。对于交叉相关图，则是计算两个尖峰序列之间的时间差直方图。
具体方法是：以一个参考尖峰序列中的每个尖峰为时间零点，统计另一个比较序列中所有尖峰相对于该零点的出现时间，并累加到直方图中。如果是对同一个序列自身进行计算，得到的就是自相关图。
在自相关图中，如果在零时间点附近（例如几毫秒内）出现一个明显的缺口，就表明存在不应期，这是单单元簇的有力证据。交叉相关图则可用于检测神经元之间的耦合关系。
相关函数可以表示为两个时间序列的卷积：
C(τ) = ∫ x(t) y(t+τ) dt
对于离散的尖峰序列，可以简化为所有尖峰时间两两之差的直方图。
分析实例
让我们观察示例数据中几个簇的相关图。
以下是关键观察结果：
- 红色簇：其自相关图相对平坦，在零时间点附近没有显示明显的缺口，表明没有观察到不应期。因此，它很可能是一个多单元簇。
- 蓝色、橙色、粉色簇：它们的自相关图在零时间点附近都显示出一个明显的缺口（凹陷），这为不应期提供了证据。因此，它们都是单单元簇的候选。
此外，观察橙色簇和粉色簇的交叉相关图会发现，它竟然也显示出一个不应期缺口。这非常引人注目，因为它表明橙色和粉色簇中的尖峰很可能来自同一个神经元。结合它们在特征空间中的重叠，这可能是由于电极漂移或尖峰幅度适应等因素，导致同一个神经元的波形随时间发生轻微变化，从而被分到了两个簇中。交叉相关图略微不对称（橙色簇的尖峰倾向于在粉色簇之前），这可能支持尖峰幅度随时间适应的假设。
量化簇分离度：线性判别分析
除了观察相关图，我们还可以量化簇之间的分离程度以辅助判断。一个有用的可视化工具是线性判别分析。
LDA 可以帮助我们找到高维空间中分离两组点的最优方向。其核心思想是找到一个投影方向，使得投影后两类数据的类间方差最大，类内方差最小。在协方差结构各向同性的简单情况下，这个方向正比于两类均值的差向量。在更一般的情况下，需要根据数据的协方差结构进行调整。
我们将对所有簇对计算这个最优的一维投影方向，并将数据投影上去进行可视化。这能直观地显示哪些簇是真正分离良好的。
从图中可以看到，蓝色簇与几乎所有其他簇都分离得很好。而橙色和粉色簇在投影上有相当大的重叠，这与我们之前认为它们可能属于同一神经元的推断一致。
尖峰分选流程总结
至此，我们完成了完整的尖峰分选流程回顾。在前三讲中，我们从原始记录开始，检测尖峰、提取特征、进行聚类，最后验证了哪些簇具有生物学意义，对应着单个神经元。
算法性能评估基准
评估尖峰分选算法的性能是一个挑战，因为获取真实数据（Ground Truth）非常困难。主要有三种策略：
- 真实数据：通过同时进行细胞内和细胞外记录来获取真实数据。这在技术上极具挑战性，数据产出率低，且数据量少可能导致算法过拟合。
- 模拟数据：通过计算机模拟神经组织来生成大量合成数据。其结论的有效性很大程度上取决于模拟的真实性，尤其是噪声特性的模拟非常困难。
- 混合数据：将已知的、来自真实神经元的尖峰波形模板，在已知的时间点注入到真实的记录数据中。这种方法结合了真实噪声背景和可控的尖峰信号，是目前进行大规模算法基准测试的较好方法。
总结
本节课中，我们一起学习了如何识别单神经元。我们了解到，不能仅凭聚类结果就断定找到了单个神经元，必须进行生物学验证。我们引入了不应期的概念，并学会了使用自相关图和交叉相关图来区分单单元簇和多单元簇。我们还介绍了使用线性判别分析来可视化簇的分离程度。最后，我们探讨了评估尖峰分选算法性能的三种基准策略：真实数据、模拟数据和混合数据。掌握这些验证和评估方法，对于确保神经数据分析结果的可靠性至关重要。
004：钙成像数据中的脉冲推断 🧠
在本节课中，我们将学习如何对双光子钙成像数据进行预处理，特别是如何从钙信号中推断出神经元的脉冲序列。我们将从钙成像的基本原理开始，逐步介绍从简单到复杂的脉冲推断模型，并比较不同方法的优劣。
上一节我们介绍了电生理记录中的脉冲排序。本节中，我们来看看其互补技术——双光子钙成像。
双光子钙成像简介
双光子钙成像是一种广泛使用的技术，用于记录神经元群体的活动。其原理是将钙指示剂引入细胞。当神经元放电产生动作电位时，细胞内的钙离子浓度会瞬时升高。钙指示剂与钙离子结合后，其荧光特性会发生改变。通过双光子显微镜测量这种荧光变化，就能间接反映神经元的电活动。
与电生理记录相比，钙成像有其优势和局限：
- 优势：可以记录不常放电的“沉默”神经元；能获取细胞的空间布局信息；结合遗传工具，可以特异性地标记和研究特定类型的神经元。
- 局限：时间分辨率较低（通常为几赫兹到几十赫兹），无法像电生理那样精确定位毫秒级的单个脉冲；钙信号动力学比电压信号慢得多。
从钙信号推断脉冲：问题定义
钙成像数据处理通常面临两个核心任务：
- 从图像中识别细胞体或感兴趣区域。
- 从记录的荧光信号中推断出原始的脉冲序列。
本节课我们聚焦于第二个任务：脉冲推断。其目标是将缓慢变化的钙信号（代理变量）转换回对神经元放电活动（脉冲序列）的估计。
方法一：简单的反卷积模型
我们从一个非常简化的模型开始。假设每个脉冲都会引起一个钙瞬变，而钙瞬变的形状可以用一个指数衰减的核函数来建模。
公式：C(t) = (S * k)(t) + ε
C(t)：观测到的钙信号。S(t)：脉冲序列（0或1的序列）。k(t)：钙核函数（例如，衰减时间常数约500毫秒的指数函数）。*：卷积操作。ε：噪声。
脉冲推断的目标是反演这个前向模型。一个直观的想法是直接对钙信号进行反卷积。
算法思路：
- 对钙信号进行低通滤波，去除高频噪声。
- 使用一个自定义的迭代局部平均程序，去除低于特定阈值的噪声，避免产生伪脉冲。
- 用指数核函数的逆核对信号进行反卷积，得到脉冲序列的估计。
局限性：
- 需要预先知道或校准钙核函数的时间常数。
- 只能重建低通滤波后的放电率，时间精度有限。
- 可能产生无意义的负放电率。
- 在真实噪声数据上，效果往往不理想。
方法二：生成式模型与凸优化
为了获得更好的结果，我们可以使用更复杂的生成式模型，并明确地加入对脉冲序列的先验假设。
模型描述：
- 观测模型：荧光信号
F_t是钙浓度C_t的线性函数加上噪声。
F_t = α * C_t + β + noise - 动力学模型：钙浓度的演化遵循一阶自回归过程，并由脉冲
n_t驱动。
C_t = γ * C_{t-1} + n_t
其中，n_t代表在时间区间Δt内发生的脉冲数，我们假设它服从泊松分布，其期望为λΔt。
直接优化这个模型（寻找最可能的脉冲序列）在计算上很困难，因为涉及离散的泊松先验。一个有效的近似方法是松弛约束：允许脉冲序列取连续的正值，并将泊松先验替换为指数先验。这使得优化问题变成凸优化问题，可以通过交替固定参数（如钙动力学参数和脉冲序列）并迭代求解。
优势：
- 强制脉冲率为非负值，更符合生理实际。
- 能更准确地定位脉冲时间。
- 比简单反卷积方法表现更好。
方法三：扩展的生成式模型 (MLspike)
我们可以进一步扩展模型，以处理钙成像数据中常见的复杂情况。
MLspike 算法的改进：
- 漂移基线：允许荧光基线存在缓慢的漂移。
- 可变振幅：允许钙瞬变的振幅随放电频率变化，可以模拟饱和或超线性响应。
该算法使用了一种类似维特比算法的动态规划反向递归方法，来联合推断最优的钙浓度轨迹和脉冲序列。它还需要一套启发式方法来估计模型参数，例如通过分析钙瞬变振幅的直方图来估计单脉冲的响应幅度。
方法四：监督学习（数据驱动）方法
除了基于生物物理知识的生成式模型，另一种思路是使用监督学习。我们可以利用那些同时记录了钙信号和膜片钳电生理信号（“金标准”）的数据集，训练一个模型直接从钙信号的特征回归到脉冲序列。
例如，可以构建一个小型神经网络，其输出参数化一个泊松分布的放电率，该放电率由钙信号的历史窗口决定。
代码思路（概念性）：
# 伪代码：一个简单的两层网络用于脉冲率推断
rate_t = neural_net(calcium_signal[t - window : t])
spike_train_estimate = sample_from_poisson(rate_t)
在2016年的一项研究中，这种简单的神经网络方法在多个数据集上取得了当时最先进的性能。
方法比较：Spikefinder 竞赛
为了系统比较不同方法，研究者组织了 Spikefinder 竞赛。竞赛收集了超过50个参赛算法，包括各种生成式模型和深度学习模型。
主要发现：
- 在25赫兹的采样频率下，推断放电率与真实放电率之间的相关系数普遍不高（约0.5），这反映了从钙信号推断脉冲的内在难度。
- 表现最好的算法组中，同时包含生成式模型（如MLspike）和深度学习模型。它们的性能在统计上无法区分。
- 不同顶级算法预测结果之间的相关性很高，表明它们学到了相似的数据结构。
- 一些深度学习模型存在对训练数据过拟合的问题。
竞赛结果表明，没有一种方法（生成式 vs. 数据驱动）具有绝对优势。未来结合两者优势的混合方法可能带来更好的结果。
预处理的前一步：运动校正与细胞分割
在能够进行脉冲推断之前，我们首先需要从原始图像数据中提取每个神经元的荧光信号。这涉及两个关键步骤：
- 运动校正：在活体成像中，动物运动或组织搏动会导致图像帧之间的错位。需要使用非刚性运动校正算法。常用方法是基于傅里叶变换的相位相关法，它在局部图像块上计算位移，并能处理不同空间区域的不同形变。
- 细胞分割（ROI提取）：识别图像中对应单个神经元胞体的区域。主流方法主要基于生成式建模思路，将电影数据分解为空间足迹（细胞形状）、时间模板（钙响应动力学）和活动序列（脉冲）的乘积。另一种方法是使用计算机视觉的深度学习模型直接进行细胞分割，但目前应用较少。
总结
本节课中我们一起学习了双光子钙成像数据的预处理，重点探讨了从钙信号中推断神经元脉冲序列的多种方法：
- 我们从简单的线性反卷积模型入手，了解了其基本思想和局限。
- 接着介绍了更复杂的生成式模型，它通过明确的概率模型和凸优化来估计脉冲，效果更好。
- 然后看到了MLspike等算法如何通过扩展模型来处理基线漂移和响应非线性等实际问题。
- 我们还探讨了监督学习方法的潜力，它直接从数据中学习映射关系。
- 通过Spikefinder竞赛，我们了解到目前生成式模型与数据驱动模型性能相当，未来二者结合是发展方向。
- 最后，我们简要回顾了脉冲推断前的必要步骤：运动校正和细胞分割。
理解这些预处理步骤对于正确解读钙成像数据、并最终理解神经环路的功能至关重要。
005：分析脉冲序列 🔬
在本节课中，我们将学习如何分析单个神经元的脉冲序列数据。我们将从可视化脉冲序列开始，介绍时间直方图，探讨如何识别简单细胞与复杂细胞，并学习如何通过调谐曲线来总结神经活动，包括方向与朝向调谐。最后，我们将为调谐的显著性建立一个统计检验方法。
脉冲序列的可视化 📊
上一节我们介绍了神经数据的预处理。本节中，我们来看看如何可视化处理后的脉冲序列数据。
分析数据时，首要步骤是观察原始数据。在试图拟合调谐函数或计算调谐指数之前，应先可视化脉冲序列本身。数据应能直观地展示你感兴趣的效果。
我们将基于前三次讲座和练习中分析过的电生理记录数据。在这些记录中，动物观看了一个移动的正弦光栅刺激，该刺激在多个试次中呈现，每个试次约2秒，光栅在16个运动方向上移动。
以下是可视化脉冲序列数据的第一种也是最简单的方法：
- 排序与对齐：根据每个试次中呈现的刺激方向对记录进行排序，并以刺激开始时间为基准对齐脉冲序列。
- 观察模式：通过这种对齐方式，可以直观地看出某个神经元并非对所有运动方向都有反应。例如，某个神经元可能只对22度、45度、202度和225度这几个方向敏感，且对不同方向的偏好程度可能不同。我们还能观察到响应中的波动模式。
时间直方图与脉冲密度函数 📈
上一节我们介绍了如何通过排序和对齐来观察脉冲序列。本节中，我们来看看如何将其转换为更定量的度量——时间直方图。
我们可以将脉冲序列转换为时间直方图。其计算步骤如下：
- 再次以刺激开始时间为基准对齐所有脉冲序列。
- 将时间划分为合理的小区间。
- 统计每个区间内的脉冲数量。
如果需要，可以进一步平滑这个直方图，得到平滑的曲线，这被称为脉冲密度函数。
这两种方法都有效，适用于不同场景，且密切相关。以下是计算PSTH的示意图：
时间轴被划分为多个区间（bin）。
每个区间内的脉冲数被计数并绘制成柱状图。
平滑后得到连续的脉冲密度函数曲线。
在实践中，选择合适的分区间隔（带宽）通常需要反复尝试，因为PSTH主要是一种视觉探索工具。然而，也存在一种优化策略，通过最小化一个成本函数来选择能最佳估计平均发放率的区间大小。对于脉冲密度函数，同样可以通过优化核宽来获得最佳平滑效果。
发放潜伏期与变点检测 ⏱️
上一节我们介绍了如何计算和优化PSTH。本节中，我们来探讨神经响应中的一个重要时间特征——发放潜伏期。
发放潜伏期是指从刺激呈现到神经元开始显著发放脉冲之间的时间间隔。从理论上定义很简单，但从数据中检测却非易事。
有多种方法可以估计潜伏期：
- 半高法：计算发放率达到其峰值一半的时间点。
- 阈值法：估计神经元的自发发放率，并设定一个阈值。
- 变点检测：这是一个统计问题，可以形式化为寻找发放率从λ1变为λ2的时间点θ。虽然理论上清晰，但在脉冲序列数据上实现效果不一定理想。
- 贝叶斯分箱法：该方法在建模PSTH时，可以同时估计出潜伏期及其不确定性。
简单细胞与复杂细胞 🧠
回顾我们之前看到的PSTH，它们显示出振荡模式。这与刺激的结构以及神经元的感受野特性有关。
在视觉皮层，我们通常有两种细胞模型：
- 简单细胞：具有Gabor滤波器状的感受野，对光栅的空间频率、朝向和相位敏感。当移动光栅经过其感受野时，其发放率会随着光栅相位与感受野最优相位的匹配程度而振荡，这与我们观察到的PSTH振荡模式一致。
- 复杂细胞：对光栅的空间频率和朝向敏感，但对相位不敏感。它们会在其偏好朝向或方向出现时持续发放，而不表现出明显的振荡。
我们可以通过计算线性指数来量化一个细胞更接近简单细胞还是复杂细胞模型。该指数比较了神经元响应在刺激频率（F1分量）上的振荡强度与其平均发放率（DC分量）的比值。
调谐曲线建模 🎯
除了观察时间动态，我们还可以总结神经元在整个试次中的平均发放率如何随刺激参数（如运动方向）变化，这称为调谐曲线。
对于简单细胞模型，其朝向调谐曲线通常可以用余弦函数描述。在几何上，我们可以将朝向视为一个二维单位向量，并将调谐建模为刺激向量与神经元偏好向量之间的点积。
余弦调谐曲线公式：
f(θ) = A + B * cos(θ - φ)
其中，φ 是神经元的偏好朝向。
在实际数据中，神经元的调谐曲线通常比余弦曲线更尖锐。我们可以通过引入非线性（如指数函数）来建模，得到冯·米塞斯调谐曲线。
冯·米塞斯调谐曲线公式：
f(θ) = A * exp(κ * cos(θ - φ))
其中，κ 控制调谐的锐度。
许多神经元不仅对朝向敏感，还对方向敏感（即对相反方向的反应不同）。为了建模方向调谐，我们可以在模型中加入第一谐波分量。
我们可以使用泊松模型来拟合这些调谐曲线。假设在给定朝向下的脉冲计数服从泊松分布，其均值由调谐函数决定，然后通过最大似然估计来拟合参数。
调谐显著性的统计检验 📊
上一节我们学习了如何拟合调谐曲线。本节中，我们来看看如何检验一个神经元的调谐是否具有统计显著性。
我们想知道，观察到的调谐模式是否显著不同于在所有方向上均匀发放的零假设。由于没有现成的标准检验方法，我们可以构建一个置换检验。
置换检验步骤：
- 选择检验统计量：例如，使用线性模型中第二傅里叶分量的系数
Q2，它反映了发放率随方向变化的强度。 - 生成零分布：在零假设（无调谐）下，所有试次是等价的。因此，我们可以随机打乱试次的刺激方向标签，生成许多（如1000次）人工数据集。
- 计算人工统计量：对每个打乱标签后的人工数据集，计算相同的检验统计量
Q2。 - 计算p值：比较实际观测到的
Q2值与零分布。p值等于零分布中大于或等于观测值的比例。 - 做出推断：如果p值小于预设阈值（如0.05），则拒绝零假设，认为该神经元的调谐是显著的。
总结 🎓
本节课中，我们一起学习了分析单个神经元脉冲序列的核心方法。
我们首先介绍了如何通过排序和对齐脉冲序列来可视化数据，然后将其转换为时间直方图或脉冲密度函数以量化时间动态。我们探讨了发放潜伏期的概念与检测方法，并区分了视觉皮层中简单细胞与复杂细胞的响应特性。接着，我们学习了如何用调谐曲线（如余弦曲线和冯·米塞斯曲线）来总结神经元对刺激参数（如朝向和方向）的编码特性，并介绍了使用泊松模型进行拟合。最后，为了判断调谐是否真实存在，我们构建了一个基于置换的统计检验来评估其显著性。
这些方法是理解神经元如何编码外部信息的基础。
006：估计神经脉冲序列中的信息 📊
在本节课中，我们将学习如何将信息论应用于神经脉冲序列的分析。我们将探讨信息论的基本概念，特别是熵和互信息，并了解如何从有限的实验数据中准确估计这些量。
信息论简介
信息论是一门关于如何度量信息的数学理论，其起源可追溯至20世纪40至50年代，由克劳德·香农等人奠定基础。最初，它被用于卫星通信或军事侦察等领域，但很快人们发现，它同样适用于分析神经元等通信系统。
在信息论模型中，我们通常有一个发送者通过一个信道向接收者发送信息。在信道中，信息会受到噪声的干扰。发送者发送符号 X，每个符号有其出现的概率 P(X)。在神经科学背景下，X 通常代表输入系统的刺激。神经系统本身被视为信道。在接收端，我们观测到被记录为符号 Y 的脉冲序列，其出现概率为 P(Y)。
熵：信息的度量
熵是信息论中一个核心且重要的度量，本节我们将介绍离散熵。熵的公式如下：
[
H(X) = -\sum_{x} P(x) \log_2 P(x)
]
其中，H(X) 是随机变量 X 的熵。让我们看几个例子。
抛硬币的熵：对于一个参数为 p 的伯努利分布（例如抛硬币），其熵为：
[
H = -p \log_2 p - (1-p) \log_2 (1-p)
]
如果绘制熵随成功概率 p 变化的曲线，可以发现，当硬币是公平的（即 p = 0.5）时，熵达到最大值 1 比特。这意味着观察一次公平抛硬币的结果，能带来1比特的信息。如果硬币不公平（p 偏离0.5），熵则会降低，因为你在实验前就能以更高概率猜中结果。
K个等可能结果的熵：例如掷一个K面的骰子，其熵为：
[
H = \log_2 K
]
一般而言，随机变量 X 的熵可以理解为，为了确定其具体结果，平均需要提问的“是/否”问题数量，其值在 H 和 H+1 之间。熵的最小值为0，最大值不超过 log₂(符号数量)。当所有符号等可能出现时，熵达到最大。
为什么这个公式是信息的好度量？早期研究者为信息度量设定了一些理想属性：信息应为非负的；概率为1或0的事件不应提供信息（因为它们是确定的）；度量应是连续且随概率单调递减的；独立事件的信息应具有可加性。满足所有这些属性的函数，正是对数函数，从而导出了熵的公式。
熵的度量可以选择不同的底数。以2为底时，单位是比特；以自然常数e为底时，单位是奈特。这是两种最常用的单位。
条件熵与互信息
上一节我们介绍了熵，本节我们来看看条件熵和互信息。
条件熵 H(X|Y) 表示在已知另一个随机变量 Y 的条件下，X 剩余的不确定性。在神经科学背景下，X 代表刺激，Y 代表观测到的脉冲序列。由于信道存在噪声和随机性，一个给定的脉冲序列可能对应多个可能的刺激。条件熵度量的就是已知脉冲序列后，刺激仍然存在的不确定性。
这个概念是对称的。我们也可以定义 H(Y|X)，即给定特定刺激时，脉冲序列的熵。
互信息 I(X;Y) 衡量的是两个变量之间共享的信息量。它可以通过熵和条件熵来定义：
[
I(X;Y) = H(X) - H(X|Y)
]
这可以理解为，观察脉冲序列 Y 后，关于刺激 X 的不确定性（熵）减少了多少。互信息也可以对称地定义为：
[
I(X;Y) = H(Y) - H(Y|X)
]
这衡量了知道刺激 X 后，关于将出现哪个脉冲序列 Y 的不确定性减少了多少。
数据处理不等式
信息论中有一个非常基本的不等式，称为数据处理不等式。它指出，如果我们对随机变量 X 或 Y 进行任何处理（例如变换、分箱、忽略某些特征），那么它们之间的互信息只会减少或保持不变，绝不会增加。
在神经科学中，这意味著许多数据处理步骤都会导致信息损失。例如，将精确的脉冲发放时间进行分箱处理，或者忽略脉冲序列之间的相关性，都会降低我们从数据中计算出的互信息估计值。因此，所有形式的数据误处理或模型误设，都必然会导致互信息估计值的降低。
信息论在神经科学中的应用与挑战
信息论在神经科学中曾受《Spikes: exploring the neural code》一书影响深远，该书为分析神经编码提供了一个理想的定量观察者框架，尤其在研究苍蝇神经元方面应用广泛。
然而，信息论如今并非分析神经数据的首选框架。原因主要有两点批评：首先，信息是一个纯粹的句法概念，生物体不一定实际“使用”这些信息；其次，尽管有其行为相关性，但从数据中准确估计信息量实际上非常困难，这也是我们接下来要重点探讨的问题。
从数据估计熵的挑战
现在，我们来看看从实际的脉冲序列数据中估计熵时会遇到的具体困难。
我们可以将脉冲序列描述为短离散序列，例如在10个时间仓内，每个仓要么有脉冲（1），要么没有（0）。这样一个“词”就是10维的二进制序列。这个离散空间共有 2¹⁰ = 1024 种可能的状态。
在所有可能的概率分布中，为这1024种状态中的每一种分配相等的概率 1/1024，会得到最大的熵。这个最大熵很容易计算：H_max = log₂(1024) = 10 比特。
现在，假设我们从一个真实分布（例如这个均匀分布）中观测到一组数据样本 {x₁, ..., x_N}。一个最直接的熵估计方法是使用最大似然估计：我们简单地统计每个不同的“词”在数据中出现的频率，用这个经验频率除以总样本数来估计概率 P(x)，然后将这些估计值代入熵的公式。
这种“即插即用”的估计量效果如何？核心问题在于未观测到的符号。由于我们只观测到有限数量的样本，总会有一些本应等概率出现的脉冲序列模式在数据集中一次也没出现。这会导致我们对熵的估计产生偏差。
估计量的偏差与方差
要理解这个问题，我们需要回顾估计量的两个基本属性：偏差和方差。
- 偏差：估计量的期望值与真实参数值之间的差异。
- 方差：估计量自身的波动程度。
估计量的总误差可以分解为偏差的平方加上方差，这被称为偏差-方差权衡。
对于熵的估计量，尤其是最大似然估计量，偏差是一个严重问题。模拟显示，即使样本数量达到数千，估计出的熵值仍然系统地低于真实熵值（向下偏差）。只有当样本数量非常大（例如接近或超过可能状态数）时，估计值才接近真实值。
为什么这在神经科学中是个问题？因为互信息 I(X;Y) = H(Y) - H(Y|X) 正是我们关心的量。其中，H(Y) 可以利用所有实验条件下的数据来估计，偏差相对较小。但条件熵 H(Y|X) 的估计则只能使用对特定刺激做出响应的那些脉冲序列数据，数据量少得多，因此其向下偏差会非常大。最终，这会导致互信息 I(X;Y) 的估计产生向上偏差，即我们可能高估了脉冲序列实际携带的关于刺激的信息量，从而得出错误的科学结论。
校正互信息偏差的技术
多年来，人们提出了大量技术来校正互信息估计中的这种偏差。以下介绍几种常见方法。
Miller-Madow 校正：这是一种简单的校正方法。通过对偏差进行泰勒展开，可以证明一阶偏差近似为 （使用的符号数 - 1） / （2N ln(2)），其中 N 是样本数。因此，我们可以在最大似然估计值上直接加上这个项来进行校正。这种方法有一定效果，但校正往往不充分。
刀切法估计量：这是一种重采样方法，可用于估计和减少偏差。其做法是：计算原始数据集的熵估计值，然后依次剔除一个样本，计算剔除后的熵估计值，利用这些值来估计总体偏差并进行校正。这种方法简单且效果显著，能大幅降低偏差。
覆盖调整估计量：这种方法的思路是将熵的公式分解为已观测符号和未观测符号两部分。它同时做两件事：1）缩小已观测符号的概率估计值；2）为罕见（未观测）符号的贡献引入一个膨胀项。这种方法能非常有效地减少平均偏差，但代价是引入了较大的方差，使得单次估计的结果波动很大。
实际上，还存在更多样化的估计量。例如，最佳上界估计量 在平均偏差和方差控制方面都表现优异，仅用约200个样本就能相当准确地估计一个拥有1024种状态的分布的熵，这令人惊讶，但得益于其内置的归纳偏置。
关键在于，从有限数据中估计熵是可行的，但必须谨慎选择适合当前数据和问题特性的估计量。
总结
本节课中，我们一起学习了如何将信息论应用于神经脉冲序列分析。我们回顾了熵、条件熵和互信息的核心概念及其公式。我们深入探讨了从有限神经数据中估计这些信息量时面临的核心挑战——由未观测事件引起的估计偏差，以及这种偏差如何导致互信息的高估。最后，我们介绍了几种校正偏差的技术，包括 Miller-Madow 校正、刀切法和覆盖调整估计量，并强调了根据具体问题选择合适的估计量的重要性。准确估计神经信息是理解神经编码的关键一步，但需要借助严谨的统计方法来实现。
007：单细胞感受野估计 🧠
在本节课中，我们将学习如何估计单个神经元的感受野。感受野描述了刺激的哪些特征会引发神经元的反应。我们将从线性模型开始，逐步深入到更复杂的非线性模型，并学习如何从数据中有效地估计这些模型参数。
感受野与线性模型
上一节我们介绍了调谐曲线的概念，它适用于一维或二维的简单参数化刺激。本节中，我们来看看一个更普遍的问题：刺激或动物行为的哪些方面会导致特定神经元放电？
例如，弱电鱼通过其电感受器官感知周围电场的波动。如果我们记录其侧线叶中某个神经元的活动，并同时记录鱼周围电场的调制信号，我们会发现放电通常发生在电场信号的下降阶段。
为了量化这种关系，我们可以在刺激和神经元电活动之间建立一个线性模型。我们可以尝试将放电率估计为刺激与感受野核的卷积。
公式表示：
rate(t) = ∫ stimulus(t - τ) * kernel(τ) dτ
如果我们计算这个核，它可能呈现特定的形状。例如，当刺激先上升后下降时，神经元会放电。
感受野的理论基础与类型
从理论上讲，这可以追溯到级数展开，类似于函数的泰勒展开，但这里我们将一个函数展开为其他函数的组合。这些核就是我们用来描述刺激与神经元反应之间关系的感受野。
对于一个在空间和时间上的简单细胞，其感受野可能看起来像一个拉长的正形区域和一个拉长的负形区域。随着时间推移，这个形状会变得更加明显，然后逐渐消失。
这种卷积不仅限于一维或二维（如电感受侧线叶的例子），实际上也可以是三维的，即同时跨越空间和时间。我们可以从中看到神经元的潜伏期，即刺激变化与相关放电变化之间的延迟，这大约在75毫秒左右，因为此时感受野的表达最强。
如果我们只沿着一条线切割这些感受野，会看到振幅最高点，并且随着时间的推移，振幅会略微反转。这样的感受野被称为时空可分离感受野。
公式表示：
receptive_field(x, t) = spatial_kernel(x) * temporal_kernel(t)
这意味着我们可以通过一个空间核和一个时间核的简单相乘，来获得二维感受野。因此，二维感受野可以分解为两个一维感受野。
方向选择性与感受野
那么，我们如何分析这种感受野的方向选择性或朝向选择性呢？为此，我们需要能够在时空图中展示移动的光栅。否则很难展示。
在时空图中，x和y是空间坐标，我们让光栅随时间沿某个方向移动。如果我们绘制x和时间T的关系，这种平行于x轴的运动看起来就像对角线条纹。
对于时空可分离的感受野，无论光栅向哪个方向移动，其反应的时间调制看起来都是一样的。这意味着运动方向无关紧要。
那么，如何使感受野具有方向选择性呢？为此，我们实际上需要一个非时空可分离的感受野。在感受野的时空图中，我们需要一个倾斜的轴。
因为如果我们有这样一个非时空可分离的感受野，那么当光栅沿一个方向移动时，对齐度高，反应强；而沿相反方向移动时，无论在任何时间点如何对齐，反应都几乎为零。
结论： 要实现方向选择性，我们需要一个非时空可分离的感受野。
数据分析：线性-非线性泊松模型
前面我们讨论了很多关于刺激与反应、时空可分离与非可分离感受野的关系。现在，我们转向数据分析。
我们通常用一个简单的预测模型来总结这种关系，即所谓的线性-非线性泊松模型。
在LNP模型中，我们假设有一个刺激（具有一维或二维空间维度和一个时间维度）。现代神经元具有感受野，即线性滤波器。将这些线性滤波器在空间和时间上应用于刺激的结果，通过一个非线性函数，然后加上一些噪声过程（在LNP神经元中通常是泊松过程），就生成了神经元的反应。
我们如何估计这样一个神经元的参数呢？我们感兴趣的一个量是所谓的脉冲触发集合。
我们观察那些在脉冲之前出现的刺激。如果这是我们的反应序列，在这些时间点有脉冲，那么我们从刺激中提取在脉冲前不久发生的部分。我们感兴趣的是寻找这个脉冲触发集合中有意义的结构。
脉冲触发平均与感受野估计
我们可以从脉冲触发集合中计算的第一个量，也是感受野的一个简单估计，就是所谓的脉冲触发平均。
具体做法是：取脉冲触发集合的各个组成部分，将它们相加，然后除以脉冲数量，得到脉冲触发平均。
公式表示：
STA = (1/N) * Σ stimulus(t_i - τ)
脉冲触发平均告诉你神经元平均偏好什么。
如果我们的刺激是白噪声刺激（意味着在空间和时间上没有相关性），那么最初通过级数展开定义的线性感受野实际上就是脉冲触发平均。对于白噪声刺激和我们的LNP神经元模型，这一点也成立。
然而，一般来说，情况未必如此，我们稍后会讨论这一点。
需要注意的是，在许多情况下，我们的符号在空间和时间上是连续的，因为空间和时间是连续量。但对于计算机科学家来说，在计算机中精确表示连续信号非常困难，因此我们只需离散化，并将所有积分转化为求和来计算这些量。
最大似然估计与正则化
如前所述，脉冲触发平均是感受野的一个估计，但在某些情况下它是有偏的。例如，如果刺激是非白噪声的（即刺激中存在相关性），比如展示自然电影而不是简单的白噪声时。
在这种情况下，我们可以进行最大似然估计。就像在调谐函数的情况下一样，我们再次假设一个泊松似然，其中给定时间仓内的脉冲按照参数为λ的泊松分布分布。这里的速率λ不是由一维调谐函数指定，而是由感受野与刺激的乘积决定。
然后我们可以写下似然函数，计算负对数似然，并计算导数（对感受野参数求导）。我们可以通过优化似然函数来找到其最大值，从而得到线性滤波器的估计（例如，在LNP神经元模型中，适用于任意刺激）。
这样做的一个问题是，你的估计在平均意义上是无偏的，但方差很高。如果这是真实的感受野，使用一定量的数据，你会得到一个看起来非常嘈杂的最大似然估计。你可以看到整体形状，但它几乎被噪声覆盖。
在统计学中，我们可以正则化我们的估计，或者在贝叶斯框架下计算最大后验估计，将我们对感受野应该是什么样子的先验假设纳入推理框架。
我们可以用这样的图模型来表示：我们的反应Y由刺激X和感受野K决定。我们对感受野K有一些超参数，决定了它的属性。这也被称为先验信息推理。
这里的新内容是带有超参数的先验。然后我们取使后验概率最大化的感受野K，即最大化似然和先验乘积的感受野。
不同的先验与正则化方法
当然，可以使用一系列不同的先验。如果你的先验是系数服从正态分布，那么你的惩罚项类似于L2收缩（也称为岭惩罚）。使用L2收缩，你已经可以得到更好的感受野系数估计，至少方差减少了很多。
你也可以使用拉普拉斯分布作为先验，这相当于使用LASSO惩罚或L1收缩。这会在你的感受野系数中引入稀疏性。在这种情况下，你可以看到这里的许多条目实际上被推到了零。总体而言，感受野看起来比以前更曲折，但它很好地遵循了整体形状。
这些方法很好，但还可以做得更好，使用两种技术：自动局部性检测和自动平滑度确定。
通过自动局部性检测，你可以得到这样的滤波器估计，其假设是感受野在时空上是局部的。因此，在先验中，你以某种方式将感受野在时空中局部化。在先验中，你不再假设系数估计是独立的，而是需要为你的先验分布设置一个协方差矩阵。
你也可以在傅里叶空间中假设局部性，这意味着只允许某些频率分量活跃。你可以在傅里叶谱上设置先验，这实际上效果出奇地好。你还可以结合这两种先验，从而得到非常好的感受野形状估计。
超参数选择与经验贝叶斯
现在，我们如何找到这些模型的超参数呢？到目前为止，我们只是简单地假设了它们。但本演示所基于的论文作者使用了称为经验贝叶斯的技术。
在经验贝叶斯中，我们在这个贝叶斯框架中计算一个称为“证据”的项，然后简单地使用使证据最大化的超参数。对于某些模型，实际上可以以相当封闭的形式或相当高效地计算证据。
由此，我们得到了一整套不同的感受野估计方法，它们都比朴素的最大似然估计要好。我们可以使用L2收缩或岭惩罚、LASSO惩罚，以及各种在空间、频率或空间-频率上局部化的先验。基于一些先验知识，这些方法已经可以局部化我们的感受野估计，并且即使使用很少的数据，也能获得非常好的感受野估计。
样条基函数：一种高效的替代方法
到目前为止，我忽略的一个问题是，为了计算超参数的经验贝叶斯解，我们需要使用线性高斯模型。但当然，我们也希望在LNP模型中能够计算这种数据高效的感受野。
此外，进行这种证据优化步骤通常非常耗时。因此，我们最近开发了一种高效的替代方法，即使用样条基函数。
我们通过选择基函数（这里是样条）来指定我们的先验信念。样条是平滑的函数，或者是在特定准则下最平滑的函数。我们在这里使用的一类样条看起来像这些带有轻微波动的“高斯帽”。
使用这些样条基，我们可以轻松地近似不同形状的感受野，也可以在二维中这样做，因为样条基可以在多个维度上高效计算。
通过这样做，我们获得了什么？如果我们处理高维刺激（因为每个像素都是一个维度，时间也是），我们对所有这些都有独立的系数。使用样条基，我们将刺激空间的维度从三维感受野中的像素数减少到我们需要的样条基分量的数量。
这样做，我们也可以得到非常好看的感受野估计。这是脉冲触发平均，这是使用样条基的最大似然估计器，这是加上L1惩罚和样条基的估计。可以看到，在没有自动局部性检测框架的情况下，我们得到了与贝叶斯框架一样好的估计。
非线性计算：复杂细胞与能量模型
那么，那些在计算上是非线性的细胞呢？在上一节课中，我们简要讨论了复杂细胞和Hubel-Wiesel复杂细胞模型。其思想是，皮层中的简单细胞由LGN神经元的投射产生，这些神经元以有序的方式排列感受野，从而在简单细胞中产生这些条纹状Gabor-like感受野。复杂细胞被认为是具有相同偏好朝向但不同偏好刺激相位的简单细胞的组合，因此复杂细胞整体上对刺激的相位不敏感。
对于这种复杂细胞，线性感受野实际上是平坦的。因此，如果你想模拟它们的反应，需要做一些更聪明的事情。
一般来说，复杂细胞的主导模型是所谓的能量模型。我们有一对滤波器（余弦和正弦滤波器），它们通过二次非线性组合并相加。通过这种二次非线性操作，我们使细胞的反应对刺激的精确位置不敏感，从而引入了一定的位置不变性。
我们可以使用一种技术（例如脉冲触发协方差）来恢复这样的滤波器。在脉冲触发协方差中，我们再次取脉冲触发集合，计算其协方差矩阵，并进行特征向量分析（很像PCA）。然后使用这些作为脉冲触发协方差分量的特征向量。
如果我们有一个根据复杂细胞模型指定的模型神经元（即具有这两个滤波器和相加的二次非线性），那么我们的脉冲触发平均将是无结构的。但如果我们观察特征值谱，会有两个滤波器或两个分量（两个特征值）明显突出。这两个分量很好地恢复了原始滤波器，尽管存在二次非线性。
在下一讲中，我们将更多地讨论这类亚单位模型以及如何拟合它们。
总结 🎯
本节课中，我们一起学习了单细胞感受野估计的核心概念与方法。我们从基础的线性模型和脉冲触发平均出发，探讨了感受野的时空特性与方向选择性。接着，我们深入了更高级的估计技术，包括最大似然估计、贝叶斯正则化（如L1/L2惩罚、自动局部性检测）以及使用样条基函数的高效方法。最后，我们简要介绍了处理非线性反应的复杂细胞能量模型与脉冲触发协方差分析。这些工具使我们能够从神经数据中更准确、更高效地揭示刺激与神经反应之间的复杂关系。
008：基于深度神经网络的系统辨识 🧠
在本节课中，我们将学习如何使用深度神经网络进行神经系统的系统辨识。我们将从简单的线性-非线性模型出发，逐步过渡到更复杂的子单元模型，并最终探讨如何利用深度神经网络来预测神经活动。通过本教程，你将理解这些模型的基本原理、优势以及在实际神经科学数据中的应用。
从线性-非线性模型到子单元模型
上一节我们介绍了感受野模型及其高效估计方法。本节中，我们将探讨更复杂的系统辨识模型，例如使用深度神经网络。
我们上周从一个线性-非线性-泊松模型开始。在该模型中，刺激信号通过一个线性感受野滤波器，然后经过一个非线性变换，并加上噪声来生成响应。
公式：响应 = 非线性(线性滤波器(刺激)) + 噪声
然而，对于执行非线性计算的神经元（例如复杂细胞），这样的模型可能不够充分。因此，我们可以将此概念扩展为所谓的线性-非线性-线性-非线性子单元模型。
在这种模型中，刺激信号通过一组感受野滤波器阵列，每个滤波器都有自己的非线性变换。然后，这些非线性变换的输出被求和，并再次通过一个输出非线性变换。
公式：响应 = 输出非线性( Σ 子单元非线性(子单元滤波器_i(刺激)) ) + 噪声
这类模型为我们模拟刺激-响应关系提供了更丰富的可能性。此外，我们还可以添加类似响应历史滤波器等组件。
子单元模型的优势与应用
子单元模型在复杂性上介于简单的LNP模型和极其复杂的视网膜或皮层回路模型之间。它提供了一个适中的复杂度水平。
在上一讲中，我介绍了使用稀疏基来推断感受野的概念。由于子单元模型有更多参数（每个子单元都需要学习一个参数），在这种情况下使用稀疏基更有帮助。
以下是使用LNP模型和子单元模型对某个神经节细胞进行建模的对比。可以看到，在可用数据量有限的情况下，子单元模型的拟合结果可能较为模糊。但使用稀疏基后，我们不仅能获得更好的LNP模型拟合，还能为各个子单元获得清晰且定位准确的感受野。
在某些系统中，这样的子单元可以提高模型性能。通过绘制预测响应率与真实数据之间的相关系数，我们发现子单元模型通常能改善预测效果。
更有趣的是，有证据表明，在某些情况下，此类子单元模型中的子单元可能对应于实际神经回路中的元素。例如，在视网膜中，神经节细胞从双极细胞接收输入。一些研究表明，通过子单元模型获得的子单元实际上可能对应于单个双极细胞。这有助于弥合统计模型与基于回路的模型之间的差距。
迈向深度神经网络
虽然子单元模型通常受到接近视网膜回路的启发，但我们也可以采取另一种方法：使子单元连续化，并使用深度神经网络来预测神经活动。
在深度神经网络，特别是卷积神经网络中，已经内置了这种子单元机制。卷积神经网络的思想是，你有一个在整幅图像上移动的滤波器，这可以对应于单个双极细胞类型的感受野，在其输出被组合之前应用于整个图像，然后由某种神经节细胞类型从这些表示中采样。
因此，可以将深度神经网络或卷积神经网络的表示视为另一种受视网膜回路启发的尝试，用于预测神经响应。事实上，深度神经网络在这方面表现相当出色。
使用深度神经网络预测神经活动的一个好处是，有大量来自深度学习领域的开源工具可供使用，这大大简化了系统辨识的实现。
代码示例（概念性伪代码）：
# 定义一个简单的卷积神经网络模型
model = Sequential([
Conv2D(filters=32, kernel_size=(3,3), activation='relu', input_shape=(height, width, channels)),
MaxPooling2D(pool_size=(2,2)),
Flatten(),
Dense(units=128, activation='relu'),
Dense(units=1, activation='linear') # 输出神经活动预测
])
model.compile(optimizer='adam', loss='mse')
深度神经网络的训练与扩展
深度神经网络参数优化的关键成分是所谓的反向传播。在前向传播中，我们将一个样本送入带有参数B和W的神经网络，计算线性表示和非线性表示，并将其传递到不同层，直到我们可以评估该特定刺激的损失函数。反向传播则是将损失信号的导数高效地传递回这个计算图的过程。
公式（链式法则）：∂损失/∂参数 = (∂损失/∂层输出) * (∂层输出/∂参数)
借助深度学习的工具箱，我们可以实现更复杂的模型。虽然之前讨论的模型是针对单个神经元拟合的，或者具有为单个神经元训练的特征，但现在我们也可以跨神经元使用特征共享。
我们可以使用从数据中学习但模拟单个神经元响应的特征空间。这在恢复细胞类型时可能很有用。例如，在模拟中，我们有两种类型（小型和大型），在我们的特征共享空间中学习到的核与这些类型非常吻合。
我们还可以为这样的模型添加更多功能。例如，可以添加时间动态。此外，如果记录了眼动或跑步等活动，我们可以通过调制器和移位器网络将其添加到模型中，自动校正并考虑这些因素。
迁移学习与预训练网络
我们不仅可以训练新的网络，还可以使用预训练网络。这意味着我们可以使用在线可用的大型神经网络架构（如AlexNet或VGG16），这些网络已经在海量图像集上训练用于物体分类。
现在，我们可以直接使用这些网络的参数来预测神经活动。其背后的希望是，网络学习到的表示在纯粹的物体分类任务之外，具有更普遍的用途。
例如，在一项研究中，研究人员向两只猴子展示了快速闪现的图像序列，并尝试建模这些图像与引发的神经活动之间的关系。他们使用了一个基于VGG的模型，其中他们简单地使用了预训练VGG的参数来计算特定特征图中的表示，然后训练线性读出器来模拟神经元的活动。在这种情况下，特征提取部分是固定的，他们只是学习了读出部分。
这种方法通常比在特定数据上从头开始训练的卷积网络效果更好，因为它从物体分类任务中继承了许多归纳偏置，因此数据效率更高。
深度神经网络的预测能力与“最兴奋图像”
近年来，研究表明，虽然深度神经网络并不特别接近底层神经组织用于计算信息的神经回路或生物物理机制，但它们提供了非常好的预测模型。一个能证明这些预测模型有多好的概念是“最兴奋图像”。
其过程是：向动物展示一组自然图像，并记录大量神经响应。然后训练一个人工神经网络来模拟这些响应。该网络可以很好地预测响应。接着，可以使用该模型找到一个应该能最优地、最强烈地兴奋某个模型神经元的图像（即“最兴奋图像”）。在后续实验中，可以向真实的神经元展示这个特定的MEI，以检验预测是否正确——即该MEI是否比其他图像更能激发神经元。
这些MEI与简单的感受野不同。在某些情况下，它们可能相似，但在其他情况下，它们可能截然不同。通常，神经元对MEI的响应比对线性感受野等控制刺激的响应要高得多。因此，通过非线性的深度神经网络，可以找到比线性感受野更能激发神经元的图像。
总结
本节课中，我们一起学习了深度神经网络在神经系统辨识中的多种应用。我们从基础的线性-非线性模型出发，探讨了更复杂的子单元模型及其生物学解释。接着，我们深入了解了如何利用深度神经网络，特别是卷积神经网络和预训练模型，来高效地预测神经活动。最后，我们看到了深度神经网络强大的预测能力，例如通过生成“最兴奋图像”来验证模型。总之，深度神经网络已成为神经科学系统辨识中 versatile 且强大的工具。
009：神经群体分析（全）
在本节课中，我们将学习如何分析同时记录的整个神经元群体的活动。我们将扩展之前讨论的单神经元模型，引入神经元自身的放电历史效应以及神经元之间的相互作用，并介绍两种用于描述群体活动模式的统计模型。
从单神经元到神经元群体
在过去的几周里，我们讨论了如何分析单个神经元。这些神经元被单独或作为一组记录，以响应不同的刺激。
今天，我们将进一步扩展我们的模型，讨论如何分析同时记录的整个神经元群体。
扩展模型：纳入放电历史效应
我们目前拥有的模型，主要是具有刺激依赖性的变体，这种依赖性通常被建模为线性滤波器，然后通过非线性传递，最后是一个随机的放电过程。
在过去的讲座中，我们也介绍了这个基本框架的扩展，其中刺激滤波器可以由子单元模型、多个并行的刺激滤波器甚至非线性深度神经网络组成。
然而，我们到目前为止没有包含或没有明确讨论的是单个神经元放电动力学的影响，例如不应期或爆发式放电，以及同时记录的神经元之间的相互作用。我们今天将讨论这些内容。
为了提醒你，不应期是动作电位之后神经元无法再次放电的一段时间。不应期的生物物理基础是钠通道在激活和失活后进入一种失活状态，在这种状态下它们不能被直接再次打开。
为了在我们的神经数据中模拟不应期，我们需要包含一个细胞反馈项。这是一个在放电后兴奋或抑制神经元的项。
因此，我们在这里扩展了我们的速率函数 r(t)，增加了一个放电历史项。我们将放电历史也建模为一个线性滤波器。你可以在下图的示意中看到，神经元自身的随机放电作为后放电滤波器反馈回模型中，与刺激滤波器的效果相加。
这种项的效果在下图中说明。每次放电后，一个抑制性的自我反馈项会降低在任何给定时刻发生放电的概率。
你可以用这个框架模拟非常不同的历史效应。例如，如果我们只输入一个模拟不应期的滤波器，这会微妙地改变我们的放电统计特性。
我们有一个滤波器，可以防止每次放电后极短时间内再次放电。如果你观察对放电间隔分布的影响，你基本上也能看到在开始处概率的下降。
从模型中可以生成或采样的放电序列看起来相当正常。
如果你想模拟爆发式放电，你可以清楚地看到，一旦神经元在这里放电一次，它通常会紧接着再发出几次放电。这可以通过一个这样的爆发滤波器来建模，它在不应期之外，在放电后不久有一个这种自我兴奋的凸起。
这非常强烈，以至于在放电间隔分布中，你可以在正的时间滞后处看到一个明显的凸起。
最后，我们也可以在放电活动中产生振荡。在原始放电数据中，使用振荡性放电历史滤波器，这些振荡有点难以发现。
这些滤波器然后也会产生振荡性的放电间隔分布。
你可以像以前一样，通过优化似然函数，从数据中估计这些滤波器。
如果你有包含多次试验的数据，并且想要同时模拟PSTH（刺激后时间直方图）或刺激依赖性以及放电历史依赖性，那么你需要同时推断这两项。因为在时间 t 发生的放电，可能是由于刺激开始导致刺激诱发速率高，也可能是由于最近的放电历史（例如，在一次爆发中）的结果。如果你在推断模型之前只是跨试验平均，这会混淆这两者。因此，如果你同时推断它们，你至少可以在一定程度上分离这些效应。
我们可以在这里说明这一点。如果我们有这种数据，其中神经元在刺激开始后不久有一段升高的放电期。那么，如果我们不考虑放电历史效应，仅根据数据估计PSTH，我们会得到这种PSTH。但如果我们同时估计这两项，那么这实际上与我们的模型PSTH匹配，并且这里的一些变异性、一些依赖性实际上被归因于自我兴奋项。
超越成对效应：神经元间相互作用
你也可以超越成对效应，模拟同时记录的神经元之间的相互作用。
为此，我们简单地添加一个耦合滤波器。耦合滤波器的工作原理与后放电滤波器类似，但它将一个神经元的放电活动作为输入，进行滤波，然后输入到另一个神经元，反之亦然。
这将允许我们模拟单个神经元之间的依赖性。
例如，在2008年的这篇论文中，他们拟合了这样一个大规模群体模型。他们推断了个体的刺激滤波器（在空间和时间上）、后放电滤波器（你可以清楚地看到不应期），并且他们还研究了耦合滤波器。例如，对于中间的这个ON细胞，你可以查看来自周围ON和OFF细胞的所有耦合滤波器。
对于这个中间的OFF细胞，与附近其他OFF细胞的耦合是正的，而与附近其他ON细胞的耦合是负的。因此，这些细胞显示出相互影响的关联活动。
模型的灵活性与非线性
我们也可以，我之前简要提到过，替换我们的非线性。通常，我们使用指数非线性，因此在统计术语中，我们基本上拥有泊松广义线性模型。但我们也可以使用，例如，S型非线性。如果我们使时间仓非常小，使得每个时间仓最多只有一个放电，那么这就是一个二项式或伯努利广义线性模型。
原则上，我们也可以使我们的非线性非常灵活，并从数据中拟合它，要么使用直方图方法，要么使用样条基。例如，在这里，你可以看到这个非线性由样条基的系数参数化，因此可以非常灵活地模拟模型线性阶段输出与放电活动之间的不同关系。
总结：预测性群体模型
总而言之，我们这里的预测性群体模型非常灵活，我们可以轻松地调整它们。有许多现有的框架可以用来拟合它们。例如，几乎每种编程语言中都有的GLM拟合工具箱。从统计学的角度来看，它们也得到了很好的理解。如果我们有想要添加到模型中的额外项，例如行为数据、眼动数据、奔跑数据，只需以与包含放电历史滤波器相同的方式，将另一个预测因子输入到同一框架中即可。
人们还研究了大量的扩展，以估计隐藏神经元或非线性相互作用，模拟共同输入，模拟状态依赖性等等。
当你将这样的模型拟合到数据时，一个普遍的问题是，最终，耦合滤波器没有一个良好的机制基础。它们模拟了神经元之间的功能或统计连接性，但它们并不是这些功能实现背后真实突触连接性的良好或必要模型。
因此，在文献中，你会经常听到“功能连接性”这个术语。实际上，如果你愿意，它并不是连接性，而更像是功能影响或有向相关性之类的东西。
第二类模型：统计群体模型
今天我想讨论的第二类群体模型是统计群体模型。在这些统计群体模型中，我们的目标不是预测神经元或神经元群体对刺激的反应，而是描述神经反应的概率分布，无论是在自发活动期间还是在响应不同刺激集合时。
因此，我们的目标是对这个分布 P(x1, ..., xN) 进行建模，它模拟了群体在特定时间窗口内的活动。例如，这里有40个神经元，时间在这个维度上。现在我们可以将时间分割成不同持续时间的窗口，并对由此产生的活动模式进行统计建模。
如何建模群体活动分布
那么我们该怎么做呢？如果我们使我们的放电计数窗口足够大，那么我们可以简单地假设我们的放电遵循高斯分布，并用简单的高斯分布来近似放电计数分布，尽管它是离散的。
例如，在这里，如果我们计算整个试验（整个12秒）的放电次数，我们会得到相当大的数字，使得我们的放电计数分布作为高斯分布是一个很好的模型。高斯分布不仅允许我们模拟平均放电计数，还可以模拟它们的协方差或相关性。因此，对于大的时间窗口，均值和协方差完全描述了放电计数分布。
实际上，在这种情况下，使用相关的泊松分布要困难得多，因为泊松分布的离散性。
但我们想要模拟的大多数有趣过程发生在更快的时间尺度上，那时放电计数显然不是高斯的。因此，让我们走到另一个极端，在非常短的时间窗口（比如20毫秒）内模拟放电计数，在那里我们通常发现要么有一个放电，要么根本没有放电。
在这种情况下，它看起来是这样的。在这里，如果我们放大那个子群体在一个小时间窗口上的活动，我们可以看到任何这些指示的时间仓要么包含0个要么包含1个放电，所以我们可以将我们的放电模式转换为如下所示的二进制模式。我们将这些模式记为 B，然后我们问，对于这些二进制模式 B，什么是好的概率分布？当然，它不是高斯分布。高斯分布是二进制分布的一个糟糕近似。那么，我们还能用什么？
我们的第一个起点将是一个独立的伯努利分布。
因此，我们说这里的概率分布，即这些二进制模式 B 的独立概率分布，仅仅是每个具有其边缘放电概率的个体神经元的伯努利分布的乘积。
例如，在这种情况下，对于下面的那个神经元，我们计算它在26个时间仓中放电两次，因此我们估计的放电率为 2/26，并将其用作我们伯努利分布的均值。如果我们这样做，该模型的参数非常容易估计。但它没有捕捉到神经元之间的任何依赖性。
另一方面，在另一个极端，我们可以直接指定完整模型，即给定群体中所有可能状态上的分布。对于下面这10个神经元，我们总共有1024个（即 2^10）不同的状态，每个不同的模式可以被唯一地分配一个状态索引。然后原则上，我们只需要计算每个状态在数据中出现的频率来估计其概率。这显示在上面。然而，由于这种指数关系（10个神经元有 2^10 种状态，20个神经元有 2^20 种状态），该模型的参数很快变得无法估计，因为它捕捉了神经元之间的所有依赖性，但过于复杂而无法实际拟合数据。对于10个神经元，我们有1000种状态；对于20个神经元，我们已经有一百万种不同的状态，其中一些只会非常罕见地出现。由于我们的模型估计总是受到记录时间的限制，这很难处理。
最大熵原理
那么我们还能做什么？我们将诉诸一个称为最大熵原理的原则。你已经从过去的一节课中知道了熵是什么。
最大熵原理背后的想法实际上非常简单。如果我们知道一些统计量（统计量只是对这个二进制向量的一种测量），我们知道活动模式的一些统计量。例如，每个神经元的平均活动、它们的相关性、群体放电计数等等。这些就是这些函数 F_k(B)。然后，我们想回答这个问题：满足这些约束条件，但不对数据的分布施加任何其他统计结构的分布是什么？因此，它只受到必要的约束，但不过度约束。这就是所谓的最大熵分布。
你可以通过最大化熵 H(B) 来找到该分布，但要满足你所拥有的统计量的这些约束条件。并且可以证明，例如，对于一组统计量 F_k(B)，满足这些约束的最大熵分布取简单的指数形式。
然而，这种形式的简单性具有欺骗性，因为难点实际上在于归一化常数 Z 的细节。正如我们稍后将讨论的，这很难计算。这里有一些统计量的例子。如果我们只设 F_k(B) 为 B，这就是均值。如果我们把它设为下面的项，那么我们就有协方差，或者我们也可以有群体放电计数。
成对模型与伊辛模型
在成对模型中，我们希望约束我们的模型以捕捉神经元之间的成对相关性。
因此，如果这些神经元之间没有相关性，而这些神经元之间有一些相关性，我们所做的本质上是使下面这个点（两者同时放电）的概率，与这里这个点相比，变得更大。
好的，那么，如果我们写下满足均值和这些相关性约束的最大熵分布，即二阶最大熵模型（也称为伊辛模型，以德国物理学家命名），它取以下形式。其中模型的参数实际上是 J 和 h。这些需要从数据中估计。你可以将其视为一个二次模型。它有 n*(n+1)/2 个参数。对于一个10维模型，学习这些参数是困难的，正如我所说，因为你需要估计这个归一化常数 Z，以便在任何时候评估概率，因为这不能以封闭形式完成，因为你需要再次对所有可能的状态求和。
你需要通过采样（例如）来近似它。因此，没有封闭形式的极大似然解。你可以做的是使用一种称为玻尔兹曼机学习的技术，在那里你只需评估似然函数关于参数的梯度，然后进行梯度下降。
有趣的是，实际上就是这个公式，它向你展示了在每个步骤中你需要如何改变成对耦合参数 J_jk。你只需评估数据中的相关性，然后减去当前模型在当前参数下预测的相关性。这就是梯度指向的方向，也是你需要迈出一步的方向。
你从模型中得到这些预测，例如，通过蒙特卡洛方法。因此，你需要从模型中采样，例如使用一种称为吉布斯采样的技术。通过这种吉布斯采样，你也可以近似这个归一化常数，如果你需要评估概率的话。
但现在有很多技术可以近似该模型的各个方面，或者使用各种技巧使推断更快。但这是你需要做的基本推断步骤。
模型性能评估
现在，所有这些模型在描述视网膜等神经群体放电活动方面有多好？例如。我们在这里可以看到独立模型的预测，下面你有观察到的模式率，以及模型预测的模式率。正如你所看到的，沿着对角线有相当多的分散点，这里有一个放电的模式，下面有两个放电的模式，等等。因此，对于那个特定群体，独立模型并不是一个特别好的模型。
如果使用伊辛模型来拟合数据，那么所有这些点都更接近对角线。事实上，你可以证明，从对角线观察到的分散程度，几乎与这个数据集原则上所能约束的程度一样接近。
然后，人们还可以评估这个模型解释了总结构的多少。如果我们尝试拟合我们的完整模型和二阶模型，在这种情况下，你可以看到这里几乎解释了数据中超过90%的方差。
替代模型：二值化高斯模型
由于伊辛模型中的采样和模型参数推断并不容易，包括我们在内的人们长期以来也在思考替代方案。
一个非常简单的、可以轻松使用的替代方案，例如用于生成具有指定均值和协方差的二进制数据，就是所谓的二值化高斯模型。在二值化高斯模型中，你假设有一个潜变量，它具有相关性和潜均值。然后你有阈值 θ1、θ2，就像我们LNP模型中的非线性，但实际上是一个阈值非线性，它将这个潜在的连续分布二值化为离散的二进制结果。
因此，如果你从你的潜高斯分布中采样，并且你的点落在下面的这个象限，那么你会分配那个 (0,0) 状态，即两个神经元都不放电。如果它落在上面的那个象限，你会说神经元2放电；如果落在下面这个象限，神经元1放电；如果落在上面这个象限，两者都放电。
因此，通过移动潜高斯的均值（或者等价地，保持均值为0并移动阈值），你可以轻松地改变每个神经元的平均活动。通过改变相关性，你可以改变在不同象限中观察到模式的概率。
在某种意义上，你可以将其视为一个膜电位，一个非常简单的膜电位电压模型，它在背景中以高斯分布波动，并在超过阈值时产生放电。
优点是，该模型的参数可以以非常直接的方式计算。通过使用逆高斯正态分布函数，我们可以简单地从观察到的二进制数据的均值计算出潜高斯变量的阈值均值。我们可以对协方差矩阵的参数做同样的事情。
二值化高斯模型的好处是，对于参数空间的大部分区域，它实际上非常接近最大熵。因此，它几乎没有做出额外的假设。这些假设也可以被很好地理解。
因此，如果我们在这里比较，例如，二值化高斯模型和伊辛模型在参数空间范围内的差异，你可以看到，只要你的相关性不是太大，两个模型的熵实际上非常接近。如果相关性变得很大且均值非常小，那么它们就会产生相当大的分歧。
这种二值化高斯模型的优点是模型更容易拟合，同时保持了接近最大熵的特性。因此，它没有做出超出严格必要范围的结构假设。你可以将其解释为阈值膜电位，并且它也可以推广到不同于0和1的放电计数。
但是存在一些二进制分布，没有对应的二值化高斯模型存在。可以证明这一点。并且也没有简单的方法来计算似然函数，因此你也必须评估高维高斯积分。
课程总结
本节课中，我们一起学习了两种方法来模拟同时记录的神经群体中的群体结构。我向你展示了如何扩展单神经元模型以包含放电历史效应和神经元间耦合，并介绍了最大熵原理下的伊辛模型以及更易处理的二值化高斯模型，用于描述群体活动的统计特性。感谢你本周的关注。
010：神经群体分析（二）
在本节课中，我们将继续学习如何分析神经群体活动。我们将重点探讨神经元反应中的“噪声相关性”，了解其来源、对信息编码的影响，以及如何通过潜变量模型来分析和解释这些相关性。
概述：神经反应的变异性与相关性
上一节我们介绍了分析神经群体活动的基本概念。本节中，我们来看看神经元反应中一个关键特性：变异性与相关性。
即使面对完全相同的刺激条件，神经元的反应也具有高度变异性。这种在不同试次中观察到的变异性，通常被称为“噪声”。这种噪声甚至可以在不同神经元之间相互关联，我们称之为“噪声相关性”。
我们将其与“信号相关性”区分开来。信号相关性源于神经元调谐曲线的相似性。例如，两个具有相同偏好方向的神经元，其平均反应模式会高度相关。
噪声相关性的来源与影响
来源：共同输入
噪声相关性是如何产生的？一个经典理论认为，它们源于神经元接收的共同输入。例如，如果两个神经元从同一组上游神经元接收大量输入，它们的反应波动就会相关。
影响：对群体编码的限制
噪声相关性对从神经群体中读取信息有重要影响。考虑一个场景：我们需要根据一群神经元的活性来区分垂直和水平光栅。
- 我们有两个“决策神经元”，一个汇总所有偏好水平方向的神经元，另一个汇总所有偏好垂直方向的神经元。
- 如果神经元之间没有噪声相关性（相关性为0），那么信噪比会随着群体规模的增加而增长。
- 然而，如果存在正的噪声相关性，信噪比的增长会达到饱和，且相关性越高，饱和得越快。
这意味着，超过一定数量后，增加更多的神经元对决策过程没有帮助。这解释了为什么在某些实验中，单个神经元的敏感度与整个生物体的行为表现惊人地接近——群体汇总信息的潜力可能被神经元之间的相关性严重限制了。
噪声相关性的结构
为了直观理解为什么噪声相关性有害，我们可以可视化群体活动。每个点代表一个神经元，其位置对应其偏好方向，颜色代表不同试次下的活性。
在无相关性的群体中，每次试次的活性分布大致围绕其平均调谐曲线（黑色曲线）波动。
在均匀相关性结构中，所有神经元具有相同的噪声相关性。这导致整个群体活性“山丘”在每次试次中整体上移或下移。这种全局波动如果能够被识别和消除，实际上可能有助于解码。
在有限范围相关性结构中，偏好方向相近的神经元具有更高的噪声相关性。这导致群体活性“山丘”的峰值在每次试次中发生漂移。从这种结构中解码刺激信息会比从独立群体中解码困难得多，即使采用最优读取方式，噪声相关性也会显著损害编码性能。
真实系统中的噪声相关性
在实际系统中，噪声相关性有多大？一项对清醒动物视觉皮层的研究发现，平均噪声相关性非常小（约0.02），即使对于物理位置很近、调谐曲线相似的神经元对也是如此。
这与早期许多研究报告的较高相关性形成了对比。那么，哪些因素会影响我们对相关性的估计呢？
以下是可能影响噪声相关性估计的关键因素：
- 全局脑状态：如睡眠或麻醉。
- 决策与注意：例如，动物注意某个区域时，该区域神经元的增益会被调制。
- 不同行为策略。
- 技术伪迹：如 spike sorting 错误或电极漂移。
案例研究：麻醉对噪声相关性的影响
本节我们聚焦于全局脑状态，特别是麻醉的影响。麻醉会诱发大脑皮层神经活动产生类似睡眠的波动。
比较清醒和麻醉状态下的数据：
- 清醒状态：神经活性波动主要与刺激相位锁定。
- 麻醉状态：神经活性会出现与刺激无关的、全脑范围的同步波动（如暂停发放），这些波动在不同试次中模式不同但在神经元间高度一致。
量化差异：
- Fano因子（方差与均值之比）：麻醉状态下更高（2.4 vs 1.1），表明神经活性波动更大。
- 平均噪声相关性：麻醉状态下更高（0.06 vs 0.01）。
使用潜变量模型分析群体活动
为了理解并量化麻醉状态下的这种变化，我们引入潜变量模型，具体是高斯过程因子分析。
GPFA 假设观测到的脉冲计数由一个未观测到的网络状态（潜变量）生成，该状态随时间平滑变化，并同时影响所有神经元。
模型的核心公式描述了观测活性 y、潜变量 x、加载权重 C 和独立噪声之间的关系。总协方差可以分解为：
Cov(y) = C * C^T + diag(ψ)
其中 C * C^T 代表由潜变量解释的协方差，diag(ψ) 代表各神经元的私有噪声。
通过期望最大化算法，我们可以从数据中估计出潜变量轨迹 x(t) 和模型参数。在麻醉数据中，推断出的潜变量清晰地捕捉到了群体兴奋性的全局波动。
模型应用：解释麻醉下的高相关性
使用训练好的 GPFA 模型，我们可以在测试集上估计潜变量，并计算残差协方差矩阵（即去除潜变量影响后的神经元间相关性）。
结果发现，去除这个代表麻醉诱导的全局脑状态的潜变量后，麻醉状态下的噪声相关性结构变得与清醒状态下的无法区分。这表明，麻醉状态下观察到的更高相关性，很大程度上可以由这个单一的、随时间波动的全局脑状态变量来解释。
该潜变量波动的时间尺度约为250毫秒（1-2 Hz），这与已知的麻醉下脑波动特征一致。
模型扩展与前沿
基础 GPFA 模型可以扩展，例如考虑脉冲计数的泊松特性，或使用更复杂的生成模型。
一个前沿方向是潜因子分析动态系统，它利用循环神经网络作为编码器和生成器，将高维神经活动映射到低维潜空间。这类模型能更好地在单试次水平上解析神经群体动力学，揭示在试次平均中掩盖的细节，对于研究决策策略等具有重要意义。
总结
本节课中我们一起学习了：
- 噪声相关性的概念、来源及其对神经群体信息编码的限制性影响。
- 噪声相关性的不同结构（均匀 vs 有限范围）如何影响解码性能。
- 在真实系统中，全局脑状态（如麻醉）是导致噪声相关性升高的重要因素。
- 如何利用高斯过程因子分析等潜变量模型来推断和分离神经活动中的全局波动成分。
- 通过建模，我们可以将麻醉下增强的相关性归因于一个单一的、缓慢波动的脑状态变量，从而更清晰地理解神经编码的本质。
潜变量模型是研究神经群体动力学的强大工具，通过结合低维状态空间和到神经反应的映射，帮助我们揭示神经群体活动背后的结构和原理。
011：单细胞转录组学 🧬
在本节课中，我们将学习神经科学中的一种新数据模态——单细胞转录组学。我们将了解其基本原理、数据获取方法、常见的数据分析任务，并重点探讨数据预处理、可视化以及细胞类型识别等核心步骤。
什么是单细胞转录组学？🧬
上一节我们介绍了课程的整体安排，本节中我们来看看单细胞转录组学的定义。
细胞中，遗传信息从DNA通过转录过程生成mRNA，mRNA再通过翻译过程生成蛋白质。单细胞转录组学专注于测量单个细胞的mRNA含量。
在许多情况下，这可以告诉我们细胞的身份、细胞类型以及细胞后续将转录或翻译的蛋白质种类，从而揭示大量关于细胞身份的信息。
数据获取方法 🧪
在过去的十年中，出现了多种从单个神经元或身体其他部位单个细胞获取这些测量的方法。这些方法大致遵循以下流程：
以下是典型单细胞RNA测序（scRNA-seq）工作流程的步骤：
- 组织解离：将组织解离，使单个细胞在溶液中真正分离。
- 细胞捕获：捕获单个细胞。
- 裂解与逆转录：裂解细胞膜，通过逆转录将mRNA反转录回cDNA。
- 扩增与准备：扩增细胞的DNA含量，并为测序做准备。
- 测序与分析：使用新一代测序技术进行快速测序，然后进行数据分析。
一个例子是微滴法：细胞从上方进入微流体装置，与条形码结合，并被浸入油滴中。在油滴中发生细胞裂解，使得单个细胞的mRNA内容与一个条形码结合，从而可以确定细胞身份，随后进行扩增和测序，最终我们获得这些细胞中许多不同基因的表达水平。
目前可用的方法范围很广，从相当慢的方法（如显微移液或激光捕获显微切割）到非常快的方法（如微滴或微流体方法）。
显微移液方法的优点是几乎可用于任何组织，并且可以真正靶向感兴趣的特定细胞。例如，如果在大脑的某个部分用荧光标记标记了一个细胞亚群，显微移液可以让你特异性地靶向这些荧光细胞，从而在测序前富集或选择该亚群。当然，另一方面，它相当耗时。
另一个极端是微滴法，如上一张幻灯片所示。这些方法速度超快，也非常便宜。可以以可承受的价格观察大量细胞，如今可达数十万。但另一方面，对观察目标没有同等的控制力。
研究现状与示例 📈
单细胞转录组学是近年来神经科学中非常热门的话题。如果你查看过去20年左右提及该术语的研究数量，可以看到近年来急剧上升。
不仅是研究数量急剧上升，每项研究测序的细胞数量也迅速增加。大约四五年前，拥有几千个细胞的研究相当常见。如今，许多研究报告的细胞数量高达数十万甚至数百万。
例如，该领域一些最突出、最美丽的可视化结果来自斯坦福大学实验室的两项研究。他们观察了小鼠神经系统的分子结构，首先是成年动物，后来也包括发育过程。他们测序了神经系统中许多不同部位的细胞，包括大脑、脊髓甚至消化道的肠神经系统。他们的研究总共约有5万个细胞，然后创建了如此美丽的可视化结果，真正展示了组织中细胞的多样性。
在最近的一项研究中，他们增加了发育维度，观察了这些不同的“岛屿”和“大陆”在动物发育过程中是如何出现的。在这个美丽的可视化中，你可以看到一个大的“神经元大陆”，并且作为该神经元大陆中的“等高线”，你可以看到不同的发育阶段。
数据结构与特点 📊
上一节我们了解了单细胞转录组学的研究概况，本节中我们来看看处理这类数据时遇到的数据结构和常见分析任务。
这与我们目前处理的数据有些不同。在锋电位序列或双光子成像数据中，时间总是一个维度，它真正排序了数据的维度或特征。在数据矩阵中索引接近的特征，在时间上也接近。因此，在时间上存在一种自然的顺序。图像也是如此。例如，计算锋电位序列与图像关系的模型时，也存在固有的空间维度来定义邻域。
对于单细胞转录组学数据，情况不一定如此。在这种情况下，特征维度是基因。在小鼠和人类中，我们大约有25000个基因可以测量表达水平，但这些基因没有自然的顺序，因此数据矩阵的列可以被置换。例如，在该数据矩阵上计算卷积没有太大意义。
数据矩阵的另一个维度通常是细胞。可能从只有1万个细胞开始，如今通常更多。但正如我们将看到的，这是一个非常稀疏的矩阵，没有很多条目。
具体数字取决于情况。如果我们处于低端，因为研究规模较小，但我们测量了25000个基因，那么我们处于样本数量确实小于特征维度的范围。但另一方面，如今如果我们测量更多的细胞，也可以进入更好的采样状态。
核心数据分析任务 🎯
以下是使用该数据矩阵进行的典型数据分析任务，这些任务在初始数据获取之后进行：
- 预处理与质量控制：执行质量控制和一些预处理步骤，以便能够很好地处理数据并去除技术伪影。
- 数据可视化与探索性分析：进行数据可视化以进行探索性分析。
- 细胞类型识别：尝试识别细胞类型，即执行某种聚类。
当然，也可以用转录组学数据做其他事情。例如，可以观察发育轨迹并根据伪时间对细胞排序，可以尝试将转录组学测量与更多生理学启发的模型联系起来。但今天，我们将专注于这些基本任务。
在许多情况下或许多实验范式中，这个数据矩阵以所谓的UMI计数矩阵的形式出现。UMI代表唯一分子标识符。这些是在逆转录后添加到DNA片段上的标签，用于识别输入的DNA分子。这是在扩增之前添加的，它可以用来减少扩增偏差。如果你查阅关于这些实验技术如何工作的技术文献，有相当多的潜在偏差来源需要处理，UMI是解决该问题的一种实验方案。
我们今天在示例中将要查看的数据来自Kenneth Harris及其同事的这篇论文，该论文研究了CA1海马神经元以及单细胞转录组学揭示的多样性。
任务一：数据预处理 🔧
现在让我们看看我们的三个数据分析任务，从第一个开始：预处理。
首先，如果我们有一个数据矩阵，或者我们从Ken Harris的研究中获取数据矩阵，那么我们发现总共有大约3600个细胞和大约18000个不同的基因。这也是数据矩阵的形状告诉你的信息。
在每个条目中，是每个基因在每个细胞中的UMI计数。即在扩增版本中，在给定细胞中检测到给定基因的频率。
整个矩阵的最小值是0。这意味着有些基因在某些细胞中根本没有被检测到。这并不奇怪。最大计数是923。这意味着至少有一个基因在一个细胞中，我们发现了923个重复。
现在，如果我们查看数据矩阵中有多少条目实际上大于0，即在多少细胞和基因组合中我们至少有一个条目，那么这个值实际上相当小，只有15%。所以这是一个非常稀疏的矩阵，对于大多数细胞中的大多数基因，我们没有发现任何表达。
最后两个测量值是我们所说的测序深度。这意味着我们在给定细胞中跨所有基因发现了多少基因计数。具有最小基因计数的细胞的测序深度大约为1111，具有最高基因计数或UMI计数的细胞的测序深度为20000。
你也可以查看这些测序深度的完整分布。在这个直方图中，你可以看到它从大约1000开始，然后也有一个很长的尾巴延伸到非常高的数字。但测序深度的主体平均在5000左右。
但正如你可以想象的，这当然是我们数据分析中一个潜在的混淆因素。因为如果对于一个给定的细胞，我们只能检测到任何基因的1000个拷贝，那么，当然，找到一个与在检测到20000个基因的细胞中表达水平相同的基因的概率将是完全不同的。
理解技术变异：过离散性
让我们更详细地看看这个混淆因素可能是什么。
因为UMI计数是计数，我们首先可以问，这些UMI计数的分布实际上是否符合泊松分布？因为泊松是计数变量的标准或首选分布。如果你查看这些，并将它们与相应的泊松分布可能的样子对应起来，那么你可能会看到，特别是在这些较大的范围内，泊松分布往往变得相当对称，而我们在UMI计数数据中发现的分布实际上有一个相当长的尾巴。
所以可能泊松分布不太适合这个UMI计数数据。
我们可以通过绘制跨基因的均值-方差关系来更精确地量化这一点。这个图中的每个点代表一个基因。我们查看所有细胞的平均表达和所有细胞的表达方差，并将它们相互绘制。
对于一个纯粹的泊松分布，这些点应该都落在一条线上，因为对于泊松分布，均值等于方差。事实上，对于许多基因，我们离得并不太远。
但有两个变异来源加入到这个图中。一个是生物变异，另一个是技术变异。生物变异实际上是我们最终感兴趣的。生物变异来自例如不同细胞类型之间基因表达的差异。所以一种细胞类型可能表达某种受体基因，因为它需要持续地将该受体放在其树突上，而另一种细胞类型不这样做，因此不表达该基因。所以这将在某些基因中导致真正的生物变异，其大于纯粹从计数统计中预期的变异。
另一方面，你也有技术变异，它来自这种测序深度效应，即不同细胞只是被更频繁地测序或扩增。通常，你可以看到生物变异体现在这些散布在线周围的单个点上，而技术变异主要影响高均值处的长尾。
我们可以用另一种方式来看待这一点。通过取方差与均值的比率，这也称为离散系数。离散系数是方差除以均值。如果对于泊松分布，方差和均值相同，那么我们期望离散系数为1。
这基本上带出了上一张幻灯片最后一个图中的相同结构。如果这完全是一个水平分布，那么这些点应该都非常接近这里的一条恒定线。我们确实看到了生物变异的证据，出现在不同基因组之间差异表达的基因中。然后我们看到了我们感兴趣去除的技术变异的特征，位于这些线上扬的斜率中。这意味着对于跨细胞基因的高平均计数，方差会增加。所以在更高的表达水平，我们也有更高的方差。因此，数据在统计学上被称为过离散。
建模与归一化
现在，这与测序深度问题有何关系？对于一个特定的基因，假设我们发现该基因的计数是两个因素的乘积。一个因素是这里的N，这是测序深度。另一个是P，这是概率或该基因在细胞表达中所占的比例。
假设在某种细胞类型中，来自该特定基因的转录本百分比是1%。如果我们有两个这种类型的细胞，对于一个细胞，测序深度为1000，那么我们的预期计数是10。如果测序深度为10000，那么我们的预期计数是100。从我们一开始已经查看过的测序深度变化这一事实，我们可以推断，与纯粹的泊松变量相比，我们将需要有一些过离散。
一个常见的用于模拟相对于泊松过离散的分布是所谓的负二项分布。在负二项分布中，你确实假设你的计数遵循泊松分布。但然后速率参数也是一个随机变量，它从伽马分布中抽取。这导致了额外的变异性，超过了纯泊松分布的变异性，并且在高均值处特别明显。
负二项分布看起来有点像这样。这里我们有一个小动画，取决于负二项分布的参数。它很适合模拟我们在数据中也看到的这些长尾。
现在，这是我之前展示给你的离散系数图。如果我们通过测序深度对计数进行归一化，在这种情况下，我们简单地将每个基因除以测序深度，然后乘以中位数测序深度，那么这里的许多弯曲形状实际上就消失了。
然后我们可以使用这样的图，例如，显示一个基因的离散系数与平均表达水平的关系，来选择我们特别感兴趣的基因，因为我们知道这些基因独立于它们显示的确切表达水平，表现出比我们通常预期更高的变异性。所以这些基因是不同细胞类别或细胞类型之间真正生物变异的良好候选者。
如果你对如何预处理单细胞转录组学数据或如何为单细胞转录组学数据计算良好分布以选择基因进行进一步分析这些主题更感兴趣，你也可以参考我们最近的预印本，该预印本可在bioRxiv上找到，我们在其中更深入地处理了这些情况。
方差稳定变换
好的，现在我们已经完成了第一个预处理步骤。我们已经通过测序深度对计数进行了归一化，并将它们重新缩放到相同的尺度，以解释数据中技术变异的一个来源。
单细胞转录组学数据经常做的另一个预处理步骤是应用非线性方差稳定变换，例如平方根变换或log(x+1)变换。这两者都具有将方差带到更相似范围的效果。这实际上是有用的。
例如，如果我们在该域中的原始数据矩阵上进行主成分分析，并简单地在原始数据上绘制前两个主成分，数据看起来会有点像这样。这里不同的颜色代表不同的细胞类别，数据重叠相当严重。
如果我们先进行平方根变换，然后再进行PCA，你已经可以看到现在这个图揭示了一些不同的细胞簇之间的明显差异。
最后，如果我们应用log(x+1)变换，这似乎更适合，那么即使在简单的PCA图中，这些不同的簇也能更清晰地分辨出来。
因此，应用这种非线性方差稳定变换或多或少是任何转录组学分析流程的标准步骤。
任务二：数据可视化与降维 📉
我们在处理单细胞转录组学数据时的下一个目标是创建一个低维嵌入，原则上我们可以使用PCA，但也可以使用其他更非线性的算法。我们为什么要这样做？因为数据的维度和缺乏固有结构，实际上很难查看转录组学数据。当然，我们可以简单地制作一个图，可视化每个细胞中每个基因的平均表达。但对于肉眼来说，除非我们以良好的方式对数据进行排序，否则很难看出太多。当然，当我们查看图像或时间轨迹时，情况就非常不同了，因为图像和时间序列具有固有的内部结构。
因此，创建这些低维嵌入对于单细胞转录组学数据的探索性数据分析尤为重要。
高维数据降维的挑战
现在，为什么这甚至是一个困难的任务？对于那些听过我们机器学习课程或了解t-SNE及相关算法的人来说，可以跳过接下来的几分钟或简要听一下。
很难对高维数据进行良好的低维嵌入，因为不可能保留高维距离。如果你在二维中查看，并从标准高斯分布中采样点，在二维中，你查看成对距离，那么这就是你会得到的直方图。如果你进入50维，那么直方图，即点对之间的平均距离已经大得多，平均为10。如果你进入更高维的状态，比如500维，那么平均距离变得更大得多。
在保留所有这些成对距离的同时将其嵌入二维根本不可行。
因此，一类流行的嵌入算法（称为邻域嵌入算法）背后的一个想法是保留邻居，而不是绝对距离。一个例子是流行的算法，称为t分布随机邻域嵌入或t-SNE。
其思想如下：对于一个给定的点i，你查看其封闭邻域内的最近邻居。对于这些点，你有一个损失函数，它基本上是高维空间中的成对相似性（即p_ij）与低维空间中的成对相似性之间的Kullback-Leibler散度。但与其他真正试图保留距离的算法相比，这里的技巧是，我们主要只关心局部邻域。
损失函数包含我刚才提到的这两个部分。我们使用这些指数核来测量高维相似性。首先，它们是有向的，但我们之后会对它们进行对称化。这个核宽度通过一个称为困惑度的参数选择，困惑度设置这个核宽度，使得局部邻域中有足够数量的点。
我们也测量低维相似性。我们使用某个核函数。在t-SNE中，这个核函数是柯西或t分布核。然后，我们也可以计算损失函数的梯度。关于这个梯度有趣的一点是，在你观察它一段时间后，它基本上由两项组成，一项是相似点被吸引，另一项是不相似点相互排斥。所以这是一个吸引-排斥的权衡。
我提到的这个参数困惑度控制核宽度。这里，我们有非常小的困惑度。这里，我们有更高的困惑度。许多实现中的默认值实际上是30左右，我们这样计算困惑度。它是某个熵值的幂，该熵值由这些成对相似性给出。
通过控制困惑度，你可以控制算法的分辨率，或者它权衡局部邻域与全局排斥的程度。
t-SNE优化与初始化
这里我们查看t-SNE的优化过程。如果我们只是针对该损失函数优化点的位置，那么你可以看到它确实需要一段时间，这里有很多次迭代，并且它在优化过程中经历不同的阶段。这里我使用的数据是MNIST数据集。所有这些不同的簇对应不同的数字。你可以看到随着优化的进行，该数据集中的结构是如何出现的。
该优化问题的问题在于它是非凸的。所以，如果你简单地使用许多标准实现，那么这些t-SNE实现可能不擅长保留全局几何结构。例如，如果我们查看这项同样非常引人注目的研究“跨新皮质区域的共享和不同的转录组学细胞类型”，那么这里有一个非常美丽的可视化。
但如果你查看GABA能神经元的位置，在这种情况下，GABA能中间神经元分布在两个大岛中。在它们之间，是谷氨酸能神经元。当然，在全局尺度上，谷氨酸能神经元与这些GABA能神经元或那些GABA能神经元（它们彼此更密切相关）的差异要大得多。但算法却将它们放在中间。这是一个不希望出现的特性。
让算法至少更好地尊重这种全局关系的一种方法是，用一个相当好的解决方案初始化算法，例如通过PCA。
如果你用PCA坐标初始化算法，那么它已经处于非线性优化之后发生的正确值范围内。所以你可以在这里看到，在这种情况下，优化是如何运行的。这里有来自MNIST数据集的三个不同的“岛”，包含不同的数字类别。因为我们通过PCA初始化，然后算法的后期阶段更详细地塑造了这些不同岛屿和子岛屿的外观。
任务三：细胞类型识别（聚类）🔍
好的，现在我们有了一个嵌入。如果你对t-SNE的更多细节感兴趣，我已经为你链接了一个视频。
所以下一步将是识别细胞类型。这基本上是聚类的一个实例。
我们在课程中已经多次做过聚类，例如，在使用高斯混合算法寻找电生理记录中的细胞时。原则上，当然也可以在这里使用这样的混合分布，并调整所使用的分布。高斯混合可能不好，因为计数是离散的，并且均值相对较低，使得高斯近似不是一个好的拟合。但是，例如，可以使用负二项混合分布，这也被提出过。
与我们在其他情况下讨论的许多情况一样，我们事先不知道簇的数量，许多人更喜欢使用非参数聚类技术，这些技术不像负二项分布那样做出明确的建模假设。一个算法需要能够处理数百万个细胞，所以它需要非常高效，并且需要能够在高维数据上运行，因为即使在降维之后，经过基因选择，你查看的可能只有几千个基因，我们确实有多达30、40、50个主成分的信息量。
该领域最流行的一类聚类算法是使用数据的K最近邻图的社区检测算法。
K最近邻图与社区检测
什么是K最近邻图？例如，对于这个人工数据集。它只是一个图，其中对于每个数据点，你简单地将其连接到最密切相关的10个邻居（取决于你的K值）。在这个图中，你已经可以看到相当多的结构。在上面这里，连接非常密集，在这个岛中，连接非常密集，这里也是如此。很少有连接超出这两个岛屿。
因此，可以想象，如果你对这个K最近邻图进行社区检测，你会将原始簇作为这些社区找到。
这类社区检测算法在图中用于检测簇时经常优化的成本函数是这样的。这称为模块度，模块度测量的是：一对节点之间存在的边数，并将其与在随机图中如果它们在同一社区中你预期会有的边数进行比较。所以这试图达到的是，平均而言，真正紧密连接或更紧密地聚集在一起的社区，比你从随机图中预期的要多。
Louvain 与 Leiden 算法
有两种流行的算法来解决此类问题。第一个提出的算法称为Louvain算法。我把伪代码放在这里。
Louvain算法背后的思想如下：你从节点的随机社区分配开始。然后在第一步，你以随机顺序访问每个节点。你检查，并将其移动到可能不同的社区，如果这改进了成本函数。你这样做的方式在这里的move_node函数中指定。这里，对于每个被访问的节点，你检查如果将其分配到与其他节点相同的社区或其中一个其他社区，成本函数的增量是多少。如果这个增量大于0，那么你就进行分配；如果不大于0，那么你就不分配。你迭代这个过程。如果你完成了所有节点，那么你执行一个聚合步骤，通过该步骤你将单个社区中的所有节点组合成一个更高级别的图。然后你一遍又一遍地重复相同的过程。当这个过程收敛时，你就完成了社区检测。
这里也是聚合步骤的伪代码，它基本上说，如果两个节点在同一个社区，那么将它们组成新图的一个元节点。
然而，这个简单算法有一些问题，其中之一是这个：可能因为你以随机顺序访问节点，最终会出现一种情况，即一个社区仅由这里的节点0维系，并且算法的工作方式可能更好。可以证明，在某些情况下，当你最后访问这个节点时，将其移动到另一个社区可能是有益的。然后你最终会为rat节点得到一个不连通、分裂的社区。
这可以通过所谓的Leiden算法来修复。Leiden算法只是有一个额外的步骤，称为细化分区，确保这种情况不会发生。
聚类结果示例
如果你将这种Leiden聚类应用到我们一直使用的标准示例数据集（我们的Harris数据集）上，那么你可以看到它在这里产生了一个相对合理的聚类。这也与我们在该数据集的低维可视化中看到的情况匹配得很好。
总结 📝
本节课中我们一起学习了单细胞转录组学的基础知识。我们从定义和原理入手，了解了如何通过实验技术获取单细胞基因表达数据。我们深入探讨了处理此类数据时的核心挑战，特别是其高维、稀疏且缺乏固有排序的特点。
我们系统性地学习了三个核心数据分析任务：
- 数据预处理：包括测序深度归一化以校正技术变异，以及应用方差稳定变换（如log(x+1)）来改善后续分析。
- 数据可视化与降维：重点介绍了t-SNE算法，它通过保留局部邻域结构来创建有意义的低维嵌入，用于探索性数据分析。
- 细胞类型识别（聚类）：介绍了基于K最近邻图和模块度优化的社区检测算法（如Leiden算法），这是一种高效且适用于高维数据的非参数聚类方法。
单细胞转录组学领域正在迅速发展，未来将在神经科学中带来更多应用，并成为现代细胞水平神经科学实验技术的基石。因此，思考我们在分析端所需的专门数据科学技术非常重要。通过本讲座以及你将在家完成的编程练习，希望我们为你在这个方向上提供了一个良好的起点。
感谢你的关注，我们下周再见。
012：神经形态分析
在本节课中，我们将学习如何分析单个神经元的神经解剖结构，即神经形态。我们将探讨用于表征神经形态的实验技术、从图像数据中重建神经元形态的挑战与方法，以及如何分析和表示这种独特的图结构数据。
概述
神经元的解剖结构自神经科学诞生之初就令科学家着迷。本节课将讨论用于表征神经解剖学的实验技术、分析单个神经元形态时需要注意的特殊性，以及从原始图像到抽象三维结构表示的处理流程。
实验技术
以下是几种用于观察神经解剖结构的主要实验方法。
传统染色技术
- 尼氏染色：使用甲苯胺蓝等染料，主要染色细胞体，无法显示构成神经元功能和特性的详细突起结构。
- 高尔基染色：该方法历史悠久，被圣地亚哥·拉蒙·卡哈尔广泛使用。组织切片先后浸入重铬酸钾和硝酸银，能稀疏地染色单个但相对完整的神经元，从而允许进行形态描绘。
- 生物胞素染色：向神经元内注射生物胞素，其本身不是染料，但后续可通过荧光团或过氧化物酶染色进行标记，从而显示单个神经元的形态。
从图像到重建
从染色图像或图像堆栈中重建单个神经元的三维结构是一个关键步骤。最终我们需要的是神经元三维结构的抽象表示，而非仅仅是染色神经元的图像。
目前在许多实验室中，手动追踪仍然是金标准。研究人员使用如NeuronJ等软件工具，在加载的图像堆栈中手动添加追踪路径。这个过程非常耗时，容易出错，且缺乏实验室间的标准化。
自动追踪的挑战
实现全自动追踪面临诸多挑战：
- 数据量巨大：为以足够分辨率捕获单个神经元，需要获取三维图像堆栈和不同成像视野的马赛克拼接。
- 缺乏训练数据：目前可用于训练的大规模数据集很少。
- 成像伪影与结构模糊：包括分辨率极限下的结构难以区分、突触棘与分支的混淆、不同细胞区室的交叉以及染色不连续（如串珠状或存在间隙）等问题。
自动追踪工作流
尽管困难，多年来已提出多种自动追踪工作流。它们通常包含两个部分：
- 预处理：校正成像伪影、抑制噪声、进行图像滤波以选择可能与神经元相关的部分。
- 分割与建树：将图像分割成与神经元解剖结构一致的合理树状图。算法可分为两类：
- 全局方法：通过二值化、骨架化等步骤生成神经元的图表示。
- 局部方法：预测树上的下一个可能点，并拟合神经元形态可能如何发展的内部模型。
神经元本质上是一种图结构数据，每个节点具有三维坐标，最适合用图来表示。
一个自动追踪算法示例：Snake Tracer
Snake Tracer是一个相对简单但效果不错的算法示例。
- 生成二值掩膜：
- 对图像进行平滑背景减去，以锐化边界并消除背景中的低频变化。
- 使用一系列“谷检测器”（实质上是高斯滤波器）与图像进行卷积。
- 对响应图像进行阈值处理，得到二值掩膜。
- 通过平滑或噪声抑制去除背景中不关心的部分。
- 生成骨架：
- 对二值掩膜计算距离变换（到非掩膜部分的最短距离）。
- 沿着距离变换的脊线（最大值）放置单一路径，作为骨架。
- 存储为SWC格式：骨架以SWC文件格式存储，这是神经元形态数据的标准格式。它使用一个简单的表格记录每个节点的：
- ID编号
- 类型（如胞体、树突、轴突）
- 三维坐标 (X, Y, Z)
- 半径
- 父节点ID（用于定义连接关系）
- 后处理与校正：检查神经元的连续性，计算结构宽度，并进行必要的校正（如连接断开的片段或删除噪声）。最后通常仍需结合手动校正步骤。
公共数据库与质量控
如果想研究形态数据，可以访问公共数据库，如NeuroMorpho.org，它包含来自不同物种、脑区和实验室的大量数据。使用这些数据时，必须进行严格的质量控制，筛选出适合自己研究目的且质量足够好的数据集。
另一个更聚焦的数据库是艾伦脑科学研究所的细胞类型数据库，它提供质量一致且较高的数据，通常还结合了电生理学或转录组学数据。
电子显微镜技术
作为替代技术，电子显微镜（EM）可用于神经形态分析，特别是近20年来发展起来的高通量EM技术，如连续块面扫描电子显微镜。
该技术将包埋的组织块进行逐层切削和成像，最终获得能分辨亚细胞结构的三维组织图像。其最大挑战在于自动重建，因为需要对密集染色的组织中的所有结构进行分割。据估计，完全重建一个100微米立方体内的所有细胞可能需要约10万小时的人工劳动，因此需要高效的算法。
深度学习分割：Flood-Filling网络
一种基于深度学习的解决方案是Flood-Filling网络。
- 工作原理：在每次迭代中，网络接收图像和一个初始掩膜（从某个种子点开始），然后处理图像和掩膜以扩展掩膜。这个过程迭代进行，直到整个对象被掩膜覆盖。
- 训练：网络在由专家标注的小型数据子集上进行训练。
- 启发式策略处理错误：
- 避免错误合并：使用集成方法，从多个起始点、不同分辨率进行分割，并保留最精细的对象。
- 避免错误分裂：检查所有相邻的分割区域对，在边界附近初始化分割，如果网络将其分割为一个对象，则合并它们。
通过这些策略，该算法能从EM数据中获得错误率很低的自动重建结果。
形态数据的分析与表示
从手动或自动重建中，我们最终获得神经元的图结构表示。分析这种图结构数据有多种方法：
形态计量学统计与分布
计算简单的度量来量化神经元的各个方面，例如：
- 神经元沿皮质轴的高度或宽度。
- 主要分支从胞体分出的角度分布。
- 树中的分支总数、节点最高阶数等。
密度图
将图结构转换为图像表示。具体方法是：在一个离散网格中，统计每个小盒子里包含的神经元长度，从而得到神经元的像素化表示。这通常是一种简单且信息量丰富的表示方法。
基于图结构的深度学习方法
为了更好地利用神经元的图结构本质，我们开发了一种变分自编码器（VAE）架构。
- 编码器：通过在图结构上进行随机游走，并使用堆叠的长短期记忆（LSTM）单元进行处理，将神经元映射到一个潜在空间表示。
- 解码器：从潜在表示解码，生成新的随机游走，用以重建神经元形态。
- 结合细胞类型标签：可以在架构中添加分类器，利用已知的细胞类型标签来微调潜在表示，提高其判别能力。
这种方法的优点是：
- 能学习到比简单统计或密度图更好的表示。
- 可以生成与给定神经元形态相似但又不完全相同的新的人工神经元图，实现了生成模型的功能。
总结
本节课我们一起学习了神经形态分析的完整流程。
我们首先回顾了用于获取神经元形态的实验技术，包括尼氏染色、高尔基染色和生物胞素染色。接着，我们探讨了从这些图像数据中重建神经元三维结构所面临的挑战，无论是耗时的手动追踪还是面临数据量、伪影和缺乏标注等困难的自动追踪。
我们介绍了一个自动追踪算法（Snake Tracer）的示例和标准的SWC数据格式。此外，我们还了解了电子显微镜技术和先进的深度学习分割方法（如Flood-Filling网络）如何助力大规模神经回路的重建。
最后，我们重点讨论了如何分析和表示这种独特的图结构数据。从传统的形态计量学统计和密度图，到我们开发的基于图结构的变分自编码器，这些方法使我们能够量化、比较甚至生成神经元的形态。该领域仍然存在许多挑战和机遇，例如处理更大规模的组织块数据和学习更强大的形态表示，这些都是神经数据科学未来令人兴奋的方向。
013：神经科学中的数据管理 🗂️
在本节课中，我们将学习神经科学研究中的数据管理。我们已经探讨了多种机器学习算法及其在数据上的应用。然而，在实验室环境及更广泛的场景中，有效使用数据的一个重要组成部分是数据管理的标准。我们将了解如何保持数据的整洁有序、易于访问，以及如何跟踪结果以确保其可重复性。
现代实验室的数据挑战
正如我们在课程中所见，在现代实验室环境中，我们记录着海量信息。例如，在一次实验中，我们可能同时监测小鼠在跑步机上的运动、呈现视觉刺激、进行钙成像、记录动物的眼动和胡须活动，甚至用电极同步记录单个神经元的膜电位或局部场电位。此外，我们还有大量元数据，例如动物编号、所属转基因品系、使用的刺激参数等等。
随着我们尝试不同的预处理参数设置或后续的推理方法，这个过程会变得极其复杂。因此，当你处理的数据超出几个简单示例时，一个关键问题就是：如何将所有信息保持得井井有条？
传统工作流程及其问题
传统的工作流程通常如下：
- 实验结束后，将原始文件存放在实验室存储服务器上。
- 在实验记录本上写下大量元数据（如果一切顺利的话），或者仅凭记忆，有时甚至会忘记记录。
- 使用实验室通常采用的参数设置对数据进行后处理。这可能是手动完成，或通过运行某个函数。这个函数通常是已离开实验室的博士生编写的，并且后续可能经过了多次修改。
- 进行分析时，在后处理数据上运行脚本，并将中间结果保存在本地或存储服务器上。
- 每次需要生成图表时重新运行分析，或者如果数据量太大，则保存结果并在发现错误时重新运行。
- 最终，手动为论文整理所有能找到的数据。
当然，这是对常规实验室工作略带夸张的描述。我知道实验和计算实验室的同事们都非常认真地跟踪数据。但处理真实数据本身就是一个复杂的过程，即使非常努力，如果全靠手动操作，最终仍会丢失大量信息。
以下是传统工作流程中常见的问题：
- 代码版本混乱：不确定计算感受野时使用了哪个版本的代码。如果需要用当前版本的代码为大量细胞重新计算，这将耗费很长时间。
- 数据丢失或混乱：不确定是否删除了保存神经元1到42感受野的文件。如果找不到，就需要从头重新计算，这在时间紧迫时（例如会议截稿前或准备报告时）是难以承受的。
- 数据顺序不匹配：不确定附加的细胞标签（例如标记哪些细胞被转基因品系标记）是否顺序正确，以及是否与感受野数据的顺序匹配。
- 元数据缺失：需要查找大脑中的确切坐标、皮层位置、视网膜位置或记录时的动物体温等信息时，可能发现实验记录本中没有相关信息，或者当时只检查了大致范围而未记录具体测量值。
解决方案：数据管理系统
那么，我们能做些什么呢？一种方法是使用数据管理系统。其中一种对神经科学实验室（以及其他领域，如深度学习模型训练、生物物理模型模拟或任何其他类型的实验数据）非常实用的系统叫做 DataJoint。
在DataJoint中，数据科学家通过编写数据管道来指定数据处理流程。有一个客户端（支持MATLAB或Python）管理与后端架构的交互，后端可以是用于大文件的二进制对象存储，也可以是用于较小数据的MySQL服务器。
这种数据管理系统的重要之处在于，一旦数据和元数据进入系统，它就能保持它们之间的依赖关系。它能自然地与MATLAB和Python接口，并拥有简单高效的查询语言，可以轻松直接地进行分布式计算。
DataJoint的核心：关系数据模型
DataJoint建立在关系数据模型之上。如果你熟悉数据库，这对你来说并不陌生。
关系数据模型需要几个基本概念：
- 元组：一组属性名-值对。例如：
(mouse_id: 1001, measure_date: ‘2023-10-10’, weight: 25)。 - 关系：一组具有相同属性名的元组的集合。关系中的每个元组必须是唯一的。你可以将关系看作一个表。
- 模式：一组逻辑上相关的关系的集合。
- 主键：用于唯一标识一个元组的属性子集。例如，在一个包含
mouse_id、measure_date和weight的关系中，mouse_id和measure_date可以共同作为主键。 - 外键：从另一个关系中继承的属性子集。这确保了数据的引用完整性。
引用完整性的重要性
让我们通过一个具体例子来理解外键和引用完整性。
假设我们有一个名为 mouse_info 的关系（表），主键是 mouse_id，它还包含老鼠的性别信息。
表：mouse_info
| mouse_id | sex |
|----------|--------|
| 1001 | Male |
| 1002 | Female |
然后我们有第二个关系 mouse_weight，它继承了 mouse_info 中的 mouse_id 作为外键，并记录在特定日期测量的体重。
表：mouse_weight
| mouse_id | measure_date | weight |
|----------|--------------|--------|
| 1001 | 2023-10-10 | 25 |
| 1002 | 2023-10-11 | 22 |
外键的作用：
- 强制元数据录入：如果你想在
mouse_weight表中为某只老鼠录入体重数据，系统会要求这只老鼠必须首先存在于mouse_info表中。这迫使你必须先录入老鼠的性别等元数据，避免了遗忘。 - 保证引用完整性：如果你从
mouse_info基表中删除了一只老鼠（例如，mouse_id: 1001），系统会自动确保mouse_weight表中对应的条目也被删除。这防止了存在引用“孤儿”数据的情况。
这有什么好处？ 映射到感受野计算的例子：如果你在数据预处理（如尖峰排序）中发现了一个错误并进行了修改，那么所有基于错误预处理结果计算的感受野数据都会被自动标记为无效或删除。这样，你就不会遇到不确定某个感受野文件是否是用有错误的代码计算出来的情况。
DataJoint的便捷操作
这样的系统通过定义一些映射到数据库操作的关键运算符，使得数据过滤和操作变得非常容易。
以下是几种核心操作：
- 匹配：如果两个元组在相同命名的属性上具有相等的值，则它们匹配。例如，
mouse_info表和mouse_weight表通过相同的mouse_id进行匹配。 - 限制：选择符合特定条件的元组，即过滤表格。例如，选择所有在
2010-10-10测量了体重的老鼠。 - 连接：连接两个关系，生成一个包含所有匹配元组及双方原始信息的新关系。例如，连接
mouse_info和mouse_weight表，得到同时包含性别和体重信息的新表。
自动化计算与数据存储
DataJoint 使得实现自动化计算变得简单。你可以将一些关系指定为“手动”关系，数据需要手动录入或由实验系统自动录入。然后，你可以定义“计算”关系，其中的数据是其他关系中信息的简单函数。
例如，定义一个 mouse_bmi 表，其数据基于 mouse_height 和 mouse_weight 表中的信息计算得出。你只需在 mouse_bmi 表中指定计算规则，然后调用一个简单的 populate() 函数。系统会自动检查需要填充的条目，并为这两个基表中所有唯一的组合执行计算。
数据存储方式：
- 每个关系映射到 MySQL 数据库中的一个表。
- 数据库中的每个表都与一个 MATLAB 或 Python 对象关联。
- 用户通过调用这些封装了对应表的 Python 或 MATLAB 对象来操作数据。
- 对象本身不包含数据，必要时从数据库中获取。你可以通过类似
fetch()的命令将数据调入工作空间进行处理。
与直接编写复杂的 SQL 查询相比，DataJoint 提供了更直观、更安全的接口。例如，在 DataJoint 中过滤所有 BMI 大于 28 的值非常简单，而对应的原生 SQL 语句则更复杂且容易出错。
此外，并行处理也很容易，因为填充某个元组中的一个条目可以在单个核心上完成，因此可以简单地使用计算集群来并行计算某个表的所有条目。系统还支持为关系或元组定义绘图函数，以便用标准化图表轻松浏览数据。
实际应用案例与总结
这是一个来自我们过去项目的简化示例。一个数据集包含多个象限的刺激呈现、每个象限的附加信息、每个双光子成像实验定义的感兴趣区域、针对每个区域和刺激组合提取的信号轨迹，以及对其中使用噪声刺激的轨迹子集基于钙信号计算的感受野。
通过将整个工作流程实现为 DataJoint 关系，我们获得了所有附加好处：引用完整性、轻松过滤、实验室成员的标准访问接口等。系统可以配置为一旦所有手动定义的组件就位，就自动执行后处理。
如果你想开始使用 DataJoint，有一个优秀的教程网站：tutorials.datajoint.org。这是一个开源系统，你可以免费下载和使用。
本节课是《神经数据科学》课程的最后一讲。我们探讨了神经科学研究中数据管理的重要性、传统工作流程的挑战，并介绍了一个强大的解决方案——DataJoint 关系数据管理系统。通过保持数据依赖性和引用完整性，该系统能帮助研究人员更可靠、更高效地组织、处理和复现复杂的实验数据。希望本课程让你对在神经科学乃至更广泛的科学领域中使用机器学习和数据科学的前景感到兴奋，并激励你在此方向继续探索。
感谢你的关注，请保持健康。
浙公网安备 33010602011771号