1"""
2Seed-based connectivity on the surface
3=======================================
4
5The dataset that is a subset of the enhanced NKI Rockland sample
6(http://fcon_1000.projects.nitrc.org/indi/enhanced/, Nooner et al, 2012)
7
8Resting state fMRI scans (TR=645ms) of 102 subjects were preprocessed
9(https://github.com/fliem/nki_nilearn) and projected onto the Freesurfer
10fsaverage5 template (Dale et al, 1999, Fischl et al, 1999). For this example
11we use the time series of a single subject's left hemisphere.
12
13The Destrieux parcellation (Destrieux et al, 2010) in fsaverage5 space as
14distributed with Freesurfer is used to select a seed region in the posterior
15cingulate cortex.
16
17Functional connectivity of the seed region to all other cortical nodes in the
18same hemisphere is calculated using Pearson product-moment correlation
19coefficient.
20
21The :func:`nilearn.plotting.plot_surf_stat_map` function is used
22to plot the resulting statistical map on the (inflated) pial surface.
23
24See also :ref:`for a similar example but using volumetric input data
25<sphx_glr_auto_examples_03_connectivity_plot_seed_to_voxel_correlation.py>`.
26
27See :ref:`plotting` for more details on plotting tools.
28
29References
30----------
31
32Nooner et al, (2012). The NKI-Rockland Sample: A model for accelerating the
33pace of discovery science in psychiatry. Frontiers in neuroscience 6, 152.
34URL http://dx.doi.org/10.3389/fnins.2012.00152
35
36Dale et al, (1999). Cortical surface-based analysis.I. Segmentation and
37surface reconstruction. Neuroimage 9.
38URL http://dx.doi.org/10.1006/nimg.1998.0395
39
40Fischl et al, (1999). Cortical surface-based analysis. II: Inflation,
41flattening, and a surface-based coordinate system. Neuroimage 9.
42http://dx.doi.org/10.1006/nimg.1998.0396
43
44Destrieux et al, (2010). Automatic parcellation of human cortical gyri and
45sulci using standard anatomical nomenclature. NeuroImage, 53, 1.
46URL http://dx.doi.org/10.1016/j.neuroimage.2010.06.010.
47"""
48
49###############################################################################
50# Retrieving the data
51# -------------------
52
53# NKI resting state data from nilearn
54from nilearn import datasets
55
56nki_dataset = datasets.fetch_surf_nki_enhanced(n_subjects=1)
57
58# The nki dictionary contains file names for the data
59# of all downloaded subjects.
60print(('Resting state data of the first subjects on the '
61       'fsaverag5 surface left hemisphere is at: %s' %
62      nki_dataset['func_left'][0]))
63
64# Destrieux parcellation for left hemisphere in fsaverage5 space
65destrieux_atlas = datasets.fetch_atlas_surf_destrieux()
66parcellation = destrieux_atlas['map_left']
67labels = destrieux_atlas['labels']
68
69# Fsaverage5 surface template
70fsaverage = datasets.fetch_surf_fsaverage()
71
72# The fsaverage dataset contains file names pointing to
73# the file locations
74print('Fsaverage5 pial surface of left hemisphere is at: %s' %
75      fsaverage['pial_left'])
76print('Fsaverage5 inflated surface of left hemisphere is at: %s' %
77      fsaverage['infl_left'])
78print('Fsaverage5 sulcal depth map of left hemisphere is at: %s' %
79      fsaverage['sulc_left'])
80
81###############################################################################
82# Extracting the seed time series
83# --------------------------------
84
85# Load resting state time series from nilearn
86from nilearn import surface
87
88timeseries = surface.load_surf_data(nki_dataset['func_left'][0])
89
90# Extract seed region via label
91pcc_region = b'G_cingul-Post-dorsal'
92
93import numpy as np
94pcc_labels = np.where(parcellation == labels.index(pcc_region))[0]
95
96# Extract time series from seed region
97seed_timeseries = np.mean(timeseries[pcc_labels], axis=0)
98
99###############################################################################
100# Calculating seed-based functional connectivity
101# ----------------------------------------------
102
103# Calculate Pearson product-moment correlation coefficient between seed
104# time series and timeseries of all cortical nodes of the hemisphere
105from scipy import stats
106
107stat_map = np.zeros(timeseries.shape[0])
108for i in range(timeseries.shape[0]):
109    stat_map[i] = stats.pearsonr(seed_timeseries, timeseries[i])[0]
110
111# Re-mask previously masked nodes (medial wall)
112stat_map[np.where(np.mean(timeseries, axis=1) == 0)] = 0
113
114###############################################################################
115# Display ROI on surface
116
117# Transform ROI indices in ROI map
118pcc_map = np.zeros(parcellation.shape[0], dtype=int)
119pcc_map[pcc_labels] = 1
120
121from nilearn import plotting
122
123plotting.plot_surf_roi(fsaverage['pial_left'], roi_map=pcc_map,
124                       hemi='left', view='medial',
125                       bg_map=fsaverage['sulc_left'], bg_on_data=True,
126                       title='PCC Seed')
127
128###############################################################################
129# Display unthresholded stat map with a slightly dimmed background
130plotting.plot_surf_stat_map(fsaverage['pial_left'], stat_map=stat_map,
131                            hemi='left', view='medial', colorbar=True,
132                            bg_map=fsaverage['sulc_left'], bg_on_data=True,
133                            darkness=.3, title='Correlation map')
134
135###############################################################################
136# Many different options are available for plotting, for example thresholding,
137# or using custom colormaps
138plotting.plot_surf_stat_map(fsaverage['pial_left'], stat_map=stat_map,
139                            hemi='left', view='medial', colorbar=True,
140                            bg_map=fsaverage['sulc_left'], bg_on_data=True,
141                            cmap='Spectral', threshold=.5,
142                            title='Threshold and colormap')
143
144###############################################################################
145# Here the surface is plotted in a lateral view without a background map.
146# To capture 3D structure without depth information, the default is to plot a
147# half transparent surface.
148# Note that you can also control the transparency with a background map using
149# the alpha parameter.
150plotting.plot_surf_stat_map(fsaverage['pial_left'], stat_map=stat_map,
151                            hemi='left', view='lateral', colorbar=True,
152                            cmap='Spectral', threshold=.5,
153                            title='Plotting without background')
154
155###############################################################################
156# The plots can be saved to file, in which case the display is closed after
157# creating the figure
158plotting.plot_surf_stat_map(fsaverage['infl_left'], stat_map=stat_map,
159                            hemi='left', bg_map=fsaverage['sulc_left'],
160                            bg_on_data=True, threshold=.5, colorbar=True,
161                            output_file='plot_surf_stat_map.png')
162
163plotting.show()
164