-
Notifications
You must be signed in to change notification settings - Fork 12
/
plot_bids_features.py
256 lines (218 loc) · 7.07 KB
/
plot_bids_features.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
"""
First level analysis of a complete BIDS dataset from openneuro
==============================================================
Full step-by-step example of fitting a :term:`GLM`
to perform a first level analysis in an openneuro :term:`BIDS` dataset.
We demonstrate how :term:`BIDS`
derivatives can be exploited to perform a simple one subject analysis with
minimal code. Details about the :term:`BIDS` standard are available at
`https://bids.neuroimaging.io/ <https://bids.neuroimaging.io/>`_.
We also demonstrate how to download individual groups of files from the
Openneuro s3 bucket.
More specifically:
1. Download an :term:`fMRI` :term:`BIDS` dataset
with derivatives from openneuro.
2. Extract first level model objects automatically
from the :term:`BIDS` dataset.
3. Demonstrate Quality assurance of Nilearn estimation against available FSL.
estimation in the openneuro dataset.
4. Display contrast plot and uncorrected first level statistics table report.
"""
from nilearn import plotting
# %%
# Fetch openneuro :term:`BIDS` dataset
# ------------------------------------
# We download one subject from the stopsignal task
# in the ds000030 V4 :term:`BIDS` dataset available in openneuro.
# This dataset contains the necessary information to run a statistical analysis
# using Nilearn. The dataset also contains statistical results from a previous
# FSL analysis that we can employ for comparison with the Nilearn estimation.
from nilearn.datasets import (
fetch_ds000030_urls,
fetch_openneuro_dataset,
select_from_index,
)
_, urls = fetch_ds000030_urls()
exclusion_patterns = [
"*group*",
"*phenotype*",
"*mriqc*",
"*parameter_plots*",
"*physio_plots*",
"*space-fsaverage*",
"*space-T1w*",
"*dwi*",
"*beh*",
"*task-bart*",
"*task-rest*",
"*task-scap*",
"*task-task*",
]
urls = select_from_index(
urls, exclusion_filters=exclusion_patterns, n_subjects=1
)
data_dir, _ = fetch_openneuro_dataset(urls=urls)
# %%
# Obtain FirstLevelModel objects automatically and fit arguments
# --------------------------------------------------------------
# From the dataset directory we automatically obtain FirstLevelModel objects
# with their subject_id filled from the :term:`BIDS` dataset.
# Moreover we obtain,
# for each model, the list of run images and their respective events and
# confound regressors. Those are inferred from the confounds.tsv files
# available in the :term:`BIDS` dataset.
# To get the first level models we have to specify the dataset directory,
# the task_label and the space_label as specified in the file names.
# We also have to provide the folder with the desired derivatives, that in this
# case were produced by the :term:`fMRIPrep` :term:`BIDS` app.
from nilearn.glm.first_level import first_level_from_bids
task_label = "stopsignal"
space_label = "MNI152NLin2009cAsym"
derivatives_folder = "derivatives/fmriprep"
(
models,
models_run_imgs,
models_events,
models_confounds,
) = first_level_from_bids(
data_dir,
task_label,
space_label,
smoothing_fwhm=5.0,
derivatives_folder=derivatives_folder,
n_jobs=2,
)
# %%
# Access the model and model arguments of the subject and process events.
model, imgs, events, confounds = (
models[0],
models_run_imgs[0],
models_events[0],
models_confounds[0],
)
subject = f"sub-{model.subject_label}"
model.minimize_memory = False # override default
from pathlib import Path
from nilearn.interfaces.fsl import get_design_from_fslmat
fsl_design_matrix_path = (
Path(data_dir)
/ "derivatives"
/ "task"
/ subject
/ "stopsignal.feat"
/ "design.mat"
)
design_matrix = get_design_from_fslmat(
fsl_design_matrix_path, column_names=None
)
# %%
# We identify the columns of the Go and StopSuccess conditions of the
# design matrix inferred from the FSL file, to use them later for contrast
# definition.
design_columns = [
f"cond_{int(i):02}" for i in range(len(design_matrix.columns))
]
design_columns[0] = "Go"
design_columns[4] = "StopSuccess"
design_matrix.columns = design_columns
# %%
# First level model estimation (one subject)
# ------------------------------------------
# We fit the first level model for one subject.
model.fit(imgs, design_matrices=[design_matrix])
# %%
# Then we compute the StopSuccess - Go contrast. We can use the column names
# of the design matrix.
z_map = model.compute_contrast("StopSuccess - Go")
# %%
# We show the agreement between the Nilearn estimation and the FSL estimation
# available in the dataset.
import nibabel as nib
fsl_z_map = nib.load(
Path(data_dir)
/ "derivatives"
/ "task"
/ subject
/ "stopsignal.feat"
/ "stats"
/ "zstat12.nii.gz"
)
import matplotlib.pyplot as plt
from scipy.stats import norm
plotting.plot_glass_brain(
z_map,
colorbar=True,
threshold=norm.isf(0.001),
title='Nilearn Z map of "StopSuccess - Go" (unc p<0.001)',
plot_abs=False,
display_mode="ortho",
)
plotting.plot_glass_brain(
fsl_z_map,
colorbar=True,
threshold=norm.isf(0.001),
title='FSL Z map of "StopSuccess - Go" (unc p<0.001)',
plot_abs=False,
display_mode="ortho",
)
plt.show()
from nilearn.plotting import plot_img_comparison
plot_img_comparison(
[z_map], [fsl_z_map], model.masker_, ref_label="Nilearn", src_label="FSL"
)
plt.show()
# %%
# Simple statistical report of thresholded contrast
# -------------------------------------------------
# We display the :term:`contrast` plot and table with cluster information.
from nilearn.plotting import plot_contrast_matrix
plot_contrast_matrix("StopSuccess - Go", design_matrix)
plotting.plot_glass_brain(
z_map,
colorbar=True,
threshold=norm.isf(0.001),
plot_abs=False,
display_mode="z",
figure=plt.figure(figsize=(4, 4)),
)
plt.show()
# %%
# We can get a latex table from a Pandas Dataframe for display and publication
# purposes
from nilearn.reporting import get_clusters_table
table = get_clusters_table(z_map, norm.isf(0.001), 10)
print(table.to_latex())
# %%
# Generating a report
# -------------------
# Using the computed FirstLevelModel and :term:`contrast` information,
# we can quickly create a summary report.
from nilearn.reporting import make_glm_report
report = make_glm_report(
model=model,
contrasts="StopSuccess - Go",
)
# %%
# We have several ways to access the report:
# report # This report can be viewed in a notebook
# report.open_in_browser()
# or we can save as an html file
from pathlib import Path
output_dir = Path.cwd() / "results" / "plot_bids_features"
output_dir.mkdir(exist_ok=True, parents=True)
# report.save_as_html(output_dir / 'report.html')
# %%
# Saving model outputs to disk
# ----------------------------
from nilearn.interfaces.bids import save_glm_to_bids
save_glm_to_bids(
model,
contrasts="StopSuccess - Go",
contrast_types={"StopSuccess - Go": "t"},
out_dir=output_dir / "derivatives" / "nilearn_glm",
prefix=f"{subject}_task-stopsignal",
)
# %%
# View the generated files
files = sorted(list((output_dir / "derivatives" / "nilearn_glm").glob("*")))
print("\n".join([str(x.relative_to(output_dir)) for x in files]))