视频来自：11_Utilities_哔哩哔哩_bilibili

pdf：​​​​​​The R-fMRI Course | The R-fMRI Network

目录

1. Utilities

1.1. Processing for Animal Data

1.2. Standardization

1.3. Utilities

1.4. Matlab作图

2. 脑影像编程

2.1. Matlab

3. r-fMRI和DPABI学习总结

3.1. 视频学习总结

3.2. 心得体会

1. Utilities

1.1. Processing for Animal Data

（1）在老师的版本中有dpabi→DPARSF 3.2 for Monkey Data

 

1.2. Standardization

（1）框架

 

1.3. Utilities

（1）DPABI→Utilities→Multiple T1 Images Averager

        ①扫了很多T1像，但是T1像单个信噪比不是很高，可以对齐到一块儿

        ②界面

 

（2）T1 Image Detacer

        ①功能：去掉T1像脸的部分

        ②界面

 

（3）ROI Signal Extractor

        ①如果得到了一个Cluster，并想提取它的ALFF值

        ②兴趣区（Region of Interest, ROI）

        ③界面

        ④类别（最多的即红框那种）

 

（4）Image Reslicer

        ①界面

        ②Reference：可以参考另外一个图像

（5）Image Calculator

        ①可以不写代码但实现计算

        ②界面

        ③Help：可提供示例计算

        ④在Expression里写了表达式，并在前缀Prefix里指定

        ⑤代码示例（y_ReadRPI是图像自动调整左右，y_ReadAll是所有文件都能读取）

 

1.4. Matlab作图

（1）被试的ALFF和MMSE散点图

        ①代码

[r p]=corrcoef(ALFF,MMSE)
plot(ALFF,MMSE,'o') %o是指画图方式，写其他的也可以

%xlabel('ALFF')
%ylabel('MMSE')

Isline

        ②结果示例

 

2. 脑影像编程

2.1. Matlab

（1）常用工作路径（Workspace）和命令行窗口（Command Window）

（2）可以在工作路径点开文件直接改参数，并使用如下代码直接按照DPABI RUN

DPARSFA_run(文件名)

（3）数组

        ①行向量

        ②列向量

 （4）取ALFF值的代码（可能也会根据自己文件不同有变化吧）

ALFF = ROISignals(:,1)

（5）创建空数组示例（此时Variables会多出一个类似文件的东西）

MMSE = []

（6）查询函数用法

doc 函数名

（7）TAB键自动补齐代码

（8）edit y_SCA.m（增加了自己的注释，但是感觉，全部看完真的有必要吗？虽然严老师带着全部过了一遍）

function [FCBrain, Header] = y_SCA(AllVolume, ROIDef, OutputName, MaskData, IsMultipleLabel, ROISelectedIndex, IsNeedDetrend, Band, TR, TemporalMask, ScrubbingMethod, ScrubbingTiming, Header, CUTNUMBER)
% [FCBrain, Header] = y_SCA(AllVolume, ROIDef, OutputName, MaskData, IsMultipleLabel, ROISelectedIndex, IsNeedDetrend, Band, TR, TemporalMask, ScrubbingMethod, ScrubbingTiming, Header, CUTNUMBER)
% Calculate Functional Connectivity by Seed based Correlation Anlyasis

