生命科学研究のためのmatlab

生命科学研究のためのmatlab ～2次元で考えるシリーズ～

一次元時系列データ (脳波、筋電など)

データを読み込む

function eeg = readmulti(fname,numchannel)

datafile = fopen(fname,'r');

eeg = fread(datafile,[numchannel,inf],'int16');

eeg = eeg';

end

トリガポイントを探す

% function Times = TriggerPoints(data,threshold,minISI)

% returns triggered points with a inter-trigger-interval of at least minISI samples

% triggered points are defined as points that exceeds 'threshold'

function Times = TriggerPoints(data,threshold,minISI)

data = data(:);

tempTimes = find(data>threshold);

if isempty(tempTimes)

Times = tempTimes;

else

Times = tempTimes([true;diff(tempTimes)>minISI]);

end

例えば1chにデータ、2chにトリガを記録したデータ eeg(time, channel) があるとしたら、data = eeg (:,2)とすれば良い。

加算平均を計算する

トリガポイントが探せたら、トリガポイント周辺でデータを切って二次元的に結合しておけば、後で便利。例えば、

eeg = readmulti(fname,2);

index = TriggerPoints(eeg(:,2),threshold,minISI);

for dd = 1:size(index,1)./NumTrials

AMP=[];

for ii = (1+(dd-1).*NumTrials:dd.*NumTrials)

AMP = cat(2,AMP,eeg(index(ii)-Tstart.*SamplingRate:index(ii)+Tend.*SamplingRate,1));

end

fileAMP = strcat(SaveDir,'/AMP',num2str(dd),'.mat');

save(fileAMP,'AMP');

end

ここまでやっておくと加算平均は一行、

Y = mean(AMP(:,Trials),2)

Evoked potentialの振幅を計算する

function [YAMP,YIndex]=CalcAMP(data,start,from,to,N)

%1:fromまでのデータの平均を最大値Ymax,

%from:toの間での最小値の周り±Nフレームの平均値をYminとし

%YmaxとYminの差からAmplitudeを計算。

Ymax = mean(data(start:from,:));

[~,index] = min(data(from:to,:));

Ymin = mean(data(from+index-N:from+index+N,:));

YAMP = Ymax - Ymin;

YIndex = from+index;

end

Evoked potentialのSlopeを計算する

function slope = CalcSlope(data,from,to,R_t1,R_t2,)

[Ymin,Imin] = min(data(from:to,:),[],1);

for ii = 1:size(data,2)

[Ymax(ii),Imax(ii)] = max(data(from:(from-1+Imin(ii)),ii));

end

Thr_t1 = Ymax+R_t1*(Ymin-Ymax);

Thr_t2 = Ymax+R_t2*(Ymin-Ymax);

for ii = 1:size(data,2)

t2(ii) = (from-1) + find(data(from:(from-1+Imin(ii)),ii)>Thr_t2(ii), 1, 'last' ); %decide minimum index;

t1(ii) = (from-1) + find(data(from:t2(ii),ii)>=Thr_t1(ii), 1, 'last' );

x = t1(ii):t2(ii);

YYY = data(t1(ii):t2(ii),ii);

