非官方中文翻译
本页译自 MNE-Python 官方文档：
查看英文原文 ↗。
文档版权归 MNE Developers 所有（BSD-3-Clause）。
说明
转到末尾下载完整示例代码。
预处理功能性近红外光谱（fNIRS）数据#
本教程介绍如何将功能性近红外光谱（fNIRS）数据从原始测量值转换为相对氧合血红蛋白（HbO）和脱氧血红蛋白（HbR）浓度，并查看平均波形以及响应在地形图上的表示。
这里我们将使用 fNIRS 运动数据。
# Authors: The MNE-Python contributors.
# License: BSD-3-Clause
# Copyright the MNE-Python contributors.
import matplotlib.pyplot as plt
import numpy as np
import mne
fnirs_data_folder = mne.datasets.fnirs_motor.data_path()
fnirs_cw_amplitude_dir = fnirs_data_folder / "Participant-1"
raw_intensity = mne.io.read_raw_nirx(fnirs_cw_amplitude_dir, verbose=True)
raw_intensity.load_data()
Loading /home/circleci/mne_data/MNE-fNIRS-motor-data/Participant-1
Reading 0 ... 23238 = 0.000 ... 2974.464 secs...
提供更有意义的注释信息#
首先，我们为以注释形式存储的触发码赋予更有意义的名称。其次，我们加入每个刺激的持续时间信息，本实验中所有条件的持续时间均为 5 秒。第三，我们移除触发码 15，它标记实验的开始与结束，与我们的分析无关。
raw_intensity.annotations.set_durations(5)
raw_intensity.annotations.rename(
{"1.0": "Control", "2.0": "Tapping/Left", "3.0": "Tapping/Right"}
)
unwanted = np.nonzero(raw_intensity.annotations.description == "15.0")
raw_intensity.annotations.delete(unwanted)
查看传感器在脑表面的位置#
这里我们验证源-探测器对和通道的位置是否符合预期。源-探测器对显示为光极之间的连线，通道（源-探测器对的中点）可选择显示为橙色圆点。源可选择显示为红色圆点，探测器则显示为黑色。
subjects_dir = mne.datasets.sample.data_path() / "subjects"
brain = mne.viz.Brain(
"fsaverage", subjects_dir=subjects_dir, background="w", cortex="0.5"
)
brain.add_sensors(
raw_intensity.info,
trans="fsaverage",
fnirs=["channels", "pairs", "sources", "detectors"],
)
brain.show_view(azimuth=20, elevation=60, distance=400)
Channel types:: fnirs_cw_amplitude: 56
选择适合检测神经响应的通道#
首先，我们移除距离过近（短通道）而无法检测神经响应的通道（光极间距小于 1 cm）。这些短通道可在上图中看到。为此，我们挑选所有不被视为短通道的通道。
picks = mne.pick_types(raw_intensity.info, meg=False, fnirs=True)
dists = mne.preprocessing.nirs.source_detector_distances(
raw_intensity.info, picks=picks
)
raw_intensity.pick(picks[dists > 0.01])
raw_intensity.plot(
n_channels=len(raw_intensity.ch_names), duration=500, show_scrollbars=False
)
从原始光强转换为光密度#
随后将原始光强值转换为光密度。
raw_od = mne.preprocessing.nirs.optical_density(raw_intensity)
raw_od.plot(n_channels=len(raw_od.ch_names), duration=500, show_scrollbars=False)
评估数据质量#
在此阶段，我们可以使用头皮耦合指数（SCI）量化头皮与光极之间耦合的质量。该方法在两次光检测信号中寻找心脏信号频段内显著的同步成分。
在本示例中数据干净，所有通道的耦合都良好，因此我们不会基于头皮耦合指数将任何通道标记为坏通道。
sci = mne.preprocessing.nirs.scalp_coupling_index(raw_od)
fig, ax = plt.subplots(layout="constrained")
ax.hist(sci)
ax.set(xlabel="Scalp Coupling Index", ylabel="Count", xlim=[0, 1])
在本示例中，我们会将所有 SCI 小于 0.5 的通道标记为坏通道（该数据集相当干净，因此没有通道被标记为坏通道）。
raw_od.info["bads"] = [raw_od.ch_names[idx] for idx in np.where(sci < 0.5)[0]]
此时适合检查数据（有关如何使用交互式数据可视化工具的说明，请参见 Raw 对象的内置绘图方法），以确保头皮耦合较差的通道已被移除。如果数据包含大量伪迹，您可能决定应用伪迹去除技术，如“可视化 NIRS 伪迹校正方法”中所述。
从光密度转换为血红蛋白浓度#
接下来，我们使用改良 Beer-Lambert 定律将光密度数据转换为血红蛋白浓度。
raw_haemo = mne.preprocessing.nirs.beer_lambert_law(raw_od, ppf=0.1)
raw_haemo.plot(n_channels=len(raw_haemo.ch_names), duration=500, show_scrollbars=False)
从信号中去除心率成分#
血流动力学响应的频率成分主要低于 0.5 Hz。数据中约 1 Hz 处可见的活动增强源于被试的心跳，属于不需要的成分。因此我们使用低通滤波器将其去除。同时还加入高通滤波器以去除数据中的缓慢漂移。
raw_haemo_unfiltered = raw_haemo.copy()
raw_haemo.filter(0.05, 0.7, h_trans_bandwidth=0.2, l_trans_bandwidth=0.02)
for when, _raw in dict(Before=raw_haemo_unfiltered, After=raw_haemo).items():
fig = _raw.compute_psd().plot(
average=True, amplitude=False, picks="data", exclude="bads"
)
fig.suptitle(f"{when} filtering", weight="bold", size="x-large")
Filtering raw data in 1 contiguous segment
Setting up band-pass filter from 0.05 - 0.7 Hz
FIR filter parameters
---------------------
Designing a one-pass, zero-phase, non-causal bandpass filter:
- Windowed time-domain design (firwin) method
- Hamming window with 0.0194 passband ripple and 53 dB stopband attenuation
- Lower passband edge: 0.05
- Lower transition bandwidth: 0.02 Hz (-6 dB cutoff frequency: 0.04 Hz)
- Upper passband edge: 0.70 Hz
- Upper transition bandwidth: 0.20 Hz (-6 dB cutoff frequency: 0.80 Hz)
- Filter length: 1291 samples (165.248 s)
Effective window size : 262.144 (s)
Plotting power spectral density (dB=True).
Effective window size : 262.144 (s)
Plotting power spectral density (dB=True).
提取 Epochs#
现在信号已转换为相对血红蛋白浓度，且不需要的心率成分已被去除，我们可以提取与各个实验条件相关的 epochs。
首先，我们提取感兴趣的事件并将其可视化，以确认它们正确无误。
events, event_dict = mne.events_from_annotations(raw_haemo)
fig = mne.viz.plot_events(events, event_id=event_dict, sfreq=raw_haemo.info["sfreq"])
Used Annotations descriptions: ['Control', 'Tapping/Left', 'Tapping/Right']
接下来，我们定义 epochs 的时间范围、拒绝标准、基线校正，并提取 epochs。我们可视化被丢弃 epochs 的日志。
reject_criteria = dict(hbo=80e-6)
tmin, tmax = -5, 15
epochs = mne.Epochs(
raw_haemo,
events,
event_id=event_dict,
tmin=tmin,
tmax=tmax,
reject=reject_criteria,
reject_by_annotation=True,
proj=True,
baseline=(None, 0),
preload=True,
detrend=None,
verbose=True,
)
epochs.plot_drop_log()
Not setting metadata
90 matching events found
Setting baseline interval to [-4.992, 0.0] s
Applying baseline correction (mode: mean)
0 projection items activated
Using data from preloaded Raw for 90 events and 157 original time points ...
Rejecting epoch based on HBO : ['S4_D4 hbo']
Rejecting epoch based on HBO : ['S4_D4 hbo', 'S8_D8 hbo']
Rejecting epoch based on HBO : ['S4_D4 hbo']
Rejecting epoch based on HBO : ['S4_D4 hbo', 'S8_D8 hbo']
Rejecting epoch based on HBO : ['S1_D1 hbo', 'S3_D3 hbo', 'S4_D4 hbo', 'S7_D6 hbo', 'S7_D7 hbo', 'S8_D8 hbo']
Rejecting epoch based on HBO : ['S4_D4 hbo', 'S6_D8 hbo', 'S8_D8 hbo']
6 bad epochs dropped
查看跨试次响应的一致性#
现在我们可以查看敲击条件的血流动力学响应。我们分别可视化氧合血红蛋白和脱氧血红蛋白的响应，可以观察到 HbO 在约 6 秒处出现符合预期的峰值，且在各试次间保持一致；同时 HbR 出现一致的下降，并相对于 HbO 峰值略有延迟。
epochs["Tapping"].plot_image(
combine="mean",
vmin=-30,
vmax=30,
ts_args=dict(ylim=dict(hbo=[-15, 15], hbr=[-15, 15])),
)
Not setting metadata
55 matching events found
No baseline correction applied
0 projection items activated
Not setting metadata
55 matching events found
No baseline correction applied
0 projection items activated
combining channels using "mean"
combining channels using "mean"
我们还可以查看对照条件的分段数据，并观察到它并未表现出预期的形态。
epochs["Control"].plot_image(
combine="mean",
vmin=-30,
vmax=30,
ts_args=dict(ylim=dict(hbo=[-15, 15], hbr=[-15, 15])),
)
Not setting metadata
29 matching events found
No baseline correction applied
0 projection items activated
Not setting metadata
29 matching events found
No baseline correction applied
0 projection items activated
combining channels using "mean"
combining channels using "mean"
查看跨通道响应的一致性#
同样，我们可以查看所选光极对之间响应的一致性。该数据中的所有通道都位于运动皮层上方，且所有通道都显示出相似的数据模式。
fig, axes = plt.subplots(nrows=2, ncols=2, figsize=(15, 6), layout="constrained")
clims = dict(hbo=[-20, 20], hbr=[-20, 20])
epochs["Control"].average().plot_image(axes=axes[:, 0], clim=clims)
epochs["Tapping"].average().plot_image(axes=axes[:, 1], clim=clims)
for column, condition in enumerate(["Control", "Tapping"]):
for ax in axes[:, column]:
ax.set_title(f"{condition}: {ax.get_title()}")
绘制标准 fNIRS 响应图像#
接下来，我们生成 fNIRS 数据最常见的可视化形式：在同一幅图中绘制 HbO 和 HbR，以展示两种信号之间的关系。
evoked_dict = {
"Tapping/HbO": epochs["Tapping"].average(picks="hbo"),
"Tapping/HbR": epochs["Tapping"].average(picks="hbr"),
"Control/HbO": epochs["Control"].average(picks="hbo"),
"Control/HbR": epochs["Control"].average(picks="hbr"),
}
# Rename channels until the encoding of frequency in ch_name is fixed
for condition in evoked_dict:
evoked_dict[condition].rename_channels(lambda x: x[:-4])
color_dict = dict(HbO="#AA3377", HbR="b")
styles_dict = dict(Control=dict(linestyle="dashed"))
mne.viz.plot_compare_evokeds(
evoked_dict, combine="mean", ci=0.95, colors=color_dict, styles=styles_dict
)
combining channels using "mean"
combining channels using "mean"
combining channels using "mean"
combining channels using "mean"
查看活动的地形图表示#
接下来，我们查看响应过程中地形活动如何变化。
times = np.arange(-3.5, 13.2, 3.0)
topomap_args = dict(extrapolate="local")
epochs["Tapping"].average(picks="hbo").plot_joint(
times=times, topomap_args=topomap_args
)
No projector specified for this dataset. Please consider the method self.add_proj.
比较左手与右手的敲击任务#
最后，我们为左、右两种条件生成地形图，以查看活动的位置。首先可视化 HbO 活动。
我们再来查看两种条件下的 HbR 活动。
epochs["Tapping/Left"].average(picks="hbr").plot_topomap(times=times, **topomap_args)
epochs["Tapping/Right"].average(picks="hbr").plot_topomap(times=times, **topomap_args)
我们还可以针对两种条件在单个时间点绘制比较图。
fig, axes = plt.subplots(
nrows=2,
ncols=4,
figsize=(9, 5),
gridspec_kw=dict(width_ratios=[1, 1, 1, 0.1]),
layout="constrained",
)
vlim = (-8, 8)
ts = 9.0
evoked_left = epochs["Tapping/Left"].average()
evoked_right = epochs["Tapping/Right"].average()
evoked_left.plot_topomap(
ch_type="hbo", times=ts, axes=axes[0, 0], vlim=vlim, colorbar=False, **topomap_args
)
evoked_left.plot_topomap(
ch_type="hbr", times=ts, axes=axes[1, 0], vlim=vlim, colorbar=False, **topomap_args
)
evoked_right.plot_topomap(
ch_type="hbo", times=ts, axes=axes[0, 1], vlim=vlim, colorbar=False, **topomap_args
)
evoked_right.plot_topomap(
ch_type="hbr", times=ts, axes=axes[1, 1], vlim=vlim, colorbar=False, **topomap_args
)
evoked_diff = mne.combine_evoked([evoked_left, evoked_right], weights=[1, -1])
evoked_diff.plot_topomap(
ch_type="hbo", times=ts, axes=axes[0, 2:], vlim=vlim, colorbar=True, **topomap_args
)
evoked_diff.plot_topomap(
ch_type="hbr", times=ts, axes=axes[1, 2:], vlim=vlim, colorbar=True, **topomap_args
)
for column, condition in enumerate(["Tapping Left", "Tapping Right", "Left-Right"]):
for row, chroma in enumerate(["HbO", "HbR"]):
axes[row, column].set_title(f"{chroma}: {condition}")
最后，我们还可以查看单个波形，以了解是什么驱动了上方的地形图。
fig, axes = plt.subplots(nrows=1, ncols=1, figsize=(6, 4), layout="constrained")
mne.viz.plot_evoked_topo(
epochs["Left"].average(picks="hbo"), color="b", axes=axes, legend=False
)
mne.viz.plot_evoked_topo(
epochs["Right"].average(picks="hbo"), color="r", axes=axes, legend=False
)
# Tidy the legend:
leg_lines = [line for line in axes.lines if line.get_c() == "b"][:1]
leg_lines.append([line for line in axes.lines if line.get_c() == "r"][0])
fig.legend(leg_lines, ["Left", "Right"], loc="lower right")
脚本总运行时间：（0 分钟 11.123 秒）