<jupyter_start><jupyter_text>Welcome to Neuroimaging in PythonIf you are running this from the SciNet jupyter hub, then use the paths in the top half of the block below.If you are running this on a local laptop, edit the bottom of the first block with the paths to your desired direcoties..<jupyter_code># We are going to use a single subject from the ABIDE dataset.
import os

# # ## use these commands if you are running this on scinet
my_home = os.environ.get('HOME')
my_scratch = my_home.replace('/home/', '/scratch/')
tutorial_dir = os.path.join(my_scratch, 'ss2017_mripython')
my_nilearn_data = os.path.join(tutorial_dir, 'data','nilearn_data')
cc200_nii = os.path.join(os.path.dirname(os.getcwd()), 'rois', 'cc200_roi_atlas.nii.gz')

## if this directory doesn't exist - make it..
rois_dir = os.path.join(tutorial_dir, 'data', 'rois')
if not os.path.exists(rois_dir):
    os.makedirs(rois_dir)

# ## edit this section if you are running this one your own laptop
# tutorial_dir ='/home/edickie/Documents/ss2017_16pythonmri/'  ## point this to the location where you want the data
# my_nilearn_data = os.path.join(tutorial_dir, 'data', 'nilearn_data') 

# ## point this to the location of the craddock atlas you downloaded
# cc200_nii = os.path.join(tutorial_dir, 'data', 'rois', 'cc200_roi_atlas.nii.gz') 

os.path.dirname(os.getcwd())
# this makes the plots show up..
import matplotlib.pyplot as plt
%matplotlib inline<jupyter_output><empty_output><jupyter_text>Using nilearn data fetchers to grab our tutorial datasetWe are grabbing data from the ABIDE dataset. A data sharing initiative with over 1000 scans. Check out the references included in the abide object for more info. Thanks to the work of the proprocessed connectome project. A lot of the data processing work has been done for us! We are going to grab the preprocessed but unfilted resting stat MR data.<jupyter_code>from nilearn import datasets

# We specify the site and number of subjects we want to download
abide = datasets.fetch_abide_pcp(data_dir=my_nilearn_data,
                                 derivatives=['func_mean', 'func_preproc', 'func_mask'],
                                 quality_checked= True, 
                                 n_subjects = 4)

print(abide.description)
print('{}'.format(abide.description))
abide.func_preproc<jupyter_output><empty_output><jupyter_text>Here we will define the files we want to work with for the rest of the tutorial We need:    + some Regions of Interest (rois) to use as masks for our analysis    + a functional MR image    + (for plotting) The mean of our functional MR image    + (for masking) a mask defining the brain vs outside of head    + to decide which roi we will focus on<jupyter_code>rois = cc200_nii
func = abide.func_preproc[3]
mask = abide.func_mask[3]
mean_func = abide.func_mean[3]
roi = 174
from nilearn import plotting

plotting.plot_roi(cc200_nii)
plotting.plot_roi(cc200_nii, bg_img=mean_func, alpha = 0.5, title = "cc200 atlas")<jupyter_output><empty_output><jupyter_text>Using nibabel to load nifti images<jupyter_code># this is for interacting with NIFTI files
import nibabel as nib

func_nib = nib.load(func)
print(func_nib)
affine = func_nib.get_affine()
affine
affine = func_nib.affine
affine
header = func_nib.get_header()
print(header.keys())
# get the dimensions for the fMRI file
dims = func_nib.shape
dims
# use get data to extract the data from it
func_data = func_nib.get_data()
func_data.shape
# reshape to voxels * timepoints (4D --> 2D)
func_data = func_data.reshape(dims[0]*dims[1]*dims[2], dims[3])<jupyter_output><empty_output><jupyter_text>Checkling that our functional and mask files have the same dimensions...<jupyter_code># now do the same thing for rois
rois_data = nib.load(cc200_nii).get_data()

print(rois_data.shape)

