教程 | 近紅外數(shù)據(jù)的預處理和平均(上)

前言
近紅外光譜(NIRS)是一種測量流經(jīng)傳感器所在組織的血液中氧合水平的方法。它基于這樣一個事實,即含氧血紅蛋白和脫氧血紅蛋白具有不同的吸收光譜,因此你會看到它有不同的顏色。大多數(shù)近紅外系統(tǒng)在每個光源光電二極管發(fā)射2個波長的光,通常約700和850nm,有些系統(tǒng)則使用3個波長。光電探測器被放置在離源光極幾厘米遠的地方。光源和探測器的光電二極管共同組成一個通道,但系統(tǒng)對于一個通道總是有兩個空間位置(大多數(shù)人認為記錄的活動主要來自光源和探測器的中間位置)和兩個波長。光散射穿過組織,一小部分光將通過光電探測器在組織的表面被探測到。在頭皮上檢測到的光強隨著位于源和探測器光電二極管之間的組織中含氧和脫氧血紅蛋白濃度的變化而有系統(tǒng)地波動。由于(或多或少的隨機)散射,大多數(shù)到達探測器的光只穿過了組織表面,只有一小部分穿過了表面下更深一點的區(qū)域。通常情況下,在頭皮上檢測到的光穿過了一個組織體積,從側(cè)面看是一個香蕉形狀。如果光源和探測器被放置得更遠,那么被探測到的光可能在組織中傳播得更深。但實際上,這是有限的,因為到達探測器的那部分光會隨著光源-探測器距離的增加而急劇下降。當測量大腦中血紅蛋白的變化時,顱骨作為一個重要的屏障,吸收了很多光。因此,對大腦深處進行掃描只有在頭骨較薄的嬰兒中才可行。
氧和脫氧血紅蛋白的濃度不能以絕對條件評估,這意味著所有的測量都將反映的是相對變化。因此,像EEG和MEG一樣,需要一個基線期。在血流動力學響應函數(shù)(HRF)上,記錄和后期處理的fNIRS信號與fMRI類似,但具有更高的時間分辨率。因此,至少在理論上,可以測量曲線的初始傾角。與fMRI相比,空間分辨率要差很多,通常我們主要從大腦表面進行記錄。因此,繪圖通常類似于EEG。EEG的時間分辨率更高,但fNIRS的一個優(yōu)勢是,測量的信號只來自光源和探測器之間的區(qū)域。fNIRS是一種非常有前景的方法,可以測量無法通過fMRI掃描儀測試的人群的大腦活動,并且允許被試做出無法在掃描儀中完成的更多動作。但是血壓會隨著你所做的任務而變化,從而也會影響測量到的信號。
單通道近紅外數(shù)據(jù)的預處理和平均
本文將演示如何分析單個通道的功能性近紅外光譜(fNIRS)數(shù)據(jù)。主要使用FieldTrip對Artinis近紅外光譜數(shù)據(jù)進行基本分析。由于本文分析的是單通道數(shù)據(jù),因而沒有展示如何處理壞通道,關于如何刪除壞通道和分析多通道的信息詳見近紅外數(shù)據(jù)的預處理和平均(下),關注茗創(chuàng)科技,不迷路哦~
本教程所使用的數(shù)據(jù)集
本文將使用到Artinis Medical Systems的Oxymon MK III系統(tǒng)采集運動皮層信息的近紅外數(shù)據(jù)集。光電極放置在運動皮層上,隨后要求被試進行手指敲擊。除了來自大腦的記錄,還記錄了一個帶有標注的事件(event)通道,在事件通道中,參與者開始敲手指時的記錄為“A”,停止敲手指的記錄為“B”。這個運動任務一共重復了12次,所有的數(shù)據(jù),包括事件通道,都保存在一個.oxy3文件中。本文中使用的數(shù)據(jù)可以從https://download.fieldtriptoolbox.org/tutorial/nirs_multichannel/上獲得。
分析步驟
根據(jù)數(shù)據(jù)和實驗設計的不同,分析可以有許多不同的方式和順序。這里將介紹分析步驟的標準順序,隨后小伙伴們可以逐步測試這些步驟。以下步驟為分析fNIRS數(shù)據(jù)提供了一個很好的標準方法,如圖1所示:1、讀取連續(xù)數(shù)據(jù);2、(優(yōu)先地)裁剪記錄開始和結(jié)尾的數(shù)據(jù),特別是當這些需要被剪切的數(shù)據(jù)是噪聲或不感興趣的片段時;3、刪除通道中的壞段(例如,運動偽影);4、將光密度轉(zhuǎn)化為含氧血紅蛋白(oxyHb)和脫氧血紅蛋白(deoxyHb)濃度的變化;5、將功能響應與系統(tǒng)反應分開,可以使用以下方法其中之一,或者組合使用;
①濾波
②減去參考通道
③每個通道都有反相關的oxyHb/deoxyHb6、在實驗任務中,定義與試次相對應的段;7、將試次數(shù)據(jù)平均并可視化。