XXX = [x',ones(length(x),1)];

AAA = XXX \ YYY;

slope(ii) = AAA(1);

P(ii,:) = AAA;

X{ii} = x;

end

波形全体の平均的な振幅と周波数を数える

function [amp,freq] = EventDetector(FFF,XX)

% FFFはデータでXXは解析したいチャネルの番号。

% 閾値を設定して、閾値を超えた部分の最大ピーク値と閾値を下回る最小ピーク値の全体の平均値によって、平均的な振幅を計算する。

% ピーク値が時間内に何回出現したかによって大まかな周波数を計算する。

xx = 0.5; % 閾値の値を決める要素

Thr1 = mean(FFF(:,:),1) + xx.*std(FFF(:,:));

Thr2 = mean(FFF(:,:),1) - xx.*std(FFF(:,:));

index1 = find(FFF(:,XX)>Thr1(:,XX));

ii1 = find(diff(index1)>1);

Time = min(index1):max(index1);

index2 = find(FFF(:,XX)<Thr2(:,XX));

ii2 = find(diff(index2)>1);

Xmax = []; Ymax = []; Xmin = []; Ymin = [];

for jj = 1:length(ii1)-1

[ymax,Maxindex] = max( FFF(index1(ii1(jj)):index1(ii1(jj+1)),XX) );

Xmax = cat(1,Xmax,Maxindex+index1(ii1(jj))-1);

Ymax = cat(1,Ymax,ymax);

end

for jj = 1:length(ii2)-1

[ymin,Minindex] = min( FFF(index2(ii2(jj)):index2(ii2(jj+1)),XX) );

Xmin = cat(1,Xmin,Minindex+index2(ii2(jj))-1);

Ymin = cat(1,Ymin,ymin);

end

T = 1./10:0.1:size(FFF,1)./10;

Thr = Thr1(:,XX)*ones(1,size(FFF,1));

figure;plot(T,FFF(:,XX),T,Thr','r',T(index1(ii1)),FFF(index1(ii1),XX),'or',T(Xmax),FFF(Xmax,XX),'xg');

figure;plot(T,FFF(:,XX),T,Thr','r',T(index2(ii2)),FFF(index2(ii2),XX),'or',T(Xmin),FFF(Xmin,XX),'xg');

amp = mean(Ymax)-mean(Ymin);

freq = length(Xmax)./(length(Time)./10);

end

例えば、自発的な脳波振動の性質をざっくりと知りたい場合はこの方法が使える。この方法は、1次元化した後の画像データの時系列解析にも適用できる。

周波数 (スペクトル)解析

data = readsingle(DataFileName);

% power spectrum plotting (old??)

Hs=spectrum.welch('Hamming',2^14,50);

hpsd = psd(Hs,data,'NFFT',2^16,'Fs',SamplingRate);

figure;

plot(hpsd.Frequencies,10*log10(abs(hpsd.Data)));

grid on

xlabel('Frequency (Hz)');

ylabel('Power Spectrum Magnitude (dB)');

axis([0 45 0 70])

HH = [hpsd.Frequencies,10*log10(abs(hpsd.Data))];

% spectrogram (current form)

% doc spectrogram, for more code examples

Window = 2^15;

[S,F,T,P] = spectrogram(data,Window,Window/2,Window*2,SamplingRate);

figure;

PP = 10*log10(P);

imagesc(T,F,PP);

axis xy

axis([0 700 0 45])

caxis([-50 10.*log10(max(P(:)))])

二次元時系列データ (画像)

例えば、512 x 512 pixelのTifフォーマットの画像列をフレームレート = SamplingRateで N枚読み込んで保存してあるとしよう。

一回一回ImageJかなんかで読み込んで解析しても良いのだが、ここではmatlabで取り扱う方法を紹介する。

さらにXYTの三次元情報を二次元に次元縮約して取り扱う、「二次元縛り」で解析してみようと思う。

画像を読み込む

img = imread(file,'tif');

ビニングする

http://mathforum.org/kb/message.jspa?messageID=6884726 参照

例えば、画像をNumBinning分の1にビニングしたいとする。

つまり周辺NumBinningのピクセルを平均して一つのピクセルと考える。

M = img; p = NumBinning; q = NumBinning;

[m,n] = size(M); %M is the original matrix

M = sum( reshape(M,p,[]) ,1 );

M = reshape(M,m/p,[]).'; %Note transpose

M = sum( reshape(M,q,[]) ,1);

M = reshape(M,n/q,[]).'; %Note transpose

img2 = M;

二次元化する

function img3 = Image2Data(file_start,file_end,NumBinning)

for fff = file_start:file_end

fprintf('%d\n',fff);

img = imread(file,'tif');

%http://mathforum.org/kb/message.jspa?messageID=6884726

%Example for pxq binning

M = str; p = NumBinning; q = NumBinning;

[m,n] = size(M); %M is the original matrix

M = sum( reshape(M,p,[]) ,1 );

M = reshape(M,m/p,[]).'; %Note transpose

M = sum( reshape(M,q,[]) ,1);

M = reshape(M,n/q,[]).'; %Note transpose

img2 = M;

img3 = cat(2,img3,img2(:));

end

これで、撮ったN枚の画像は全て512x512行N列の行列に変換された。

今後はこの作製した行列の演算を行えばいいということになる。

ファイルの書き出し/読み込み

一回一回画像を読み込むのは時間がかかるので、上で二次元化した行列は名前をつけてsaveしておく。

filename = strcat(SaveDir,'/img3.mat');

save(filename,'img3');

行列化した画像列がどういう内容だったか簡単に記述しておくために、画像全体で平均化して時間方向にプロットした画像も同じフォルダに保存しておく。

mimg3 = mean(img3,1);

file2 = strcat(SaveDir,'/Event.tif');

h = figure;

plot(mimg3);

saveas(h,file2,'tif');

次回からは、この.matファイルを読み込んで使用する。

img3 = strcat(ReadDir,'/img3.mat');

importfile(img3);

蛍光輝度変化率 (ΔF/F)を計算する

撮った画像がCa²⁺イメージングなどの蛍光画像だとしよう。

蛍光画像の評価には、基準値からの変化の割合 (蛍光輝度変化率)用いる。

蛍光輝度から蛍光輝度変化率への変換は ΔF/F₀ = (F - F₀)/F₀ で与えられる。

ここで、F₀は平均輝度なのだが、平均は任意の時間で設定する。例えば1から600フレーム目までとすると、

avg = nanmean(FFF(:,1:600),2)*ones(1,size(FFF,2));

DF = (FFF - avg)./avg.*100;

例えば、intrinsic optical imagingやtranscranial Ca²⁺ imagingなどで繰り返し同じ刺激を与えた時の蛍光変化の時系列データの加算平均を計算する場合、

F₀は、刺激前の数フレームで良いのだが、長時間露光などの影響で褪色してしまう場合も考慮に入れて、刺激を入れるその都度ΔF/Fを計算して加算平均を計算すると良い。

疑似カラー表示画像を出力する

function outputColorMap(str6,fstart,fend,N,dir,min,max)

% Matrix to ColorMap and make tiffs

for fff = fstart:fend %1:size(str6,2)

close all;

fprintf('%d\n',fff);

str = str6(:,fff);

str2 = reshape(str,N,N);

h = figure('visible','off');

set(h, 'PaperPositionMode', 'manual');

set(h, 'PaperUnits', 'inches');

%set(h,'papersize',[3 3],'PaperPosition', [0, 0, 3, 3]);

ti(1) = 0; ti(2) = 0; ti(3) = 0; ti(4) = 0;

pos(3) = 3; pos(4) = 3;

set(gca,'Position',[ti(1) ti(2) 1-ti(3)-ti(1) 1-ti(4)-ti(2)]);

set(gcf, 'PaperSize', [pos(3)+ti(1)+ti(3) pos(4)+ti(2)+ti(4)]);

set(gcf, 'PaperPosition',[0 0 pos(3)+ti(1)+ti(3) pos(4)+ti(2)+ti(4)]);

set(h,'InvertHardcopy','off')

b = imagesc(str2);

axis ij

axis square

axis off

caxis([min max])

file = strcat(dir,'/Image_',num2str(fff),'.tif');

saveas(h,file,'tif');

end

ImageJで作ったROIを読み込む

ImageJで作ったROIをmatlabに読み込む関数は以下のリンクから得ることができます。

http://www.dylan-muir.com/articles/read_ij_roi/

http://www.mathworks.com/matlabcentral/fileexchange/32479-import-imagej-rois

[sROI] = ReadImageJROI(RoiSetFile);

使い方は簡単で、RoiSetFileに'.roi'ファイルかまたは'.zip'ファイルのフルパスを入力すれば読み込めます。

.zipの場合はセルを自動で作ってその中に全roi情報が格納されるようです。

セルの中に、mnCoordinatesという項目があってそこに、x座標とy座標の情報が登録されています。よって、

xv = sROI{1,m}.mnCoordinates(:,1);

yv = sROI{1,m}.mnCoordinates(:,2);

とすれば、m番目のROIの座標情報が取得できます。

さて、問題は、取得できる情報は輪郭だけなので、囲んだ内部の情報が得たい場合は少し工夫が必要なようです。

今回は愚直に、512x512全ピクセルをスキャンして、その輪郭の内部にあるかないかを判定する関数、inpolygonを用いました。

あまり賢い方法とは言えないですが、これなら確実にROIの内部の座標が全部、indexに返ってきます。

indexの情報は新たに作ったROIindexというセルに全部放り込んで、後で取り出せるようにしておきます。

何かもっといい方法があったら教えて下さい。

例) 512x512画像の場合

function [ROIindex,numROI] = ReadROIindex(RoiSetFile)

% RoiSetFile = '0044-0257-0253.roi'

% RoiSetFile = 'RoiSet.zip'

[sROI] = ReadImageJROI(RoiSetFile);

numROI = size(sROI,2);

ROIindex = cell([1,numROI]);

for m = 1:numROI

xv = sROI{1,m}.mnCoordinates(:,1);

yv = sROI{1,m}.mnCoordinates(:,2);

index = [];

for x = 1:512

for y = 1:512

B = (x-1).*512+y;

in = inpolygon(x,y,xv,yv);

if in == 1

index = cat(1,index,B);

end

ROIindex{1,m} = index;

end

RoiSetFile = strcat(RoiSavedDir,'/RoiSet.zip');

[ROIindex,numROI] = ReadROIindex(RoiSetFile);

file = strcat(RoiSavedDir,'/ROIindex.mat');

save(file,'ROIindex');

これで二次元の時系列画像データを一次元の時系列データに縮約できました。

あとはこのデータを、一次元時系列データの時と同様に計算してゆけば良いというわけです。

以上。