if rois_data.shape[0:2] == dims[0:2]:
    print("ROIs and func file dimensions match, Hooray!!")
else:
    print("FAIL, they do not match")<jupyter_output><empty_output><jupyter_text>Using nilearn to sample the rois nifti file so that they match in dimensions<jupyter_code>from nilearn import image


    
resampled_cc200 = image.resample_to_img(cc200_nii, mean_func, interpolation = 'nearest')

plotting.plot_roi(resampled_cc200, mean_func)

resampled_cc200.to_filename(os.path.join(rois_dir, "resample_cc200.nii.gz"))
# now do the same thing for rois
new_rois = os.path.join(tutorial_dir, 'data', 'rois', "resample_cc200.nii.gz")
rois_data = nib.load(new_rois).get_data()

print(rois_data.shape)

if rois_data.shape[0:2] == dims[0:2]:
    print("ROIs and func file dimensions match, Hooray!!")
else:
    print("FAIL, they do not match")<jupyter_output><empty_output><jupyter_text>Calculating our own seed correlation map using numpy<jupyter_code>rois_data = rois_data.reshape(dims[0]*dims[1]*dims[2], 1)<jupyter_output><empty_output><jupyter_text>Calculating our mean timeseries<jupyter_code>import numpy as np

# get the seed time series
idx = np.where(rois_data == roi)[0]
ts = np.mean(func_data[idx, :], axis=0)

# look at the timeseries
plt.plot(ts)
plt.show()
len(ts) # looks like it's the right length (i.e., our axis is correct)<jupyter_output><empty_output><jupyter_text>Correlate the mean timeseries with the timeseries from every voxel in the func file<jupyter_code># make an output matrix
output = np.zeros(dims[0]*dims[1]*dims[2])

# correlate seed against all voxels
for i in range(dims[0]*dims[1]*dims[2]):
    output[i] = np.corrcoef(ts, func_data[i, :])[0][1]<jupyter_output><empty_output><jupyter_text>When will do a faster and with no "true divide" warning if we incorporate the brain mask<jupyter_code>## add a mask to get rid of the weird bits
mask_data = nib.load(mask).get_data()
mask_data = mask_data.reshape(dims[0]*dims[1]*dims[2], 1)
idx_mask = np.where(mask_data > 0)[0]

## recalculate the ts taking the mask into account
idx_masked = np.intersect1d(idx,idx_mask)
ts = np.mean(func_data[idx_masked, :], axis=0)

# make an output matrix
output = np.zeros(dims[0]*dims[1]*dims[2])

# correlate seed against all voxels
for i in np.arange(len(idx_mask)):
    output[idx_mask[i]] = np.corrcoef(ts, func_data[idx_mask[i], :])[0][1]<jupyter_output><empty_output><jupyter_text>Using nibabel to write our result out to a nifti file<jupyter_code>np.unique(rois_data)
# get back to 4D
output_filename = os.path.join(tutorial_dir, 'seed_correlation.nii.gz')
output = np.reshape(output, (dims[0], dims[1], dims[2], 1))

# write the results into a NIFTI file
output = nib.nifti1.Nifti1Image(output, affine)
output.to_filename(output_filename)<jupyter_output><empty_output><jupyter_text>nilearn's glass brain and statmap plots help us visualize our new result<jupyter_code>plotting.plot_glass_brain(output_filename, threshold = 0.5, 
                          colorbar = True)
plotting.plot_stat_map(output_filename, threshold = 0.3, 
                          colorbar = True)<jupyter_output><empty_output><jupyter_text>instead of calculating one mean timeseries, let's do it for the whole atlas<jupyter_code># create output matrix
output_ts = np.zeros((len(np.unique(rois_data))-1, dims[3]))

# rois_data_masked = np.multiply(rois_data,mask_data)
# load in all mean time series from atlas
for i, roi in enumerate(np.unique(rois_data)):
    if roi > 0:
        idx = np.where(rois_data == roi)[0]
        ts = np.mean(func_data[idx, :], axis=0)
        output_ts[i-1, :] = ts