% Input(这是每个输入参数可以接收的东西):
% 	AllVolume		-	4D data matrix (DimX*DimY*DimZ*DimTimePoints) or the directory of 3D image data file or the filename of one 4D data file
%   ROIDef              ROI definition, cells. Each cell could be:
%                       -1. 3D mask martrix (DimX*DimY*DimZ)
%                       -2. Series matrix (DimTimePoints*1)
%                       -3. REST Sphere Definition
%                       -4. .img/.nii/.nii.gz mask file
%                       -5. .txt Series. If multiple columns, when IsMultipleLabel==1: each column is a seperate seed series
%                                                             when IsMultipleLabel==0: average all the columns and take the mean series (one column) as seed series
%	OutputName  	-	Output filename
% 	MaskData		-   The Brain Mask matrix (DimX*DimY*DimZ) or the Brain Mask file name
%   IsMultipleLabel -   1: There are multiple labels in the ROI mask file. Will extract each of them. (e.g., for aal.nii, extract all the time series for 116 regions)
%                       0 (default): All the non-zero values will be used to define the only ROI
%   ROISelectedIndex -  Only extract ROIs defined by ROISelectedIndex. Empty means extract all non-zero ROIs.
%   IsNeedDetrend   -   0: Dot not detrend; 1: Use Matlab's detrend
%   Band            -   Temporal filter band: matlab's ideal filter e.g. [0.01 0.08]
%   TR              -   The TR of scanning. (Used for filtering.)
%   TemporalMask    -   Temporal mask for scrubbing (DimTimePoints*1)
%                   -   Empty (blank: '' or []) means do not need scrube. Then ScrubbingMethod and ScrubbingTiming can leave blank
%   ScrubbingMethod -   The methods for scrubbing.
%                       -1. 'cut': discarding the timepoints with TemporalMask == 0
%                       -2. 'nearest': interpolate the timepoints with TemporalMask == 0 by Nearest neighbor interpolation 
%                       -3. 'linear': interpolate the timepoints with TemporalMask == 0 by Linear interpolation
%                       -4. 'spline': interpolate the timepoints with TemporalMask == 0 by Cubic spline interpolation
%                       -5. 'pchip': interpolate the timepoints with TemporalMask == 0 by Piecewise cubic Hermite interpolation
%   ScrubbingTiming -   The timing for scrubbing.
%                       -1. 'BeforeFiltering': scrubbing (and interpolation, if) before detrend (if) and filtering (if).
%                       -2. 'AfterFiltering': scrubbing after filtering, right before extract ROI TC and FC analysis
%   Header          -   If AllVolume is given as a 4D Brain matrix, then Header should be designated.
%   CUTNUMBER       -   Cut the data into pieces if small RAM memory e.g. 4GB is available on PC. It can be set to 1 on server with big memory (e.g., 50GB).
%                       default: 10
% Output:
%	FCBrain         -   the FC of the final ROI (Only Use it correctly with only one seed)
%   Header          -   The NIfTI Header
%   All the FC images will be output as where OutputName specified.
%-----------------------------------------------------------


if ~exist('IsMultipleLabel','var')
    IsMultipleLabel = 0;
end

if ~exist('CUTNUMBER','var')
    CUTNUMBER = 10;
end

if ~exist('ROISelectedIndex','var')
    ROISelectedIndex = [];
end

%显示电脑时间，matlab里直接输入cputime就好了，也可以分别输入tic和toc，记录之间的时间差
theElapsedTime =cputime;
fprintf('\n\t Calculating Functional Connectivity by Seed based Correlation Anlyasis...');

%可以判断是不是文件，是文件就读进来
if ~isnumeric(AllVolume)
    [AllVolume,VoxelSize,theImgFileList, Header] =y_ReadAll(AllVolume);
end

%是矩阵的话读出矩阵的属性
[nDim1 nDim2 nDim3 nDimTimePoints]=size(AllVolume);
BrainSize = [nDim1 nDim2 nDim3];
VoxelSize = sqrt(sum(Header.mat(1:3,1:3).^2));

%看是不是Mask矩阵，但是多数为文件名
if ischar(MaskData)
    fprintf('\nLoad mask "%s".\n', MaskData);
    %如果是空字符串就全脑都要计算，如果不为空就读，变成同一个朝向
    if ~isempty(MaskData)
        [MaskData,MaskVox,MaskHead]=y_ReadRPI(MaskData);
        %看数值对不对应，不对应会报错
        if ~all(size(MaskData)==[nDim1 nDim2 nDim3])
            error('The size of Mask (%dx%dx%d) doesn''t match the required size (%dx%dx%d).\n',size(MaskData), [nDim1 nDim2 nDim3]);
        end
        %logical()是二值化
        MaskData = double(logical(MaskData));
    else
        MaskData=ones(nDim1,nDim2,nDim3);
    end
end