本文分析的數(shù)據(jù)集名為“motor_cortex.oxy3”。首先,需要將FieldTrip指向這個文件。注意,F(xiàn)ieldTrip沒有圖形用戶界面(GUI),而是通過腳本來實現(xiàn)最終目標。需要將FieldTrip文件夾添加到MATLAB路徑中,然后執(zhí)行ft_defaults函數(shù),輸入:
ft_defaults
通常,F(xiàn)ieldTrip會將函數(shù)的所有輸入?yún)?shù)存儲在一個名為“configuration”或“cfg”的本地變量中。雖然MATLAB變量可以取任何名稱,但為了方便起見,這里將繼續(xù)使用名稱“cfg”。分析流程中的一個(關鍵)步驟是讀入數(shù)據(jù)并對數(shù)據(jù)進行預處理。為此,這里將使用FieldTrip函數(shù)ft_preprocessing。因為FieldTrip是一個開源工具箱,可通過輸入以下代碼來查看代碼細節(jié):
edit ft_preprocessing
讀取和裁剪數(shù)據(jù)
在變量cfg.dataset中輸入數(shù)據(jù)的文件名。
cfg = [];
cfg.dataset = 'motor_cortex.oxy3';
通常需要指定完整路徑,包括驅(qū)動器號。這可以確保代碼將獨立于當前的MATLAB位置執(zhí)行。將原始數(shù)據(jù)、處理過的數(shù)據(jù)和腳本保存在三個不同的位置是一個很棒的做法(推薦)。在本例中,將使用當前工作目錄中的數(shù)據(jù)。
現(xiàn)在,嘗試執(zhí)行ft_preprocessing:
[data] = ft_preprocessing(cfg);
在處理.oxy3文件時,需要加載包含光極布局的光極模板。如果MATLAB無法在你的路徑上找到模板文件,它將會彈出一個圖形用戶界面對話框,要求你定位XML文件。默認情況下,該文件位于C:\Program Files(x86)\Artinis Medical Systems BV\Oxysoft 3.0.103。需要選擇的文件是optodetemplates.xml。鑒于FieldTrip經(jīng)常會需要搜索此函數(shù),最好將此函數(shù)復制到存放MATLAB分析腳本或fNIRS數(shù)據(jù)的文件夾中。關于如何選擇最佳模板的信息,請參閱Artinis網(wǎng)頁上的這篇博客文章(網(wǎng)址:https://www.artinis.com/blogpost-all/2017/6/27/how-do-i-choose-the-correct-fibers-and-template-for-my-oxymon)。
當你在屏幕上看到類似這樣的內(nèi)容時,F(xiàn)ieldTrip就完成了。
>> the call to "ft_preprocessing" took 9 seconds
仔細觀察workspace,就會注意到添加了一個名為“data”的新變量。這個名稱與ft_preprocessing前面輸入的名稱相同。如果寫入[something],那么變量'something'就會被添加到workspace。
如果按照以上描述進行預處理,數(shù)據(jù)將具有以下字段:
data =
? ? ? ?hdr: [1x1 struct]
???? ?label:?{48x1?cell}
????? ?time:?{[1x4462?double]}
? ? ?trial: {[48x4462 double]}
?? ?fsample:?5
sampleinfo:?[1?4462]
? ? ? opto: [1x1 struct]
? ? ? ?cfg: [1x1 struct]
“hdr”字段包含關于數(shù)據(jù)的所有最高級信息,例如原始采樣率、通道數(shù)量等。字段“l(fā)abel”列出了你決定讀入的所有通道的名稱。字段“time”,表示數(shù)據(jù)集的時間軸?!皌rial”字段包含所有通道的數(shù)據(jù),它被稱為“試次”,因為通常情況下,數(shù)據(jù)被分割成與實驗中的不同試次相對應的片段。字段“fsample”描述了“trial”字段中數(shù)據(jù)的采樣率?!皊ampleinfo”字段描述了相對于磁盤上原始測量值的每個試次的樣本編號?!皁pto”包含了有關通道和光電極組成的信息,例如使用了什么波長,光點極位置等等。最后,“cfg”字段與剛才使用的cfg相同,只是通過一些默認值進行了擴展。
為了快速查看數(shù)據(jù),可以使用ft_databrowser函數(shù)。databrowser不僅僅是一個簡單的“數(shù)據(jù)瀏覽器”,不過現(xiàn)在將利用這個功能來實現(xiàn)查看數(shù)據(jù)的目的。在配置中添加'ylim ='maxmin'來調(diào)整y軸,使圖中的最小值決定y軸上的最低點,圖中的最大值決定y軸上的最高點。
cfg = [];
cfg.ylim = 'maxabs';
ft_databrowser(cfg, data);

使用ft_databrowser,還可以刪除不需要的數(shù)據(jù)片段。例如,如果你在放置光電極的同時開始記錄,可能會在記錄的開始有一大段不需要分析的數(shù)據(jù),這些數(shù)據(jù)包含與大腦無關的值,并且快速波動。把這些片段剪掉更有利于后續(xù)的分析。本示例數(shù)據(jù)不需要裁剪。
此外,這里先只選擇一對通道,以減少分析過程的復雜性。
cfg = [];
cfg.ylim = 'maxmin';
cfg.channel = {'Rx4b-Tx5 [860nm]', 'Rx4b-Tx5 [764nm]'}; % you can also use wildcards like 'Rx4b-Tx5*'
ft_databrowser(cfg, data);

去除偽影近紅外光譜數(shù)據(jù)不僅包含了氧和脫氧血紅蛋白濃度的變化,而且還包含了其他生理信號、隨機噪聲和來自測量環(huán)境的變化。因此,可能存在較為明顯的運動偽影,它是由光電二極管和頭皮之間接觸的臨時變化產(chǎn)生的,通常是由頭部運動引起。這些運動偽影可能存在于所有通道中,但也可能只有一個或幾個通道受到影響。比如,數(shù)據(jù)中短的、意想不到的峰值被認為由運動導致,可以通過ft_artifact_zvalue檢測和刪除這些偽影。
指定一些預處理選項來檢測偽影,如下所示:
cfg?=?[];
cfg.artfctdef.zvalue.channel?=?{'Rx4b-Tx5?[860nm]',?'Rx4b-Tx5?[764nm]'};
cfg.artfctdef.zvalue.cutoff?=?5;
cfg.artfctdef.zvalue.hpfilter?=?'yes';
cfg.artfctdef.zvalue.hpfreq?=?0.1;
cfg.artfctdef.zvalue.rectify?=?'yes';
cfg.artfctdef.zvalue.artpadding?=?2;
%?cfg.artfctdef.zvalue.interactive?=?'yes';?%?the?interactive?display?makes?more?sense?after?segmentating?data?in?trials
[cfg, artifact] = ft_artifact_zvalue(cfg, data);
運行上述代碼可以發(fā)現(xiàn),F(xiàn)ieldTrip在此過程中識別了8個偽影。請注意,只是標識出來了這些偽影,還沒有被刪除。將其緩存,并在數(shù)據(jù)濾波分段之后調(diào)用ft_rejectartifact。
轉(zhuǎn)化為oxyHB/deoxyHB濃度的變化
此時的數(shù)據(jù)是光密度(OD)值,而不是含氧和脫氧血紅蛋白濃度??梢允褂胒t_nirs_transform_ODs將數(shù)據(jù)轉(zhuǎn)換為濃度變化。在使用ft_nirs_transform_ODs時,要做出的一項選擇是差分路徑長度因子(DPF),它取決于參與者的年齡和組織類型(例如,當分析的是肌肉組織而不是大腦的血氧變化時)。
cfg = [];
cfg.dpf = 5.9;
cfg.channel = {'Rx4b-Tx5 [860nm]', 'Rx4b-Tx5 [764nm]'};
data_conc = ft_nirs_transform_ODs(cfg, data);
將功能響應與系統(tǒng)反應分開
除了出現(xiàn)在數(shù)據(jù)中的運動偽影,fNIRS測量還包含系統(tǒng)反應,這是研究中不感興趣的,并希望從數(shù)據(jù)中能夠過濾掉這些反應。有多種方法可以提取功能響應(大腦的血流動力學響應),這里將采用一種簡單但相對有效的方法,即應用帶通濾波器來消除頻帶以外的信號值,該頻帶捕獲了大腦血流動力學響應中的動態(tài)。由于實驗任務和血流動力學響應函數(shù),這里預計活動低于0.1Hz(對應~10s的波動)。此外,將較低的范圍設置為0.01Hz(波動~100s),以消除緩慢漂移。
cfg = [];
cfg.bpfilter = 'yes';
cfg.bpfreq = [0.01 0.1];
data_filtered = ft_preprocessing(cfg, data_conc);
定義興趣段
如果想知道大腦是否會對實驗中發(fā)生的事件做出特定的反應,所以需要把分析的重點放在事件發(fā)生的時間窗上(這些時間窗通常被稱為epoch或trials)。
此前,通過調(diào)用ft_preprocessing時不指定cfg.trl,讀取了所有epoch。在幫助文檔中,cfg.trl指向ft_definetrial,現(xiàn)在將使用它來定義試次。接下來先看一下函數(shù):
help ft_definetrial
假設我們對數(shù)據(jù)中的triggers沒有任何線索。因此,現(xiàn)在將利用隱藏的功能來檢索這些信息。
cfg?=?[];
cfg.dataset?=?'motor_cortex.oxy3';
cfg.trialdef.eventtype = '?';
問號表示不確定事件triggers,因此函數(shù)ft_trialfun_general將輸出在數(shù)據(jù)集中找到的所有事件。接下來,調(diào)用ft_definetrial
ft_definetrial(cfg);
在命令窗口中,將看到發(fā)現(xiàn)了24個事件,并且使用了兩種類型的事件,即事件'A'和事件'B'。此外,需要注意的是,現(xiàn)在沒有定義輸出變量'cfg',因此變量'cfg'將保持不變。如果我們想看看從刺激前10s(事件'A'之前)開始到'A'之后35s(刺激后)結(jié)束的epoch/trials??梢远x試次,然后使用它來調(diào)用ft_preprocessing:
cfg.trialdef.eventtype = 'event';
cfg.trialdef.eventvalue = 'A';
cfg.trialdef.prestim = 10;
cfg.trialdef.poststim = 35;
cfg = ft_definetrial(cfg);
cfg.channel = {'Rx4b-Tx5 [860nm]', 'Rx4b-Tx5 [764nm]'};
data_epoch = ft_redefinetrial(cfg, data_filtered);
現(xiàn)在已經(jīng)選擇了一對通道,并在12個試次中剪切了數(shù)據(jù)。使用databrowser進行檢查,使用以下這些設置,可以使繪制圖形看起來更整潔,同時也可視化之前被識別出的偽影:
cfg = [];
cfg.ylim = [-1 1];
cfg.viewmode = 'vertical';
cfg.artfctdef.zvalue.artifact = artifact;
ft_databrowser(cfg, data_epoch);

現(xiàn)在可以刪除確定為偽影的試次。
cfg = [];
cfg.artfctdef.zvalue.artifact = artifact;
cfg.artfctdef.reject = 'complete';
data_epoch = ft_rejectartifact(cfg, data_epoch);
鎖時分析
使用ft_timelockanalysis計算試次的平均值。隨后,可以使用FieldTrip繪圖函數(shù)來查看數(shù)據(jù),也可以簡單地使用MATLAB繪圖函數(shù)進行可視化。
cfg = [];
data_timelock = ft_timelockanalysis(cfg, data_epoch);
輸出是帶有以下字段的data_timelock數(shù)據(jù)結(jié)構
data_timelock =
? ? ?avg: [2x225 double]
? ? ?var: [2x225 double]
? ? time: [1x225 double]
? ? ?dof: [2x225 double]
? ?label: {2x1 cell}
? dimord: 'chan_time'
? ? opto: [1x1 struct]
? ? ?cfg: [1x1 struct]
接下來,繪制A-10s到A+35s的平均O2Hb和HHb軌跡。根據(jù)fNIRS慣例,O2Hb為紅色,HHb為藍色。
time = data_timelock.time;
O2Hb = data_timelock.avg(1,:);
HHb = data_timelock.avg(2,:);
figure;
plot(time,O2Hb,'r'); hold on;
plot(time,HHb,'b');
legend('O2Hb','HHb'); ylabel('\DeltaHb (\muM)'); xlabel('time (s)');

參考來源:Preprocessing and averaging of single-channel NIRS data.https://www.fieldtriptoolbox.org/tutorial/nirs_singlechannel/https://www.artinis.com/blogpost-all/2017/6/27/how-do-i-choose-the-correct-fibers-and-template-for-my-oxymon
https://download.fieldtriptoolbox.org/tutorial/nirs_multichannel/
小伙伴們點星標關注茗創(chuàng)科技,將第一時間收到精彩內(nèi)容推送哦~
