<jupyter_start><jupyter_text>Welcome to Neuroimaging in PythonFor a good intro to the jupyter notebook interface and general python check out the following web tutorials.https://miykael.github.io/nipype_tutorial/  Edit and run this chunck if you are using your own laptop<jupyter_code># We are going to use a single subject from the ABIDE dataset.
import os
## edit this section if you are running this one your own laptop
tutorial_dir ='/home/USERNAME/Documents/ss2018_pythonmri/'  ## point this to the location where you want the data
nilearn_data = os.path.join(tutorial_dir, 'data', 'nilearn_data')<jupyter_output><empty_output><jupyter_text>If you are running this on the SciNet Jyputer Hub Run the next section..If not you can skip down  if the "coss2018_mri" environment not availableRun the next block to link the built conda env into your own home. (You might need to restart the kernel)<jupyter_code>! mkdir -p ~/conda/envs
! ln -s /scinet/course/ss2018/3_bm/1_mripython/conda_envs/coss2018_mri ~/.conda/envs/coss2018_mri
# 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, 'coss2018_mripython')

# make the tutorial dir if it does not exist
if not os.path.exists(tutorial_dir):
    os.makedirs(tutorial_dir)

# we are all going to read from the data in the courses folder (on the burst buffer!)
nilearn_data = '/scinet/course/ss2018/3_bm/1_mripython/data/nilearn_data'

cc200_nii = os.path.join(os.path.dirname(os.getcwd()), 'rois', 'cc200_roi_atlas.nii.gz')<jupyter_output><empty_output><jupyter_text>make the output directories if they don't exist<jupyter_code>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)<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
adhd = datasets.fetch_adhd(data_dir=nilearn_data, 
                                 n_subjects = 3)

print(adhd.description)
print('{}'.format(adhd.description.decode()))<jupyter_output><empty_output><jupyter_text>Concept Question 1:Explore the adhd object. What is the path to the functional files?  Concept Question 2:How was the data processed?<jupyter_code>adhd.phenotypic[0]<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 = adhd.func[0]
func_confounds = adhd.confounds[0]
roi = 174<jupyter_output><empty_output><jupyter_text>remember this line..it is magik!!! and important for seeing plots<jupyter_code># this makes the plots show up..
import matplotlib.pyplot as plt
%matplotlib inline  
from nilearn import plotting
cc200_nii
plotting.plot_roi(cc200_nii)
from nilearn import image as img

mean_img = img.mean_img(func)

plotting.plot_epi(mean_img)
plotting.plot_roi(cc200_nii, bg_img=mean_img, 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)<jupyter_output><empty_output><jupyter_text>We extract the affine matrix (we will use this to make a new image_Note_: in older version of nibabel the syntax is `affine = func_nib.get_affine()`<jupyter_code>affine = func_nib.affine
affine
header = func_nib.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<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>resampled_cc200 = img.resample_to_img(cc200_nii, mean_img, interpolation = 'nearest')

plotting.plot_roi(resampled_cc200, mean_img)

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**Step 1:** To make our math easier, we will reshape our data from 4D to 2D**Now :**+ axis 0 (rows): space (all of the voxels stacked)+ axis 1 (columns): time (timepoints)<jupyter_code># reshape to voxels * timepoints (4D --> 2D)
func_data = func_data.reshape(dims[0]*dims[1]*dims[2], dims[3])
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
import seaborn as sns

# get the seed time series
idx = np.where(rois_data == roi)[0]

# look at the timeseries
sns.tsplot(func_data[idx, :], err_style="unit_traces")

len(ts) # looks like it's the right length (i.e., our axis is correct)
ts = np.mean(func_data[idx, :], axis=0)

# look at the timeseries
sns.tsplot(ts)

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 maskfirst lets look at our main img to make one..<jupyter_code>mean_img_reshaped = mean_img.get_data().reshape(dims[0]*dims[1]*dims[2],1)
np.unique(mean_img_reshaped)[1:5]
## add a mask to get rid of the weird bits
mask_data = np.where(mean_img_reshaped > 1)[0]
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># 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.1, 
                          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.decode())<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, confound = adhd.confounds[0]) # 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><jupyter_text>Now let's repeat the sphere's with confound cleaning<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, confounds=func_confounds)

print("the shape of the timeseries is {}".format(timeseries.shape))
correlation_measure = ConnectivityMeasure(kind='covariance')
dosenbach_matrix = correlation_measure.fit_transform([timeseries])[0]
seaborn.clustermap(dosenbach_matrix)
## now let's try cleaning the whole image with image.clean
# image clearning inside nilearn

cleaned_func = img.clean_img(func, confounds=func_confounds, low_pass=0.1, high_pass=0.01, t_r=2.0,)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.09681787406905891, 'fraction_numerical': 0.01647483638004965, 'mean_word_length': 5.796523517382413, 'pattern_counts': {'":': 0, '<': 93, '<?xml version=': 0, '>': 97, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '23', 'dir': 'jupyter-structured-clean-dedup', 'id': '127241', 'license': 'MIT', 'n_tokens_mistral': 4671, 'n_tokens_neox': 4450, 'n_words': 1326, 'path': 'notebook/mripython_tutorial.ipynb', 'repo_name': 'edickie/nibabel_nilearn_tutorial_2017', 'size': '13293'}