% Convert into 2D，转了的话matlab方便做矩阵运算
AllVolume=reshape(AllVolume,[],nDimTimePoints)';
%把时间转成第一维，效率更高
% AllVolume=permute(AllVolume,[4,1,2,3]); % Change the Time Course to the first dimention
% AllVolume=reshape(AllVolume,nDimTimePoints,[]);

%变成一维
MaskDataOneDim=reshape(MaskData,1,[]);
%找出里面不为0的
MaskIndex = find(MaskDataOneDim);
%对于找出的体素，只保留脑内的
AllVolume=AllVolume(:,MaskIndex);

% Scrubbing
if exist('TemporalMask','var') && ~isempty(TemporalMask) && ~strcmpi(ScrubbingTiming,'AfterFiltering')
    if ~all(TemporalMask)
        AllVolume = AllVolume(find(TemporalMask),:); %'cut'
        if ~strcmpi(ScrubbingMethod,'cut')
            xi=1:length(TemporalMask);
            x=xi(find(TemporalMask));
            AllVolume = interp1(x,AllVolume,xi,ScrubbingMethod);
        end
        nDimTimePoints = size(AllVolume,1);
    end
end


% Detrend
if exist('IsNeedDetrend','var') && IsNeedDetrend==1
    %AllVolume=detrend(AllVolume);
    fprintf('\n\t Detrending...');
    SegmentLength = ceil(size(AllVolume,2) / CUTNUMBER);
    for iCut=1:CUTNUMBER
        if iCut~=CUTNUMBER
            Segment = (iCut-1)*SegmentLength+1 : iCut*SegmentLength;
        else
            Segment = (iCut-1)*SegmentLength+1 : size(AllVolume,2);
        end
        AllVolume(:,Segment) = detrend(AllVolume(:,Segment));
        fprintf('.');
    end
end

% Filtering：前面的不存在直接跳过了，存在且不是空就滤波
if exist('Band','var') && ~isempty(Band)
    fprintf('\n\t Filtering...');
    %内存有限，一次只操作一部分
    SegmentLength = ceil(size(AllVolume,2) / CUTNUMBER);
    for iCut=1:CUTNUMBER
        if iCut~=CUTNUMBER
            Segment = (iCut-1)*SegmentLength+1 : iCut*SegmentLength;
        else
            Segment = (iCut-1)*SegmentLength+1 : size(AllVolume,2);
        end
        AllVolume(:,Segment) = y_IdealFilter(AllVolume(:,Segment), TR, Band);
        fprintf('.');
    end
end



% Scrubbing after filtering
if exist('TemporalMask','var') && ~isempty(TemporalMask) && strcmpi(ScrubbingTiming,'AfterFiltering')
    if ~all(TemporalMask)
        AllVolume = AllVolume(find(TemporalMask),:); %'cut'
        if ~strcmpi(ScrubbingMethod,'cut')
            xi=1:length(TemporalMask);
            x=xi(find(TemporalMask));
            AllVolume = interp1(x,AllVolume,xi,ScrubbingMethod);
        end
        nDimTimePoints = size(AllVolume,1);
    end
end


% Extract the Seed Time Courses
%做功能链接需要提取某一行的时间序列

SeedSeries = [];
MaskROIName=[];

