MATLAB の基礎とスパイクデータ解析 MATLAB tutorial: Analysis of Extracelluar recordings この演習では、モデルによって生成されたスパイクデータをどのようにしてMATLABで取り扱うのかについて説明します。 MATLABでは行列計算と微分方程式が扱えるところまでお話しましたが、ここでは特に簡単な統計解析についてみていきましょう。 その前に、ここで扱うMATLABの基本事項をいくつか確認しておきましょう 1. 基本事項 MATLAB は基本的に行列を１つの変数として扱う言語ですが、スカラーはもちろん行ベクトルや列ベクトルを扱えます。 スカラーの場合 >> s = 3 行ベクトルの場合 >> v = [8 3 5] 列ベクトルの場合 >> u = [4.2; 5; 7] のようにセミコロンで分割します。 3かけ3行列をこの書式でかいてみれば >> A = [3 4 0; 1 4 6; 9 7 2] のようになりますね。 さてこのような行列を転置 $A^t(i,j) = A(j,i)$ する場合 >> A' また、これを使えば行ベクトルを列ベクトルに、またその逆に変換できます。 >> v' これらの変数から新たな3かけ3行列を作ってみましょう >> B = [v; u'; 1 0 4]; これで２つの行列ができました。これらの和 >> A+B 積 >> A*B が計算できます。 ここで気をつけなければならないのは、行列演算であるということです。とくに積に関しては、行列の積のルールに従っていますから、行列の次元に気をつける必要があります。また一方、行列積ではなく、同じサイズの行列の同一要素同士の積 ($A(i,j)B(i,j) = answer(i,j)$)を行いたい場合もあります。このような場合には、 >> A.*B とドット(.)を用いると同一要素間の演算が行えます。このドットは他の演算たとえば累乗 >> A.^2 などにも使えます。これは各要素の2乗をおこなうのであってもちろん A^2 とは違います。和や差にもつかえますが、これはつけなくても同じですね。 MATLABでは特別な行列に対して行列生成関数を用意しています。 たとえば >> ones(3,4) は、すべての要素が1になるような3行4列の行列を作ります。とすると >> zeros(5,7) はどうなるかわかりますね。5次元（5かけ5) 単位行列は >> eye(5) として作ることができます。 等間隔の要素を持つベクトルを作りたい時には、 >> x = 1:10 と打ち込めば1から10まで1刻みのベクトルが作られるます。 刻み幅を変えたい場合には、 >> y = 0:0.1:8 とすれば0から8まで0.1刻みでベクトルを作ってくれます。 このようにMATLABで変数を作ればこれはworkspace上に保管されています。 どんな変数を作ったのかを確かめたいときには >> who または >> whos を使います。whosでは変数のサイズやクラスも表示してくれます。 使い方やMATLABに組み込まれている関数やコマンドについての情報を得たいときには >> help と打ち込めばテキストベースのヘルプを見ることができます。 また特定のコマンドについての情報はhelp に続けて関数名やコマンド名を書きます。 たとえば、 >> help rand を試してみてください。ちょっと長いですが、[0 1]間の一様乱数生成関数randの使い方が表示されます。他にも、関係するトピックを探したい場合には、lookforを使います。 >> lookfor complex しばらくサーチに時間がかかりますが、複素数に関連する関数を表示してくれます。 他にも >> helpdesk と打てばブラウザーが立ち上がってhypertextのヘルプを見ることができます。 キーワード検索や関数検索、GUIオブジェクトのヘルプや特定のtoolboxのヘルプなど がほしい場合に使います。 1-2. 添え字 行列やベクトルの特定のようをを取り出すときには、添え字を使いますベクトルの場合には次元は１なので行ベクトルでも列ベクトルでもi番目の要素を取り出すときには >> v(i) とすれば得られます。 行列の場合には`A(i,j)'と打てば i行 j列の要素を取り出すことができます >> A(2,3) ならAの2行3列の１要素だけを見る事ができます。 行列から特定の行や列を取り出したい場合には`:'を使います。たとえば >> A(1,:) なら１行目の全要素を含む行ベクトルが >> A(:,2) なら２列目の全要素を含む列ベクトルが取り出せます。同様に >> A(:,1:2) とすると１から２列目の部分行列を取り出すことができます。 行列のサイズはわからないけれどとにかく最終行の行ベクトルが必要な場合 >> A(end,:) でとりだせます。 これらは特定の部分だけの総和が計算したい場合などに便利です。 たとえば、３行目の和を求めたければ >> sum(A(3,:)) ですし、３列目の全要素の和なら >> sum(A(:,3)) です。１部分の要素だけを特定の値にセットしたい場合もあります。たとえば >> B(1:2,2:3) = 0 をためしてみてください。どの部分が0になりますか？ 1-3. 論理表現とデータ抽出 >> A>3 と打ってみてください。Aの各要素が3より大きいところに1が立ち、 それ以外は0の行列が表示されたと思います。これを論理表現とよびます。 さて、この行列を使うと要素のうち必要な部分だけに演算を加えたり操作し たりするのが簡単にできます。次を試してみましょう。 >> (A>3).* A どうですか、結果はどうなるかわかりますか？ ちょっと脱線 ここまでは演算や行列を作るごとにコンソールに結果が表示されてきました。 表示させないようにするのには文末に`;’をつけます。`;'から先は別の文として 処理されます。 >> x=3; y=[1 x ones(1,5); 5*rand(1,7); 1:7]; と打っても結果は表示されませんね。 >> y とすればyに何が入っているのか表示されます。 先ほどは論理表現を使った行列との演算を示しましたが、論理表現を使って要素を取り出すことができます。 >> L = A>3 とするとLにAの要素が3以上である要素が１となるAと同じサイズの論理表現行列 ができます。これを使って >> A(L) とするとどうなるでしょうか。Lが1である要素だけを取り出した列ベクトルを生成 することができます。 MATLABでは特別な数でない要素'NaN'と'Inf' を扱うことができます。測定器の不調でデータを見失ったり、演算の結果0割りで無限大になったりした場合にこのラベルを使うことができます。 >> 0/0 >> 1/0 >> -1/0 をやって結果がどうなるのか見てください。 たとえば, >> x = [ 2.1 1.7 1.6 1.5 NaN 1.9 1.8 1.5 5.1 NaN 1.8 1.4 2.2 1.6 1.8] というデータがあったとします。ここから観測できなかったデータを取り除きましょう。それには特別な論理関数を使います。データが有限である(NaN やInfでない） 論理表現を出力する関数'finite'を使います。 >> x = x(finite(x)) さて、5.1はどうも外れ値のようですね。このような外れ値を取り除く処理はどのように したら良いでしょうか、以下を試してください。 >> x = x(abs(x-mean(x)) < 3*std(x)) abs は符号を取り除く絶対値関数、`mean', `std'は平均と標準偏差を求める 統計関数です。そのほかMATLABでは基本的な統計処理のための関数をいくつか 用意しています統計関数に関しては、`help datafun' を参照してください。 うまく使えばこのようにデータの偏差が標準偏差の３倍以下の要素だけを新たに xに入れるという操作が１行で表現できます。 'find' 関数もこういった要素の抽出に便利です。特に要素のインデックスを 取り出すことができます。くわしくは >> help find 参照のこと。 1-4. 可視化　(Visualization) MATLAB では非常に強力な描画コマンドが用意されています。もっとも簡単なものの１つは'plot' 関数です。使い方はいたって簡単、[-5 , 5] の間のsin関数を描きたければ、 >> x = [-5:0.1:5]; >> y = sin(x); >> plot(x, y); のように描画したい点列のx軸、y軸をベクトルとして渡してやれば描いてくれます。 ここでは、並んでいる点間を線で結んでいますが、データ点のみを書きたければ >> plot(x,y,'.'); とすれば点のみを書いてくれます。 もちろん xは単調増加する点列でなくてもあつかえるので、たとえば >> plot([0 1 0 1], [0 0 1 1], '.') なら　正方形の四隅に点を打ったグラフを表示します。 plotの３つめの入力に指定されている'.'はプロットシンボルの種類を指定してます。 >> plot(x,y,'o'); を試してみましょう。他にも色や点間の線の種類なども指定できます。その種類については >> help plot 参照のこと。 plot のような２次元のグラフのみでなく、３次元のデータも簡単に書くことができます。もし使いたい場合は >> help graph2d >> help graph3d 参照のこと。 2. データの読み込み 道のりは長かったですが、やっと基礎の説明が終わったところで、神経細胞活動の解析に移りましょう。 神経活動データを録っているシステムによって記録形式は様々です。統一形式でも あれば問題ないんでしょうけれども、なかなか録る側になってみると、リアルタイムで 取りこぼし無くデータを簡単に取れるようにするということも非常に重要なことなので 他人に合わせるというとこまでなかなか行かないのではないでしょうか。 ここでは、京都府立医大の南本氏が開発したシステムで採られたデータをサンプルとしてもっとも基本的なヒストグラムを表示するプログラムをMATLABを使って作ってみます。 スパイクタイミングのデータはイベントコードシステムを使って記述されています。イベントコードとは、実験中におきたタスクのイベント、LEDの点灯消灯や行動のタイミングや報酬などを特定のイベントコードを割り当てて記録しています。もちろんスパイクも１つのイベントとしてイベントコードが割り振られます。 このデータを取った時のタスクは、パネル上に3つボタン付のLEDがあり、点灯したLEDを順に３つ押すと報酬信号（ビープ）が出て報酬（水）がもらえるというタスクになっています。点灯するLEDの順番は、まず下のLED（LED1）が点灯し、時計回り(CW)と反時計回り(CCW)の２種類があります。よって CCWではLED1→LED3→LED2となります． ○LED2 ○LED3 ○LED1 このタスクのイベントコードは以下のようになっています。 ----------------------イベントコード表---------------------------------- 1, 2 Neural Spikes 13 REW Signal (CW) 14 REW Signal (CCW) 17 REW 18 REW off 21 Button1 Push 22 Button1 Release 23 Button2 Push 24 Button2 Release 25 Button3 Push 26 Button3 Release 64 LED 1 ON(CW) 65 LED 1 OFF(CW) 66 LED 2 ON(CW) 67 LED 2 OFF(CW) 68 LED 3 ON(CW) 69 LED 3 OFF(CW) 76 LED 1 ON(CCW) 77 LED 1 OFF(CCW) 78 LED 2 ON(CCW) 79 LED 2 OFF(CCW) 80 LED 3 ON(CCW) 81 LED 3 OFF(CCW) ----------------------イベントコード表------------------------------ データファイルはたいていバイナリーファイルや、テキストファイルで保存されています。この場合にはバイナリーファイルとして保存されており、最初の100word にヘッダ情報が書き込まれており、続いてイベントコード、時刻が並んで記録されています。 このデータを workspace に読み込んでみましょう。 最初の100wordのヘッダファイルのフォーマットには --------------------------------------------------------------------- データファイルのバイナリ構造ですが， 初めの100個（Short)がヘッダーになっています． ヘッダーの中身は [1] ; 次のファイルの有無（0;無, 1;有） [2] ; ファイル名 [3] ; extention(1,2,3,.....)1 [10] ; DI data数+４ [11] ; A/Dの測定開始時のカウンタ となっています．（ほかにもいくつか入ってますがとりあえず必要な情報です) このあとイベントコードがDI Data数だけ続き，終了コード（-1）4つを挟んで 時刻(単位msec)がDI Data数続きます． --------------------------------------------------------------------- とありました。 まず最初の100wordのデータを読み込みましょう。 ここでは、サンプルデータとしてファイル`2244.1'を読み込んでみます。 MATLABでデータを読み込むには、指定するファイルを開く fopen 関数 バイナリーデータを、指定した個数だけ読み込む fread が使えます。 >> fid = fopen('2244.1','r','l'); >> Header =fread(fid,100,'int16'); これでヘッダの情報が読み込めました。次にイベントコードを読み込みますが、 その前にイベントコードのデータ数は >> EventNum = Header(10); 次にイベントコード >> Event = fread(fid,EventNum,'int16'); さらに時刻 >> Time = fread(fid,EventNum,'int16'); を読み込みます。 Timeはそのイベント発生時刻をmsecであらわします。ただし、このTime は16bitカウンタの値そのままなので、32768msecを超えると-32767に戻ってしまいます。 そこで、まず前のイベントからのイベント間隔 TimeDif をこのデータから計算します。 それを元に絶対時間を求めたデータを作成します。 TimeDifはTimeデータの直前の要素間の差によってこのデータは作れますが、直前との差が負になってしまうときだけ例外的な処理を入れます。この場合、MATLABのif 文とfor文を使って >>TimeDif(1) = 0; >>for i = 2:EventNum-1 >> if Time(i) - Time(i -1) < 0 >> TimeDif(i) = 2^16 + Time(i) - Time(i - 1); >> else >> TimeDif(i) = Time(i) - Time(i-1); >> end >>end とすれば　TimeDif に時間差データを作ることができます。 同じことをfor文を使わずにやろうとすれば >> TimeDif = [0; Time(2:end)-Time(1:end-1)]; >> TimeDif(TimeDif<0) = TimeDif(TimeDif<0) + 2^16; となります。 次に、Event と TimeDif を使って相対時間データ行列 datadif を用意します。 >> datadif = [Event(1:end-4) TimeDif(1:end-4)]; Event TimeDifの最後の４個のデータは終了コードが入っているため、ここでは使いません。 さて、TimeDifは相対時間ですので、絶対時間に直しましょう。そのためには、相対時間を順に和をとっていくようにすれば良いので >> data = [Event(1:end-4) cumsum(TimeDif(1:end-4))]; これで、データを取り始めてからの絶対時間になります。 安全のためにできあがったデータ行列を保存しておきましょう。いつ電源が落ちる、OSがハングアップする。Math works に連絡しろというメッセージを出してMATLABの動作がおかしくなるかもしれません。 workspace 全体を保存するには単純に >> save とすれば、現在のディレクトリにmatlab.mat というファイルを作成しデータをそこにストアします。 特定のファイル名を指定したいときは >> save ファイル名 特定の変数行列のみを保存したいときは >> save ファイル名 変数名 で保存できます。 >> save '2244.mat' data とすれば出来上がったデータを保存することができます。 3. 特定のイベントにあわせたラスターの表示 興味のあるイベントの周りでのスパイクをラスター表示してみましょう。 ここでは、LED1が点灯したときのラスターを作ってみます。 まず、LED1が点灯した時間を求めます。反時計回りでLED1が点灯した時間を 抽出するには、イベントコード64の時刻を知ればよいので、dataの1列目が64 であるところの２列目 >> align = data(data(:,1)==64,2) で知ることができます。LED1の点灯時刻の１秒前から1秒後までのスパイクタイミングを知るためには、まずスパイクをデータから取りだしておきましょう。スバイクのイベントコードは１または２ですから >> allspikes = data(data(:,1) ==1 | data(:,1)==2 ,2) '|' は or 演算子です。１回目のLED1の点灯前後のスパイクの相対時刻を計算しましょう。 >> spikes = allspikes-align(1) 必要となるのは-1000msecから1000msecのデータだけですから >> spikes = spikes(spikes > -1000 & spikes < 1000) これをすべてのLED1点灯時について計算します。 それにはMATLABでのforループを使って >> buf = []; >> for i = 1:length(align) >> spikes = allspikes-align(i); >> spikes = spikes(spikes > -1000 & spikes < 1000); >> buf=[buf ; ones(size(spikes))*i spikes]; >> end buf はすべてのスパイクの相対時間をためておくバッファだと思ってください。 buf の1列目には何番目のイベント（ここでは試行数）を示すインデックスが入っていて ２列目にはそのイベントからの相対時間(msec)が入ります。 データが用意できたので描画に移りましょう。plot関数を使ってラスタープロットを 描いて見ましょうスパイクの相対時間を横軸に試行数を縦軸にプロットすると >> plot(buf(:,2), buf(:,1), '.') となりますね。 4. PSTHの作成 ラスターが取れるようになれば、同じデータからヒストグラムを作ることも簡単に できます。buf にためられたすべてのスパイクについてのヒストグラムを作って見ましょう。MATLAB では簡単にヒストグラムを作る関数 'hist'が用意されています。 >> hist(buf(:,2)) これでは自動的にスケールとビンを決められてしまうので、 >> hist(buf(:,2), [-995:10:995]) とすれば10msec の幅でのヒストグラムを作ることができます。ただし、この場合は 縦軸はスパイク数になります。縦軸をスバイク数ではなく平均発火頻度にしたい場合 には、histで直接ヒストグラムを描画してしまわずに一旦スパイク数のデータを取り出して計算する必要があります。histを使ってデータを抽出し >> [spcount, bins] = hist(buf(:,2),[-995:10:995]); >> sprate = spcount/(10*10^-3)/length(align) ヒストグラムのような棒グラフをデータから書くには >> bar(bins, sprate,1); とすればLED1の点灯を0とするスパイクヒストグラムを書くことができます。 これでヒストグラムとラスターの両方を描くことができるようになりましたが、それぞれ独立に見るのではなく２つの情報を１つのfigureの中に収めたいですね。上にヒストグラム、下にラスター表示すると見やすそうです。このような場合には subplot を使います。１つのfigure の中に2つのaxisと呼ばれるグラフ領域を作成し、それぞれを描画します。 >> subplot(2,1,1) >> plot(buf(:,2), buf(:,1), '.') >> subplot(2,1,2) >> bar(bins, sprate,1); のようにグラフを描画する関数の直前にsubplotを入れ、特定のaxisを指定します。 subplot(m,n,i)は１つのfigureに　m行n列の axis を用意し、i番目のaxis に描画し なさいという意味です。 5. m-file スクリプトと関数 さてここまでのルーチンを１つの流れとして一発のコマンドでやってのけたいという場合には、これまで書いてきた作業を１つのファイルに書き込むことによってできます。 Windows 版であればMATLAB main window から　ファイル -> 新規作成 -> M-file とすることによってm-file editorを立ち上げることができますし、既存のエディタを 使ってm-fileを作ることができます。 ここではデータファイルを読み込んだ後の作業を１つのm-file Psth.mとして作ります。 ----------------------------Psth.m----ここから-------------------------- % % Psth.m % コメントは% を使います、各行の% 以降 はコメントとして無視されます。 % % またスクリプトや関数の最初の行からのコメントは %　>> help 関数名 % を実行したときに表示されます。 % ここに関数の使い方などを書き込んでおけば、他人や未来の自分（しばしば作 % った本人が関数の使い方を忘れます）がスクリプトを直接読まずに使い方がわ % かって便利です。 % % Psth.m : dataからラスター、PSTHを描画する。 % % Input : data ; column1- event code, column2 - Time(msec) % align = data(data(:,1)==64,2); allspikes = data(data(:,1) ==1|data(:,1)==2 ,2) buf = []; for i = 1:length(align) spikes = allspikes-align(i); spikes = spikes(spikes > -1000 & spikes < 1000); buf=[buf ; ones(size(spikes))*i spikes]; end [spcount, bins] = hist(buf(:,2),[-995:10:995]); sprate = spcount/(10*10^-3)/length(align); subplot(2,1,1) plot(buf(:,2), buf(:,1), '.'); subplot(2,1,2) bar(bins, sprate,1); ---------------------------ここまで------------------------------------ このm-file を使えばある任意のイベントに対するヒストグラムが作れるはずです。 10行目を少し改良して、 align = data(data(:,1) == aevent,2); として >> aevent = 64 >> Psth さらに、LED1の点灯は反時計回り(64)だけでなく時計回り(76)でも点灯します。 複数のイベントコードに対応するように改良します。 10行目をさらに４行増やして、 align = []; for e = aevent align = [align; data(data(:,1) == e, 2)]; end として >> aevent = [64 76]; >> Psth を試してみてください。 さて、このようなスクリプトとしてm-file を使うだけでなくある変数を渡して ヒストグラムを描くような関数を作ることもできます。それには、m-fileの1行目に function sprate=Psth(data,aevent) というのを書き加えるだけです。 これで >> Psth(data, [66 78]) とすれば、イベントコード66と78にそろえたラスターとヒストグラムを作成してくれるはずです。 スクリプトと関数の違いは、現在のワークスペースを使うかどうかです。関数の場合は、現在のワークスペース内にある変数や行列を直接参照したり、変更することはできません。そのために入力と出力を特定してやる必要があるのです。 演習 1. サンプルデータ2244.1は何に反応しているといえるでしょうか？ さまざまなイベントにそろえたPSTHで検討してみてください。 2. 作ったPsth.m をもっと一般的に使えるようにしたくなりました。 2-1. ここでは1秒前から1秒後までとしましたが、もっと一般的にstartw から stopw までとするにはどう改良したらよいでしょうか？ 2-2 bin幅を10msec ではなく可変にするには？