# correlation matrix
cc200_correl = np.corrcoef(output_ts)

# view output
plt.imshow(cc200_correl, interpolation='nearest', cmap=plt.cm.RdBu_r)
plt.colorbar()<jupyter_output><empty_output><jupyter_text>The seaborn package has some really awesome option..uncluding and easy clustered heatmap plot<jupyter_code>import seaborn

seaborn.clustermap(cc200_correl)<jupyter_output><empty_output><jupyter_text>Getting fancier with nilearnnilearn was built by people doing machine learning of resting state fMRI data using sci-kit learn.It has A LOT of functionallity beyond what we have time to teach here. Including:+ image manipulation (smoothing, resampling)+ algorithms (ICA and dictionary learning) for parcellating your own dataset (i.e. defining your own ROIS) + timeseries extraction with many many options for filtering, detrending + multiple methods for covariance esimtaiton (i.e. partial correlation, covariance, lasso..)+ machine learning and cross validationThey also did a beautiful job of building a large set of ipython notebook tutorials to teach you about all these methods.Go to the [the nilean website](http://nilearn.github.io/index.html) for more info.Let's use another example to see some of the more complex things nilearn can do...We'll start by using the data fetcher to grab a different atlas. This one (dosenbach_2010) if a set of coordinates in the brain to sample from..<jupyter_code>from nilearn import datasets

dosenbach = datasets.fetch_coords_dosenbach_2010()
print(dosenbach.description)<jupyter_output><empty_output><jupyter_text>Using nilearn's sphere masker to extract the timeseriesnilearn has a built in function for extracting timeseries from functional files and doing a little extra signal processing at the same time!<jupyter_code>from nilearn import input_data

spheres_masker = input_data.NiftiSpheresMasker(
    seeds=dosenbach.rois, #the seeds are the dosenbach roi atlas
    smoothing_fwhm=4, radius=4.5, # set the radius of a sphere around the roi you want extracted
    standardize=True, # the time-series are centered and normed (mean 0, variance 1 in the time dimension)
    detrend=True, low_pass=0.1, high_pass=0.01, t_r=2.5) # additional signal cleaning and filtering params

timeseries = spheres_masker.fit_transform(func)

print("the shape of the timeseries is {}".format(timeseries.shape))<jupyter_output><empty_output><jupyter_text>Using nilearn's ConnectivityMeasure to calculate our correlation matrixAvalable options are “correlation”, “partial correlation”, “tangent”, “covariance”, “precision” or other utilites in sci-py could be plugged in ([see here for an example](http://nilearn.github.io/auto_examples/03_connectivity/plot_multi_subject_connectome.htmlsphx-glr-auto-examples-03-connectivity-plot-multi-subject-connectome-py))Let's do partial correlation this time<jupyter_code>from nilearn.connectome import ConnectivityMeasure
from nilearn import plotting

correlation_measure = ConnectivityMeasure(kind='covariance')
dosenbach_matrix = correlation_measure.fit_transform([timeseries])[0]
seaborn.clustermap(dosenbach_matrix)
from nilearn import plotting

coords = np.vstack((
    dosenbach.rois['x'],
    dosenbach.rois['y'],
    dosenbach.rois['z'],
)).T

plotting.plot_connectome(dosenbach_matrix, node_coords = coords, edge_threshold='98%')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.098256735340729, 'fraction_numerical': 0.015102078866411859, 'mean_word_length': 5.427801078490114, 'pattern_counts': {'":': 0, '<': 69, '<?xml version=': 0, '>': 72, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '17', 'dir': 'jupyter-structured-clean-dedup', 'id': '280124', 'license': 'MIT', 'n_tokens_mistral': 3752, 'n_tokens_neox': 3582, 'n_words': 1052, 'path': 'notebook/ss2017_16mripython.ipynb', 'repo_name': 'edickie/nibabel_nilearn_tutorial_2017', 'size': '10727'}