for iROI=1:length(ROIDef)
    IsDefinedROITimeCourse =0;
    if strcmpi(int2str(size(ROIDef{iROI})),int2str([nDim1, nDim2, nDim3]))  %ROI Data
        MaskROI = ROIDef{iROI};
        MaskROIName{iROI} = sprintf('Mask Matrix definition %d',iROI);
    elseif size(ROIDef{iROI},1) == nDimTimePoints %Seed series %strcmpi(int2str(size(ROIDef{iROI})),int2str([nDimTimePoints, 1])) %Seed series
        SeedSeries{1,iROI} = ROIDef{iROI};
        IsDefinedROITimeCourse =1;
        MaskROIName{iROI} = sprintf('Seed Series definition %d',iROI);
    elseif strcmpi(int2str(size(ROIDef{iROI})),int2str([1, 4]))  %Sphere ROI definition: CenterX, CenterY, CenterZ, Radius
        %生成一个球
        MaskROI = y_Sphere(ROIDef{iROI}(1:3), ROIDef{iROI}(4), Header);
        MaskROIName{iROI} = sprintf('Sphere definition (CenterX, CenterY, CenterZ, Radius): %g %g %g %g.',ROIDef{iROI});
    elseif exist(ROIDef{iROI},'file')==2	% Make sure the Definition file exist
        [pathstr, name, ext] = fileparts(ROIDef{iROI});
        %.txt后缀表明可能给的时间序列
        if strcmpi(ext, '.txt'),
            TextSeries = load(ROIDef{iROI});
            if IsMultipleLabel == 1
                for iElement=1:size(TextSeries,2)
                    MaskROILabel{1,iROI}{iElement,1} = ['Column ',num2str(iElement)];
                end
                SeedSeries{1,iROI} = TextSeries;
            else
                SeedSeries{1,iROI} = mean(TextSeries,2);
            end
            IsDefinedROITimeCourse =1;
            MaskROIName{iROI} = ROIDef{iROI};
        elseif strcmpi(ext, '.img') || strcmpi(ext, '.nii') || strcmpi(ext, '.gz')
            %The ROI definition is a mask file
            
            MaskROI=y_ReadRPI(ROIDef{iROI});
            %检验masksize和输出数据一不一样
            if ~strcmpi(int2str(size(MaskROI)),int2str([nDim1, nDim2, nDim3]))
                error(sprintf('\n\tMask does not match.\n\tMask size is %dx%dx%d, not same with required size %dx%dx%d',size(MaskROI), [nDim1, nDim2, nDim3]));
            end

            MaskROIName{iROI} = ROIDef{iROI};
        else
            error(sprintf('Wrong ROI file type, please check: \n%s', ROIDef{iROI}));
        end
        
    else
        error(sprintf('File doesn''t exist or wrong ROI definition, please check: %s.\n', ROIDef{iROI}));
    end

    %如果还没有提取时间序列，则提取
    if ~IsDefinedROITimeCourse
        % Speed up! YAN Chao-Gan 101010.
        MaskROI=reshape(MaskROI,1,[]);
        MaskROI=MaskROI(MaskIndex); %Apply the brain mask
        
        if IsMultipleLabel == 1

            if ~isempty(ROISelectedIndex) && ~isempty(ROISelectedIndex{iROI})
                Element=ROISelectedIndex{iROI};
            else
                Element = unique(MaskROI);
                %把NaN和0扔掉，它们不会是ROI
                Element(find(isnan(Element))) = []; % ignore background if encoded as nan. Suggested by Dr. Martin Dyrba
                Element(find(Element==0)) = []; % This is the background 0
            end

            SeedSeries_MultipleLabel = zeros(nDimTimePoints,length(Element));
            for iElement=1:length(Element)
                
                %求均值
                SeedSeries_MultipleLabel(:,iElement) = mean(AllVolume(:,find(MaskROI==Element(iElement))),2);
                
                MaskROILabel{1,iROI}{iElement,1} = num2str(Element(iElement));

            end
            SeedSeries{1,iROI} = SeedSeries_MultipleLabel;
        else
            SeedSeries{1,iROI} = mean(AllVolume(:,find(MaskROI)),2);
        end
    end
end


%Merge the seed series cell into seed series matrix
SeedSeries = double(cell2mat(SeedSeries)); %Suggested by H. Baetschmann.  % SeedSeries = cell2mat(SeedSeries);


%Save the ROI averaged time course to disk for further study
[pathstr, name, ext] = fileparts(OutputName);

save([fullfile(pathstr,['ROI_', name]), '.mat'], 'SeedSeries')
save([fullfile(pathstr,['ROI_', name]), '.txt'], 'SeedSeries', '-ASCII', '-DOUBLE','-TABS')

