精华内容
下载资源
问答
  • nilearn资源包

    2018-04-10 11:07:11
    nilearn开发医疗影像处理,大数据领域的爆发和硬件性能的飞速提升给医疗人工智能企业的发展提供了“燃料”
  • nibabel_nilearn_tutorial python教程中的神经成像入门材料,由多伦多Compute Ontario Summer School正式提供,2017年7月 连接到SciNet jupyter集线器。 首先,您需要ssh到niagara节点之一以克隆此存储库 ssh ...
  • 最近在学习Nilearn包的使用,现将学习笔记和思考发布在这里,供大家参考。 以下训练的数据均来自于官网的练习数据 haxby2001的数据 一、导入核磁数据 这个没啥特别的,指定图像地址就好了,唯一要注意的是,需要是4D...

    最近在学习Nilearn包的使用,现将学习笔记和思考发布在这里,供大家参考。
    以下训练的数据均来自于官网的练习数据 haxby2001的数据

    一、导入核磁数据

    这个没啥特别的,指定图像地址就好了,唯一要注意的是,需要是4D的nii.gz文件格式,这个可以通过其他软件转换。

    #导入数据就是命名nii.gz的数据地址
    file_name = r'D:\jupyter\nilearn\haxby2001\subj2\bold.nii.gz'
    

    可选教程:如何查看图像

    #查看数据的方法,用 plotting.view_img函数
    from nilearn import plotting
    plotting.view_img(file_name)
    

    二、定义一个mask

    这个定义的mask,就是ROI。
    需要注意的是:是不是一定需要mask?
    答案是:不是必须,但是建议。

    定义mask的好处在于,能容易出结果,范围小了,也更适合讲故事了。
    如果不要mask,那就是全脑的体素进行分析,那矩阵量就太大了,但是也可以分析,这是另一种方法,后面会介绍。这里先介绍简单的。

    mask_filename = r'D:\jupyter\nilearn\haxby2001\subj2\mask4_vt.nii.gz'
    

    可选教程:如何查看ROI

    #如何查看ROI
    anatomical_image = r'D:\jupyter\nilearn\haxby2001\subj2\anat.nii.gz'
    # cmap:The colormap to use
    plotting.plot_roi(mask_filename, bg_img=anatomical_image, cmap='Paired')
    

    在这里插入图片描述

    mask_filename的图像:
    plotting.view_img(mask_filename)
    在这里插入图片描述
    anatomical_image的图像:
    plotting.view_img(anatomical_image)
    在这里插入图片描述
    补充:
    1、这个anatomical image图像,其实就是T1像,就是患者的结构像
    2、plot_roi函数的用法:plot_roi(roi_img, bg_img=, cut_coords=None, output_file=None, display_mode=‘ortho’, figure=None, axes=None, title=None, annotate=True, draw_cross=True, black_bg=‘auto’, threshold=0.5, alpha=0.7, cmap=<matplotlib.colors.LinearSegmentedColormap object at 0x000002144F1BD490>, dim=‘auto’, vmin=None, vmax=None, resampling_interpolation=‘nearest’, view_type=‘continuous’, linewidths=2.5, **kwargs)

    3、加载标签 Label

    把标签存储在 CSV 文件中,以空格分隔,然后使用 Pandas 将它们加载到数组中。

    import pandas as pd
    # Load behavioral information
    target = r'D:\jupyter\nilearn\haxby2001\subj2\labels.txt'
    behavioral = pd.read_csv(target, delimiter=' ')
    print(behavioral)
    

    在这里插入图片描述
    labels.txt文件:
    在这里插入图片描述

    conditions = behavioral['labels']
    

    输出的conditions是panda的序列:pandas.core.series.Series,一共有1452个数字,刚好对应bold.nii.gz的1452个时间点。
    这里我需要跟新手提一嘴,一般常规的功能像分析,是要去除前10个时间点的,所以使用自己的数据的时候,就需要和去除后的时间点对应上。

    4、将分析限制为猫和人脸

    上面的labels里面包含了很多conditions。因此,数据相当大。并非所有这些数据都对我们感兴趣进行解码,因此我们将仅保留与人脸或猫对应的fmri信号。我们创建属于条件的样本的掩码;然后将此掩码应用于fmri数据,以将分类限制为面部与猫的区分。

    condition_mask = conditions.isin(['face', 'cat'])
    

    pandas库Series的函数的isin函数(A的元素是否在B当中),输出True和False.
    由于数据在一张大的4D图像中,我们需要使用index_img来轻松进行拆分。

    from nilearn.image import index_img
    fmri_niimgs = index_img(fmri_filename, condition_mask)
    

    label_mask输出的是True和False,需要得到具体的label

    conditions = conditions[condition_mask] #conditions剩下216行
    # Convert to numpy array
    conditions = conditions.values #把conditions的值提取出来,就是216个face或cat的字符串,numpy.ndarray格式
    

    5、定义评估器

    前面我们得到了目标label,得到了label对应的激活时间点,也得到了局部的ROI体素。首先,有一个很基础知识需要知道,就是在bold信号里,大脑某个体素值增高,就代表该区域被激活,反之就是抑制。
    因此,事实上,以上这么多步骤,其目的就是:把某区域的体素的灰度值增高与降低同标签label进行建模预测。如果能预测成功,那这块区域激活或者抑制就与这类标签具有很强的关系,也就达到了我们的目的。
    上面是废话,下面开始继续:
    1.定义评估器:

    from nilearn.decoding import Decoder
    decoder = Decoder(estimator='svc', mask=mask_filename, standardize=True)
    

    standardize=True , 对每个体素上的时间序列数据做0-1归一化。
    2.拟合数据

    decoder.fit(fmri_niimgs, conditions)
    

    3.检测预测

    prediction = decoder.predict(fmri_niimgs)
    print(prediction)
    

    在这里插入图片描述
    当然,这种流程的预测是没有意义的,原因是没有划分测试集和验证集,相当于提高告诉答案,然后再考试,准确率当然是1了。

    6、使用交叉验证测量预测分数

    让我们忽略训练期间的最后 30 个数据点,并测试对这最后 30 个点的预测:

    fmri_niimgs_train = index_img(fmri_niimgs, slice(0, -30))
    fmri_niimgs_test = index_img(fmri_niimgs, slice(-30, None))
    conditions_train = conditions[:-30]
    conditions_test = conditions[-30:]
    decoder = Decoder(estimator='svc', mask=mask_filename, standardize=True)
    decoder.fit(fmri_niimgs_train, conditions_train)
    prediction = decoder.predict(fmri_niimgs_test)
    # The prediction accuracy is calculated on the test data: this is the accuracy
    # of our model on examples it hasn't seen to examine how well the model perform
    # in general.
    print("Prediction Accuracy: {:.3f}".format(
        (prediction == conditions_test).sum() / float(len(conditions_test))))
    

    在这里插入图片描述
    0.767的准确率看起来还不错,接下来需要进行手动的K折交叉验证,然后输出AUC曲线下面积来正确评估这个ROI mask来预测行为的准确率。

    from sklearn.model_selection import KFold
    cv = KFold(n_splits=5)
    # The "cv" object's split method can now accept data and create a
    # generator which can yield the splits.
    fold = 0
    for train, test in cv.split(conditions):
        fold += 1
        decoder = Decoder(estimator='svc', mask=mask_filename, standardize=True)
        decoder.fit(index_img(fmri_niimgs, train), conditions[train])
        prediction = decoder.predict(index_img(fmri_niimgs, test))
        print(
            "CV Fold {:01d} | Prediction Accuracy: {:.3f}".format(
                fold,
                (prediction == conditions[test]).sum() / float(len(
                    conditions[test]))))
    

    在这里插入图片描述
    AUC曲线下面积的代码日后补充进来。

    6.2 与解码器的交叉验证

    这一节还是说交叉验证的事儿,上一节是手动进行K折交叉验证,这里是用评估器默认的功能来实现交叉验证。
    解码器还默认实现交叉验证循环并返回一个形状数组(交叉验证参数,n_folds),其中将准确度定义为 评分参数来使用准确度分数来衡量其性能。

    n_folds = 5
    decoder = Decoder(
        estimator='svc', mask=mask_filename,
        standardize=True, cv=n_folds,
        scoring='accuracy'
    )
    decoder.fit(fmri_niimgs, conditions)
    

    以上是根据时间点来计算的,事实上,静息态功能磁共振在设计上,经常会有block设计,block就是组块的意思,通俗来讲,就是静息-任务-静息-任务-静息。。。这样的实验。
    那么来说,对于这种组块设计的实验,其实应该加上组块的信息进行多个事件组的分析,也就是group(在这里应该是block)分析,其结果会更好点。
    详细解释:
    试想一下,把bold.nii.gz里面单个时间点当成1个人,那么就有216个人,我们去预测哪些人是cat,哪些人是face。
    在传统临床上,一般是把216个人使用随机抽样的方法,变成训练集和验证集(有人会提到应该分3个集合,但是216的数据量还达不到那个程度),然后在训练集上使用K折交叉验证,最后用验证集去验证。
    但是在核磁里,时间点是连续的,我们没办法确定具体这个时间点,这个患者的大脑状态就对应着cat或者face,时间分辨率就不是核磁的长处,所以这里就产生了误差。
    所以更好的方法是,因为时间点是连续的,静息-任务-静息的block是设计好的,那么研究组间差异,会比研究单个时间点更有价值。即这里的216个人是按号码排好的,face和cat是固定出现的,那么我们去用单个block去机器学习,然后统计单个block的预测准确率,这种方法会更适用于核磁。所以这里的交叉验证方法也就相应的使用留一法,即留下一个block。但是要记住,对单个时间点进行留一法是非常错误的。
    官网数据的block重复出现了12次,每次出现的时候face和cat包含了9个时间点。如果是自己的数据的话,block重复的次数肯定要2和更高,不然就会报错。
    block留一法:

    session_label = behavioral['chunks'][condition_mask]
    from sklearn.model_selection import LeaveOneGroupOut
    cv = LeaveOneGroupOut()
    decoder = Decoder(estimator='svc', mask=mask_filename, standardize=True, cv=cv)
    decoder.fit(fmri_niimgs, conditions, groups=session_label)
    print(decoder.cv_scores_)
    

    在这里插入图片描述
    可以看出,在第八次的block中,预测有点问题,其他block都非常棒。
    下面是12个block中,cat和face的分布:
    在这里插入图片描述
    暂时没看出第八次block有啥特别的地方,不负责盲猜可能和头动有一定关系。

    8、检查模型权重,并转换成nii图像

    coef_ = decoder.coef_
    print(coef_)
    

    这是一个 numpy 数组,每个体素只有一个系数,所以就是[1, 464]的shape。

    coef_img = decoder.coef_img_['face']
    decoder.coef_img_['face'].to_filename('haxby_svc_weights.nii.gz')
    

    保存带有权重的nii图像

    plotting.view_img(
        decoder.coef_img_['face'], bg_img=anatomical_image,
        title="SVM weights", dim=-1
    )
    

    在这里插入图片描述
    使用虚拟分类器(dummy_classification)来作为底线(假设基线)检验上面的SVC分类器性能。

    dummy_decoder = Decoder(estimator='dummy_classifier', mask=mask_filename, cv=cv)
    dummy_decoder.fit(fmri_niimgs, conditions, groups=session_label)
    # Now, we can compare these scores by simply taking a mean over folds
    print(dummy_decoder.cv_scores_)
    

    在这里插入图片描述
    可以看出,虚拟分类器的预测率很低,这就是随便使用一个分类器能达到的底线预测率,而SVC确实提高了很多。
    关于虚拟分类器(dummy_classification)的解释:
    虚拟分类器为您提供了一种“基准”性能测量 - 即,即使只是猜测,也应该达到预期的成功率。
    假设您想确定某个给定的对象是拥有还是不具有某种属性。如果您已经分析了大量这些对象并且发现90%包含目标属性,那么猜测该对象的每个未来实例都拥有目标属性,这使您有90%的正确猜测的可能性。以这种方式构建您的猜测等同于您在引用的文档中使用most_frequent方法。由于许多机器学习任务试图增加(例如)分类任务的成功率,所以评估基线成功率可以为分类器应该超出的最小值提供底限值。在上面讨论的假设中,您会希望分类器的准确率达到90%以上,因为即使是“虚拟”分类器,成功率也可达到90%。
    如果使用上面讨论的数据训练带有stratified参数的虚拟分类器,则该分类器将预测其遇到的每个对象都有90%的概率拥有目标属性。这与训练具有most_frequent参数的虚拟分类器不同,因为后者会猜测未来对象具有目标属性。

    理解以上,简单的Nilearn就算入门了,后续会继续更新我的学习笔记。我认为,这个方法还是很不错的,简单不复杂,而且契合现在机器学习的热潮,具有应用潜力,比如在我的领域,我准备用这个方法找分类针刺刺激的真假穴的脑区。

    我学习完这个初步教程,我还有两个问题:
    1、对于机器学习方法来说,使用原始未经处理的数据还是使用转换后的数据,对于结果预测率来说没有区别,重要的是,输入的变量有没有包含决定性变量。那么,在核磁领域,是经过slice timing、head motion、realign、segmentation等等步骤后的数据用来预测好点,还是原始数据好点?
    个人见解:其实没太大区别,未处理的数据和处理后的数据,只要他们都是nii.gz格式的图像,那就能用上述方法分析,然后看那种分析准确率更高。
    2、需不需要做slice timing,把时间点基线统一?
    个人见解:这个等我通过学习了更高阶的nilearn方法,阅读了他们团队的文章后,再回答这个问题。
    3、对于非block设计,如病人干预前后,那么1个病人有2个核磁状态,假如有100个病人,对于这种实验,该怎么使用nilearn去分类不同干预措施呢?
    个人见解:等我后续学完更高阶的nilearn,学习完后,我再来更新这个问题。
    提供一个思路:对于严格遵循随机对照方法的随机对照试验,因为其基线是差异很小的,如果这个时候,基线的差异小到可以忽略,那么就可以把干预后的结果作为label,然后抽样就在100个病人里随机抽,这样就回到了传统的机器学习方法。值得一试。

    Over.
    clancy wu

    展开全文
  • from nilearn import datasets dataset = datasets.fetch_atlas_harvard_oxford('cort-maxprob-thr25-2mm') atlas_filename = dataset.maps labels = dataset.labels print('Atlas ROIs are located in nifti image ...
    from nilearn import datasets
    dataset = datasets.fetch_atlas_harvard_oxford('cort-maxprob-thr25-2mm')
    atlas_filename = dataset.maps
    labels = dataset.labels
    print('Atlas ROIs are located in nifti image (4D) at: %s' %
          atlas_filename)  # 4D data
    # One subject of brain development fmri data
    data = datasets.fetch_development_fmri(n_subjects=1)
    fmri_filenames = data.func[0]
    

    在这里插入图片描述
    在这里插入图片描述

    解释一下:首先是获取了哈佛的地图所有信息(dataset),atlas_filename就是获取哈佛地图,其次获取label,label就是地图的标签。len(lables)=49,哈佛地图一共有49个标签。但是,需要注意的是,第0个标签是background,背景呈现用的,所以第0个标签是不作为ROI去提取时间序列的,所以masker默认是48。

    from nilearn.input_data import NiftiLabelsMasker
    masker = NiftiLabelsMasker(labels_img=atlas_filename, standardize=True,
                               memory='nilearn_cache', verbose=5)
    # Here we go from nifti files to the signal time series in a numpy
    # array. Note how we give confounds to be regressed out during signal
    # extraction
    time_series = masker.fit_transform(fmri_filenames, confounds=data.confounds)
    

    NiftiLabelsMasker 是个重要的函数,该函数可以根据地图标签,制作mask。可惜的是,这个mask无法被预览,是个特殊的格式,NiftiLabelsMasker数据格式。
    fmri_filenames是一个4D的时间序列。
    在这里插入图片描述
    time_series导师很普通的,是一个(168,48)的矩阵。一共有有48个标签,每个标签有168个时间点组成的时间序列,所以就是(168,48)。
    在这里插入图片描述

    计算并显示相关矩阵

    from nilearn.connectome import ConnectivityMeasure
    correlation_measure = ConnectivityMeasure(kind='correlation')
    correlation_matrix = correlation_measure.fit_transform([time_series])[0]
    
    # Plot the correlation matrix
    import numpy as np
    from nilearn import plotting
    # Make a large figure
    # Mask the main diagonal for visualization:
    np.fill_diagonal(correlation_matrix, 0)
    # The labels we have start with the background (0), hence we skip the
    # first label
    # matrices are ordered for block-like representation
    plotting.plot_matrix(correlation_matrix, figure=(10, 8), labels=labels[1:],
                         vmax=0.8, vmin=-0.8, reorder=True)
    

    在这里插入图片描述
    在这里插入图片描述

    这段的意思是:计算的方法是相关(即dparsf的full correlation),然后生成correlation_matrix,相关性矩阵,矩阵维度是(48,48),最后呈现出的就是矩阵图。

    同样的事情没有混淆,强调混淆的重要性

    time_series = masker.fit_transform(fmri_filenames)
    # Note how we did not specify confounds above. This is bad!
    
    correlation_matrix = correlation_measure.fit_transform([time_series])[0]
    
    # Mask the main diagonal for visualization:
    np.fill_diagonal(correlation_matrix, 0)
    
    plotting.plot_matrix(correlation_matrix, figure=(10, 8), labels=labels[1:],
                         vmax=0.8, vmin=-0.8, title='No confounds', reorder=True)
    
    plotting.show()
    

    这个就不展示了,图像基本都是强相关,没啥意义。
    这一步的主要意思是,计算时间序列相关性时,需要把白质、脑脊液等信号给回归掉,这个在dparsf里就是nuisance这一步,是自动处理的,就是不用太关注这个,了解就行。

    总结:总结一下提取时间序列的过程:
    1、导入4D的图像(因为3D的图像没有时间维度,所以必须是4D的volumes)
    2、选择带有标签的地图集,map and labels
    3、制作labels的masks,用NiftiLabelsMasker函数
    4、使用masker.fit_transform功能在4D的volumes上提取出masks的时间序列
    5、计算并呈现相关性矩阵
    重要补充:对于举一反三来说,假如我想提取masks中特定的几个labels,其实很简单,就是把labels和time_series根据索引值去掉那几项,因为这两个都是numpy格式,所以去掉很简单,这样就能得到期望中的几个labels的时间序列。很完美。

    补充代码:

    # 导入labels和altas
    from nilearn import datasets
    dataset = datasets.fetch_atlas_harvard_oxford('cort-maxprob-thr25-2mm')
    atlas_filename = dataset.maps
    labels = dataset.labels
    # 制作ROI的masker
    from nilearn.input_data import NiftiLabelsMasker
    masker = NiftiLabelsMasker(labels_img=atlas_filename, standardize=True)
    # 根据masker提取时间序列,混淆矩阵就是csv_file
    # 当然,对于已经去除掉混淆矩阵的nii_files,就不用啦
    time_series = masker.fit_transform(frmi_files, confounds=csv_file)
    

    如果是概率图谱提取时间序列,只需要把NiftiLabelsMasker换成NiftiMapsMasker就行了,很简单记住,确定的就是labels,因为规定好了,概率的就是maps,因为只是地图,没有规定好。

    from nilearn.input_data import NiftiMapsMasker
    masker = NiftiMapsMasker(maps_img=atlas_filename, standardize=True)
    time_series = masker.fit_transform(frmi_files, confounds=csv_file)
    

    功能连接

    # 导入地图
    from nilearn import datasets
    yeo = datasets.fetch_atlas_yeo_2011() #yeo 是3D图像,有17个labels
    # 导入数据
    data = datasets.fetch_development_fmri(n_subjects=10)
    # 提取yeo地图块的坐标
    from nilearn.input_data import NiftiLabelsMasker
    from nilearn.connectome import ConnectivityMeasure
    # ConenctivityMeasure from Nilearn uses simple 'correlation' to compute
    # connectivity matrices for all subjects in a list
    connectome_measure = ConnectivityMeasure(kind='correlation')
    # useful for plotting connectivity interactions on glass brain
    from nilearn import plotting
    # create masker to extract functional data within atlas parcels
    masker = NiftiLabelsMasker(labels_img=yeo['thick_17'], standardize=True,
                               memory='nilearn_cache')
    # extract time series from all subjects and concatenate them
    time_series = []
    for func, confounds in zip(data.func, data.confounds):
        time_series.append(masker.fit_transform(func, confounds=confounds))
    # calculate correlation matrices across subjects and display
    correlation_matrices = connectome_measure.fit_transform(time_series)
    # Mean correlation matrix across 10 subjects can be grabbed like this,
    # using connectome measure object
    mean_correlation_matrix = connectome_measure.mean_
    # grab center coordinates for atlas labels
    coordinates = plotting.find_parcellation_cut_coords(labels_img=yeo['thick_17'])
    # plot connectome with 80% edge strength in the connectivity
    plotting.plot_connectome(mean_correlation_matrix, coordinates,
                             edge_threshold="80%",
                             title='Yeo Atlas 17 thick (func)')
    

    在这里插入图片描述
    在这里插入图片描述

    在这里插入图片描述
    导入的yeo是一个bunch数据,有很多信息。导入的data同理,也是bunch数据。
    在这里插入图片描述
    在这里插入图片描述

    time_series是一个(10,168,17)的三维矩阵。10代表的是10个被试,168代表的是168个时间点,17代表的是yeo[‘thick_17’]里面有17个labels,这是一个皮层的数据。我记得yeo模板有好多,有皮层模板,有7大脑功能网络模板,这里的是皮层模板。如下展示,包含17个ROI。
    在这里插入图片描述
    记住:对于多个被试,就用for循环提取labels的时间序列,并制作成list。
    在这里插入图片描述
    有意思的来了,用correlation_measure计算相关矩阵时,对于list,就是list的每一条进行相关性分析,所以就是(10,17,17),很有意思。
    在这里插入图片描述
    mean_correlation_matrix就是计算10个被试的平均时间序列,记得要在mean后面加一个_,特殊记住就行了,按照这个格式来。
    在这里插入图片描述
    前面都没有啥,这一个非常重要,plotting.find_parcellation_cut_coords函数,只需要输入labels的图像,就可以自动提取labels图像里面的位置坐标,XYZ格式,非常实用!!!对于所有的有labels的图像,都可以使用这个方法。
    在这里插入图片描述
    这个就没啥说的,看代码就行。

    import numpy as np
    # Define a custom function to compute lag correlation on the time series
    def lag_correlation(time_series, lag):
        n_subjects = len(time_series)
        n_samples, n_features = time_series[0].shape
        lag_cor = np.zeros((n_subjects, n_features, n_features))
        for subject, serie in enumerate(time_series):
            for i in range(n_features):
                for j in range(n_features):
                    if lag == 0:
                        lag_cor[subject, i, j] = np.corrcoef(serie[:, i],
                                                             serie[:, j])[0, 1]
                    else:
                        lag_cor[subject, i, j] = np.corrcoef(serie[lag:, i],
                                                             serie[:-lag, j])[0, 1]
        return np.mean(lag_cor, axis=0)
    # Compute lag-0 and lag-1 correlations and plot associated connectomes
    for lag in [0, 1]:
        lag_correlation_matrix = lag_correlation(time_series, lag)
        plotting.plot_connectome(lag_correlation_matrix, coordinates,
                                 edge_threshold="90%",
                                 title='Lag-{} correlation'.format(
                                     lag))
    

    在这里插入图片描述
    如果是对称的矩阵,就会画一个无向的图,如果是用lag-1作为连接度量,就会画一个有向的图。
    尝试学习lag-1这个指标。可以先做lag0,就是普通的(同上面的功能连接),再作lag1,就是滞后1个时间点的。两个图作为补充提交。

    对于概率图谱的坐标提取,比较难,仅作了解就行,代码如下:

    from nilearn.input_data import NiftiMapsMasker
    # create masker to extract functional data within atlas parcels
    masker = NiftiMapsMasker(maps_img=difumo.maps, standardize=True,
                             memory='nilearn_cache')
    # extract time series from all subjects and concatenate them
    time_series = []
    for func, confounds in zip(data.func, data.confounds):
        time_series.append(masker.fit_transform(func, confounds=confounds))
    # calculate correlation matrices across subjects and display
    correlation_matrices = connectome_measure.fit_transform(time_series)
    # Mean correlation matrix across 10 subjects can be grabbed like this,
    # using connectome measure object
    mean_correlation_matrix = connectome_measure.mean_
    # grab center coordinates for probabilistic atlas
    coordinates = plotting.find_probabilistic_atlas_cut_coords(maps_img=difumo.maps)
    # plot connectome with 85% edge strength in the connectivity
    plotting.plot_connectome(mean_correlation_matrix, coordinates,
                             edge_threshold="85%",
                             title='DiFuMo with {0} dimensions (probabilistic)'.format(dim))
    plotting.show()
    
    展开全文
  • nilearn, 在 python 中,面向NeuroImaging的机器学习 nilearnNilearn是一个 python MODULE,用于快速简便的NeuroImaging数据统计学习。它利用了 scikit学习多变量统计,并应用了预测建模。分类。解码或者连接性分析...
  • GT-MVPA-nilearn 在法国马赛的协调的MVPA工作组(GT MVPA)中集中讨论和专业知识共享的存储库。 它主要集中在 python模块的使用上。 请使用“问题”部分作为问题/答案论坛! 另外,对于代码的初始上载,您可以...
  • 这是使用Nilearn解码的教程,它以Haxby 2001研究中猫辨别任务的数据为基础。fMRI解码入门教程功能性磁共振成像(FMRI,functional magnetic resona...

    97f8332581841a50a036c5407945a431.png

    这是使用Nilearn解码的教程,它以Haxby 2001研究中猫辨别任务的数据为基础。

    fMRI解码入门教程

    功能性磁共振成像(FMRI,functional magnetic resonance imaging)是一种新兴的神经影像学方式,其原理是利用磁振造影来测量神经元活动所引发之血液动力的改变。由于fMRI的非侵入性、没有辐射暴露问题与其较为广泛的应用,从1990年代开始就在脑部功能定位领域占有一席之地。目前主要是运用在研究人及动物的脑或脊髓。

    主要内容包括:

    1.从Haxby研究中检索并加载fMRI数据

    2.利用SVM解码

    3.使用交叉验证测量预测分数

    4.检查模型权重

    1

    下载数据

     
    from nilearn import datasets
    import warnings
    warnings.filterwarnings("ignore")
    """
    调用fetch_haxby()下载Haxby研究数据集,
    如果数据集已经存在本地了,就直接加载本地的,
    否则会从网上下载
    """
    haxby_dataset = datasets.fetch_haxby()
    #由于数据集有很多组,我们只是用其中一组
    fmri_filename = haxby_dataset.func[0]
    # 打印数据集的基本信息
    print('First subject functional nifti images (4D) are at: %s' %
          fmri_filename)  # 4D data

    2

    将fMRI数据转换为数据矩阵

     
    """
    以受检者的解剖图像为背景对其进行可视化
    """
    mask_filename = haxby_dataset.mask_vt[0]
    from nilearn import plotting
    plotting.plot_roi(mask_filename, bg_img=haxby_dataset.anat[0],
                     cmap='Paired')

    b8984ebd6be82cd84f009614a64a59cf.png

    """
    利用 nilearn.input_data.NiftiMasker来提取mask上的fMRI数据,
    并将其转换为数据序列
    """
    from nilearn.input_data import NiftiMasker
    masker = NiftiMasker(mask_img=mask_filename, standardize=True)
    
    # 根据 filename检索二维数据
    fmri_masked = masker.fit_transform(fmri_filename)
    
    # 打印frmi_masked
    print(fmri_masked)

    d9aa59adea3a3cfce87e1ac811f366e2.png

    # 查看其形状
    print(fmri_masked.shape)

    (1452, 464)

    3

    加载行为标签

     

    行为标签存储在csv文件中,用空格分隔,这里用pandas将其排列成一个数组

    import pandas as pd
    # 加载行为信息
    behavioral = pd.read_csv(haxby_dataset.session_target[0], sep=" ")
    print(behavioral

    d84449a32df6e08902bf01bcf75b985a.png

    conditions = behavioral['labels']
    print(conditions

    99c4be5890d2049a54a0e6e7024ae842.png


    4

    只对猫和面孔脸部进行分析

     

    根据上面的labels可以看出,有多种条件,本次实验只需猫和脸部的数据

    """
    只需要face和cat数据
    """
    condition_mask = conditions.isin(['face', 'cat'])
    
    fmri_masked = fmri_masked[condition_mask]
    # 打印其形状
    print(fmri_masked.shape)

    (216, 464)

    conditions = conditions[condition_mask]
    print(conditions.shape
    (216,)


    5

    利用SVM进行解码

     

    我们使用scikit-learn机器学习工具对fmri_masked数据进行分析。

    这里使用线性核的支持向量机分类器。

    """
    创建线性核的SVM
    """
    from sklearn.svm import SVC
    svc = SVC(kernel='linear')
    print(svc

    41d152595962086844e0dcbfd9041450.png

    """
    训练数据
    训练集:fmri_masked
    标签为:conditions
    """
    svc.fit(fmri_masked, conditions)
    """
    预测数据
    对fmri_masked数据进行预测,
    得到其对应的标签
    """
    prediction = svc.predict(fmri_masked)
    print(prediction)

    d487051080e03ecd5069cd647c3a8c74.png


    e58fc0eb64210bdf486c5379d26edf7e.png

    文章仅用于学术交流,不用于商业行为,

    若有侵权及疑问,请后台留言,管理员即时删侵!

    893a5317c64301bf3e6c56ce9bd91c69.png

    更多阅读

    美国“脑计划”取得重大进展,从小鼠到猴子再到人类

    深圳大学梁臻博士提出EEGFuseNet高维脑电图

    混合无监督深度特征表征与融合模型及其在情绪识别中的应用

    用于情绪识别的生物信号数据集汇总

    临港实验室、上海脑科学与类脑研究中心 脑机接口平台联合招聘公告

    a930199233941a78084d1d623f9b53a3.png

    28871ee4c77a565d7cb9252a8a3e2677.png

    你的每一次在看,我都很在意!

    展开全文
  • from nilearn import plotting n_subjects = 4 # subjects to consider for group-sparse covariance (max: 40) def plot_matrices(cov, prec, title, labels): """Plot covariance and precision matrices, for a ...

    总代码

    import numpy as np
    from nilearn import plotting
    n_subjects = 4  # subjects to consider for group-sparse covariance (max: 40)
    def plot_matrices(cov, prec, title, labels):
        """Plot covariance and precision matrices, for a given processing. """
        prec = prec.copy()  # avoid side effects
        # Put zeros on the diagonal, for graph clarity.
        size = prec.shape[0]
        prec[list(range(size)), list(range(size))] = 0
        span = max(abs(prec.min()), abs(prec.max()))
        # Display covariance matrix
        plotting.plot_matrix(cov, cmap=plotting.cm.bwr,
                             vmin=-1, vmax=1, title="%s / covariance" % title,
                             labels=labels)
        # Display precision matrix
        plotting.plot_matrix(prec, cmap=plotting.cm.bwr,
                             vmin=-span, vmax=span, title="%s / precision" % title,
                             labels=labels)
    

    第一步 获取数据集

    from nilearn import datasets
    n_subjects = 4  
    msdl_atlas_dataset = datasets.fetch_atlas_msdl()
    rest_dataset = datasets.fetch_development_fmri(n_subjects=n_subjects)
    # print basic information on the dataset
    print('First subject functional nifti image (4D) is at: %s' %
          rest_dataset.func[0])  # 4D data
    

    n_subjects就是多少个被试
    在这里插入图片描述msdl_atlas_dataset就是地图的标签,格式是bunch格式。
    在这里插入图片描述
    rest_dataset功能像图像,另外还附带了混淆矩阵。

    提取区域信号

    from nilearn import input_data
    # A "memory" to avoid recomputation
    from joblib import Memory
    mem = Memory('nilearn_cache')
    masker = input_data.NiftiMapsMasker(
        msdl_atlas_dataset.maps, resampling_target="maps", detrend=True,
        high_variance_confounds=True, low_pass=None, high_pass=0.01,
        t_r=2, standardize=True, memory='nilearn_cache', memory_level=1,
        verbose=2)
    masker.fit()
    subject_time_series = []
    func_filenames = rest_dataset.func
    confound_filenames = rest_dataset.confounds
    for func_filename, confound_filename in zip(func_filenames,
                                                confound_filenames):
        print("Processing file %s" % func_filename)
        region_ts = masker.transform(func_filename,
                                     confounds=confound_filename)
        subject_time_series.append(region_ts)
    

    mem设置内存运行文件
    在这里插入图片描述
    nilearn.input_data.NiftiLabelsMasker和我们平常使用的mask(掩码)差不多,这里的设置的意思是:使用msdl_atlas_dataset.maps地图集。
    resampling_target{“data”, “mask”, “maps”, None}, optional.
    表示对输入的图像重采样,这里是变成和maps地图集一样的尺寸形状。它有多个选择,可以是data或者mask,默认是none。
    detrend去除线性趋势。
    high_variance_confounds,在指定的maps图像上计算高方差混淆并回归。
    standardize标准化。
    subject_time_series 空列表
    在这里插入图片描述
    func_filenames就是功能像,4D的
    在这里插入图片描述
    confound矩阵文件
    然后通过for循环,提取了每一个功能像的ROI的时间序列,形成被试的时间序列。
    在这里插入图片描述
    在这里插入图片描述
    这就是subject_time_series文件,可以看出原始文件是python原始数据格式list,然后它的形状是(4,168,39),4个被试,每个ROI有168个时间序列,一共有39个ROI。
    默默吐槽一点,python自带的list还是不如numpy好用,好像nilearn以后会自全部变成使用numpy格式,这样更方便点。

    计算群稀疏精度矩阵

    from nilearn.connectome import GroupSparseCovarianceCV
    gsc = GroupSparseCovarianceCV(verbose=2)
    gsc.fit(subject_time_series)
    try:
        from sklearn.covariance import GraphicalLassoCV
    except ImportError:
        # for Scitkit-Learn < v0.20.0
        from sklearn.covariance import GraphLassoCV as GraphicalLassoCV
    gl = GraphicalLassoCV(verbose=2)
    gl.fit(np.concatenate(subject_time_series))
    

    gsc就是 具有交叉验证的参数选择稀疏的逆协方差矩阵估算器。
    group 组 sparse 稀疏的 covariance 协方差,cv就是交叉验证。如果不指定cv,默认是3. verbose就是打印出处理内容的方式,没啥用,可以设置可以不设置。
    在这里插入图片描述
    np.concatenate就是连接、拼接函数,就是把3维数据拼接成二维数据,所以被试1有168个时间序列,被试2有168个时间序列…4个被试就是672个时间序列,行是时间序列,列是39个地图的ROI。
    下面这个是计算协方差矩阵估算器fit的用法:
    estimator.fit([time_series_1, time_series_2, …])
    estimator.covariances_[0] 协方差
    estimator.precisions_[0] 逆协方差

    在这里插入图片描述
    注意,这里有两个估算器,lasso那个其实是参考,真正要用的就是gsc估算器,这个估算器是可以输入三维数据的,不需要去降维。

    显示结果

    def plot_matrices(cov, prec, title, labels):
        """Plot covariance and precision matrices, for a given processing. """
        prec = prec.copy()  # avoid side effects
        # Put zeros on the diagonal, for graph clarity.
        size = prec.shape[0]
        prec[list(range(size)), list(range(size))] = 0
        span = max(abs(prec.min()), abs(prec.max()))
        # Display covariance matrix
        plotting.plot_matrix(cov, cmap=plotting.cm.bwr,
                             vmin=-1, vmax=1, title="%s / covariance" % title,
                             labels=labels)
        # Display precision matrix
        plotting.plot_matrix(prec, cmap=plotting.cm.bwr,
                             vmin=-span, vmax=span, title="%s / precision" % title,
                             labels=labels)
    
    atlas_img = msdl_atlas_dataset.maps
    atlas_region_coords = plotting.find_probabilistic_atlas_cut_coords(atlas_img)
    labels = msdl_atlas_dataset.labels
    plotting.plot_connectome(gl.covariance_,
                             atlas_region_coords, edge_threshold='90%',
                             title="Covariance",
                             display_mode="lzr")
    plotting.plot_connectome(-gl.precision_, atlas_region_coords,
                             edge_threshold='90%',
                             title="Sparse inverse covariance (GraphicalLasso)",
                             display_mode="lzr",
                             edge_vmax=.5, edge_vmin=-.5)
    plot_matrices(gl.covariance_, gl.precision_, "GraphicalLasso", labels)
    
    title = "GroupSparseCovariance"
    plotting.plot_connectome(-gsc.precisions_[..., 0],
                             atlas_region_coords, edge_threshold='90%',
                             title=title,
                             display_mode="lzr",
                             edge_vmax=.5, edge_vmin=-.5)
    plot_matrices(gsc.covariances_[..., 0],
                  gsc.precisions_[..., 0], title, labels)
    plotting.show()
    

    首先是定义了一个plot_matrices的函数,这个函数的功能就是画相关性矩阵的图。先看下面的代码。
    在这里插入图片描述
    atlas_img就是一个4D的地图集,前3维是XYZ坐标,第四维是ROI的label。
    在这里插入图片描述
    在这里插入图片描述
    altas_region_coords就是ROI的坐标,需要跟地图空间保持一致。
    在这里插入图片描述
    plotting.plot_connectome这个函数就是用来画功能连接图的,如下
    在这里插入图片描述

    nilearn.plotting.plot_connectome(adjacency_matrix, node_coords, node_color=‘auto’, node_size=50, edge_cmap=<matplotlib.colors.LinearSegmentedColormap object>, edge_vmin=None, edge_vmax=None, edge_threshold=None, output_file=None, display_mode=‘ortho’, figure=None, axes=None, title=None, annotate=True, black_bg=False, alpha=0.7, edge_kwargs=None, node_kwargs=None, colorbar=False)
    第一个参数adjacency_matrix是表示图的链接强度。矩阵可以是对称的,这将导致无向图,也可以是不对称的,这将导致有向图。在这里就是用的-gsc.precisions_[…, 0]。
    补充说明:sklearn.covariance.GraphicalLassoCV是稀疏的逆协方差估计器,来自sklearn。需要注意的是,一般使用它的负值,就是-gl.precision_,这样就得到了偏相关,也就是直接相关。这比full correlation更加的精确。而对于这里使用的组被试,就需要用GroupSparseCovarianceCV来计算逆协方差,计算出来的值和Lasso是一样的。所以它也需要负值,-gsc.precisions_[…, 0]。covariance_是普通的协方差矩阵值。
    当然,最好的还是建议使用tangent space correlation,就是切线空间的相关。
    在这里插入图片描述
    用法就是:
    1.先定义估算器,measure = ConnectivityMeasure(kind=‘tangent’)
    2.开始估算,connectivities = measure.fit([time_series_1, time_series_2, …]) 。注意这里是time_1, time_2,所以需要用concatenate函数降维。
    (注意,GraphicalLassoCV—gl是专门用于单个主题的,对于多个主题,就需要先数据降维,变成time_1, time_2…的格式。而GroupSparseCovarianceCV—gsc是专门用于组计算的,所以不需要降维,直接3维就可以计算。这里的切线空间相关,得先降维)
    3.组内ROI的连接强度,使用measure.mean_ 平均值。

    edge_threshold=‘90%’:如果是数字,则仅显示值大于 edge_threshold 的边。如果是字符串,则必须以百分号结尾,例如“25.3%”,并且只会显示 abs(value) 高于给定百分位数的边。

    plot_matrices就是画下面这种相关性图:
    在这里插入图片描述
    plot_matrices(cov, prec, title, labels):
    输入协方差矩阵,精确矩阵,标题,脑区标签

    先来看看gsc.precisions_和gl.covariance_
    在这里插入图片描述
    和我猜的一样,就是39个ROI两两做相关,但是因为gsc是3维去做相关的,所以出来的值也是三维的,所以就是(39,39,4).
    而同理,gl是二维的,它把3维数据展平成2维,所以出来的矩阵是二维的。但是有一个问题是,它出来的(39,39),对于同一个label,它是进行平均了嘛?这个函数manual上没有说,一般是这么做。
    在这里插入图片描述
    1-复制精确矩阵-prec;2-获取精确矩阵第一个维度的长度-size;3-先生成0-38的值合并成列表,随后XX; 4-计算精确矩阵的最小值和最大值,把它们取绝对值,然后取两者的最大值-span; 5-

    展开全文
  • MRI机器学习工具箱nilearn: masker

    千次阅读 2019-04-25 12:15:13
    masker 对象的概念 对于任何基于神经影像的研究来说, 第一步都是... nilearn中 masking data 本质上是将4D的fmri数据变形成2D(voxel * timepoints). 但是如何将4D数据转为2D的数据, 对于不同的问题来说, 要选择的脑...
  • Nilearn学习笔记1-Nilearn库介绍

    万次阅读 2016-12-05 15:07:48
    nilearn是一个将机器学习、模式识别、多变量分析等技术应用于神经影像数据的应用中,能完成多体素模式分析(MVPA:mutli-voxel pattern analysis)、解码、模型预测、构造功能连接、脑区分割、构造连接体等功能。...
  • Nilearn中的基本操作和查看

    千次阅读 2019-12-14 22:14:57
    目录Nilearn简介Nilearn操作第一步:查看数据第二步:平滑操作第三步:保存结果到文件中 本分享为脑机学习者Rose整理发表于公众号:脑机接口社区(微信号:Brain_Computer).QQ交流群:903290195 Nilearn简介 Nilearn...
  • 知乎:Nilearn 入门、简介 Windows下安装 第一步:安装依赖项。 使用分发包管理器安装或要求系统管理员安装以下包:ipython、scipy、scikit learn(有时称为sklearn或python sklearn)、joblib、matplotlib(有时...
  • Nilearn 小记

    2018-11-19 21:23:00
    4.绘制脑图像 4.1 绘图功能 当打开了太多图像而不关闭时,会出现如下问题: 每次调用绘图函数都会创建一个新图像。当在非交互式设置(例如脚本...from nilearn import plotting display = plotting.plot_stat_map(...
  • pip install nilearn 提取海马:215-218 from nilearn import image import numpy as np area_name = 'hipp' ROI_label = range(215,219) img = image.load_img('BN_Atlas_246_1mm.nii') data = img.get_fdata() # ...
  • Nilearn教程系列(4)-脑部地图集绘制

    千次阅读 2019-12-26 19:12:30
    Harvard-Oxford 脑部地图集 第一步:下载数据集 数据集文件格式为: HarvardOxford-cort-maxprob-thr25-2mm.nii.gz from nilearn import datasets import warnings warnings.filterwarnings("ignore") """ 构建数据集...
  • 对于这个问题,nilearn提供了群体估计器:nilearn.connectome.GroupSparseCovarianceCV estimator. 它的用法如下: Estimator =nilearn .connectome .GroupSparseCovarianceCV Estimator .fit ([timeseries1,...
  • 目录3D和4D niimgs:处理和可视化第一步:加载数据第二步:可视化可视化4D文件 本分享为脑机学习者Rose整理发表于公众号:脑机接口社区(微信号:...from nilearn import datasets import warnings warnings.filterw...
  • Nilearn学习笔记2-从FMRI数据到时间序列

    万次阅读 多人点赞 2017-01-01 11:31:00
    (1) nilearn.masking.compute_background_mask for brain images where the brain stands out of a constant background. This is typically the case when working on statistic maps output after a brain ...
  • 将CanICA应用在核磁共振图像 from nilearn.decomposition import CanICA canica = CanICA(n_components=20, smoothing_fwhm=6., memory="nilearn_cache", memory_level=2, threshold=3., verbose=10, random_state=0...
  • Nilearn学习笔记3-提取时间序列建立功能连接体

    万次阅读 热门讨论 2017-01-02 14:31:47
    nilearn库中,提供了两种从fmri数据中提取时间序列的方法,一种基于脑分区(Time-series from a brain parcellation or “MaxProb” atlas),一种基于概率图谱(Time-series from a probabilistic atlas)。...
  • nilearn里的ICA功能生成fMRI激活区域 一、加载图像 import nilearn from nilearn import datasets func_filenames = nilearn.image.load_img("data/4D/Ontario_sub94652_rest.nii") print(func_filenames.shape) ...
  • 使用nilearn里的Dictlearning函数实现激活区域显示 说明:Dictlearning提取的信号比ICA更干净,所以要注意控制步骤二中alpha参数的大小 一、加载数据 import nilearn from nilearn import datasets nilearn.image....
  • Nilearn开发日的网站,2021年5月 安装 此主题需要Jekyll 3.8,因此它与GitHub Pages兼容。 将此行添加到您的Jekyll网站的Gemfile : gem "bulma-clean-theme" 并将此行添加到您的Jekyll网站的_config.yml : ...
  • 上次介绍了有mask的nilearn计算方法。 独特的mask(或者又称为ROI,种子点)是论文的亮点,但是其实有一半的情况下,尤其是对我们学生来说,没有充足的文献积累和前瞻,只是为了能发文章毕业,因此这种时候不要mask...
  • nilearn:用于Python中NeuroImaging的机器学习
  • 本分享为脑机学习者Rose整理发表于...可视化Nilearn附带的HCP connectome工作台颜色图,可用于在表面上绘制大脑图像。 import numpy as np import matplotlib.pyplot as plt from nilearn.plotting.cm import _...
  • 应用程序网络 dPys / PyNets的brainlife.io包装器
  • <div><p>on macOS: <pre><code> >>> import matplotlib >>> assert matplotlib.get_backend() =... import nilearn.plotting ...<p>I have <code>backend : Qt5Agg\...nilearn/nilearn</p></div>
  • 前面学习完带ROI的mask的预测(入门教程),以及group mask的预测(anova-SVM),把这些前置内容融汇贯通后,就可以开始看我下面的个人理解了。 首先回顾下前面的内容: 对于机器学习的基础:机器学习,说白了,就是...

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 657
精华内容 262
关键字:

nilearn