%Write the order key file as .tsv
fid = fopen([fullfile(pathstr,['ROI_OrderKey_', name]), '.tsv'],'w');
if IsMultipleLabel == 1
    if size(MaskROILabel,2) < length(ROIDef) %YAN Chao-Gan, 131124. To avoid if the labels of the last ROI has been defined.
        MaskROILabel{1,length(ROIDef)} = []; % Force the undefined cells to empty
    end
    fprintf(fid,'Order\tLabel in Mask\tROI Definition\n');
    iOrder = 1;
    for iROI=1:length(ROIDef)
        if isempty(MaskROILabel{1,iROI})
            fprintf(fid,'%d\t\t%s\n',iOrder,MaskROIName{iROI});
            iOrder = iOrder + 1;
        else
            for iElement=1:length(MaskROILabel{1,iROI})
                fprintf(fid,'%d\t%s\t%s\n',iOrder,MaskROILabel{1,iROI}{iElement,1},MaskROIName{iROI});
                iOrder = iOrder + 1;
            end
        end
    end
else
    fprintf(fid,'Order\tROI Definition\n');
    for iROI=1:length(ROIDef)
        fprintf(fid,'%d\t%s\n',iROI,MaskROIName{iROI});
    end
end
fclose(fid);



% FC calculation
AllVolume = AllVolume-repmat(mean(AllVolume),size(AllVolume,1),1);
AllVolumeSTD= squeeze(std(AllVolume, 0, 1));
AllVolumeSTD(find(AllVolumeSTD==0))=inf;

SeedSeries=SeedSeries-repmat(mean(SeedSeries),size(SeedSeries,1),1);
SeedSeriesSTD=squeeze(std(SeedSeries,0,1));

Header.pinfo = [1;0;0];
Header.dt    =[16,0];

for iROI=1:size(SeedSeries,2)
    
    FC=SeedSeries(:,iROI)'*AllVolume/(nDimTimePoints-1);
    FC=(FC./AllVolumeSTD)/SeedSeriesSTD(iROI);
    
    %把二维矩阵变成三维
    FCBrain=zeros(size(MaskDataOneDim));
    FCBrain(1,MaskIndex)=FC;
    FCBrain=reshape(FCBrain,nDim1, nDim2, nDim3);
    
    %Also produce the results after Fisher's r to z transformation
    zFCBrain = (0.5 * log((1 + FCBrain)./(1 - FCBrain))) .* (MaskData~=0);

    [pathstr, name, ext] = fileparts(OutputName);
    if size(SeedSeries, 2)>1
        %Save every maps from result maps
        y_Write(FCBrain,Header,[fullfile(pathstr,['ROI',num2str(iROI),name, ext])]);
        y_Write(zFCBrain,Header,[fullfile(pathstr,['zROI',num2str(iROI),name, ext])]);
    elseif size(SeedSeries, 2)==1,
        %Save one map
        y_Write(FCBrain,Header,[fullfile(pathstr,[name, ext])]);
        y_Write(zFCBrain,Header,[fullfile(pathstr,['z',name, ext])]);
    end
end


theElapsedTime = cputime - theElapsedTime;
fprintf('\n\t Calculating Functional Connectivity by Seed based Correlation Anlyasis finished, elapsed time: %g seconds.\n', theElapsedTime);

3. r-fMRI和DPABI学习总结

3.1. 视频学习总结

（1）感觉DPABI略小众，很多时候报错在网上也搜不出解决办法

（2）视频没有字幕导致初学者完全不知道专业术语或英文缩写是在说什么

（3）可能讲课范围涉及得不是很全面，不知道处理出来的脑影像有什么意义

（4）没有什么基础知识学习，直接就开论文讲解了

（5）最后一节的代码也太...断崖式了，有一种“明明是学计算机的学生却听不懂医学老师讲的代码”的感觉

3.2. 心得体会

（1）功能按钮确实是学到了不少

（2）为什么现在的软件普遍都要做得那么复杂？

（3）静息脑影像会不会一定程度上还是有点鸡肋？（纯纯猜测）但是确实动态下头动可能太严重了

（4）感觉也没有教会我处理数据？好像只说了导出来然后画个ALFF和S什么什么的回归图然后我也不知道这俩是啥

（5）路漫漫其修远兮...