Commit 7f7b1ab6 authored by Sean Fitzgibbon's avatar Sean Fitzgibbon
Browse files

improved notes

parent 26fe0bef
% MIGP - MELODIC's Incremental Group-PCA
%% MIGP - MELODIC's Incremental Group-PCA
% Steve Smith, FMRIB, 2012-2013
% not yet finally tweaked/evaluated, or published - please do not pass on.
% adapted by Sean Fitz, 2020
%%%%%%%%%%%%%%%%%%%%%%%%%% USER OPTIONS
%% User Options (equivalent to melodic defaults)
INlist=dir('data/cobre/fmri_*.nii.gz'); % list of input 4D NIFTI standard space datasets
INmask='data/brain_mask.nii.gz'; % 3D NIFTI mask image
GO='matMIGP'; % output filename
dPCAint=550; % internal number of components - typically 2-4 times number of timepoints in each run (if you have enough RAM for that)
dPCAout=100; % number of eigenvectors to output - should be less than dPCAint and more than the final ICA dimensionality
INlist=dir('~/nilearn_data/cobre/fmri_*_smooth.nii.gz'); % list of input 4D NIFTI standard space datasets
INmask='~/nilearn_data/brain_mask.nii.gz'; % 3D NIFTI mask image
GO='matMIGP'; % output filename
dPCAint=299; % internal number of components - typically 2-4 times number of timepoints in each run (if you have enough RAM for that)
dPCAout=299; % number of eigenvectors to output - should be less than dPCAint and more than the final ICA dimensionality
sep_vn=false; % switch on separate variance nomalisation for each input dataset
%%%%%%%%%%%%%%%%%%%%%%%%%% END OF USER OPTIONS
%% Run MIGP
Nsub=length(INlist); [~,r]=sort(rand(Nsub,1)); % will process subjects in random order
mask=read_avw(INmask);
......@@ -27,26 +26,33 @@ for i = 1:Nsub
% read data
filename=[INlist(r(i)).folder, '/', INlist(r(i)).name];
fprintf('Reading data file %s\n', filename);
grot=double(read_avw(filename));
grot=reshape(grot,prod(masksize),size(grot,4));
% demean
fprintf('\tRemoving mean image\n');
grot=demean(grot(mask~=0,:)');
% var-norm
[uu,ss,vv]=ss_svds(grot,30);
vv(abs(vv)<2.3*std(vv(:)))=0;
stddevs=max(std(grot-uu*ss*vv'),0.001);
grot=grot./repmat(stddevs,size(grot,1),1);
% var-norm separately for each subject
if sep_vn==true
fprintf('\tNormalising by voxel-wise variance\r');
[uu,ss,vv]=ss_svds(grot,30);
vv(abs(vv)<2.3*std(vv(:)))=0;
stddevs=max(std(grot-uu*ss*vv'),0.001);
grot=grot./repmat(stddevs,size(grot,1),1);
end
if (i==1)
W=demean(grot); clear grot;
else
% concat
W=[W; demean(grot)]; clear grot;
% reduce W to dPCAint eigenvectors
if size(W,1)-10 > dPCAint
fprintf('\tReducing data matrix to a %i dimensional subspace\n', dPCAint);
[uu,dd]=eigs(W*W',dPCAint);
W=uu'*W;
clear uu;
......@@ -56,8 +62,11 @@ for i = 1:Nsub
end
%% Reshape and save
grot=zeros(prod(masksize),dPCAout);
grot(mask~=0,:)=W(1:dPCAout,:)'; grot=reshape(grot,[masksize ,dPCAout]);
fprintf('Save to %s.nii.gz\n', GO)
save_avw(grot,GO,'f',[1 1 1 1]);
system(sprintf('fslcpgeom %s %s -d',filename,GO));
......@@ -16,7 +16,7 @@
"* [Clean the data](#clean-the-data)\n",
"* [Run `melodic`](#run-melodic)\n",
"\n",
"Firstly we will import the necessary packages for this notebook: "
"Firstly we will import the necessary packages for this notebook: "
]
},
{
......@@ -69,7 +69,9 @@
"<a class=\"anchor\" id=\"download-the-data\"></a>\n",
"## Download the data\n",
"\n",
"We use a method from [`nilearn`](https://nilearn.github.io/index.html) called `fetch_cobre` to download the fMRI data:"
"Create a directory in the users home directory to store the downloaded data:\n",
"\n",
"`expanduser` will expand the `~` to the be users home directory:"
]
},
{
......@@ -78,13 +80,25 @@
"metadata": {},
"outputs": [],
"source": [
"# Create a directory in the users home directory to store the downloaded data\n",
"data_dir = op.expanduser('~/nilearn_data')\n",
"\n",
"if not op.exists(data_dir):\n",
" os.makedirs(data_dir)\n",
" \n",
"# Download the data (if not already downloaded)\n",
" os.makedirs(data_dir)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Download the data (if not already downloaded). We use a method from [`nilearn`](https://nilearn.github.io/index.html) called `fetch_cobre` to download the fMRI data:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"d = datasets.fetch_cobre(data_dir=data_dir) "
]
},
......@@ -158,7 +172,7 @@
"outputs": [],
"source": [
"# generate melodic command line string\n",
"melodic_cmd = f\"melodic -i {','.join(smooth)} --mask={op.join(data_dir, 'brain_mask.nii.gz')} -d 10 -v --Oall -o cobre.gica \"\n",
"melodic_cmd = f\"melodic -i {','.join(smooth)} --mask={op.join(data_dir, 'brain_mask.nii.gz')} -d 10 -v -o cobre.gica \"\n",
"print(melodic_cmd)\n",
"\n",
"# run melodic\n",
......@@ -183,6 +197,13 @@
"ics = nb.load('cobre.gica/melodic_IC.nii.gz')\n",
"fig = map_plot(ics)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
......
%% Cell type:markdown id: tags:
# Fetch Data
This notebook will download an open fMRI dataset (~50MB) for use in the MIGP demo. It also regresses confounds from the data and performs spatial smoothing with 10mm FWHM.
This data is a derivative from the COBRE sample found in the International Neuroimaging Data-sharing Initiative (http://fcon_1000.projects.nitrc.org/indi/retro/cobre.html), originally released under Creative Commons - Attribution Non-Commercial.
It comprises 10 preprocessed resting-state fMRI selected from 72 patients diagnosed with schizophrenia and 74 healthy controls (6mm isotropic, TR=2s, 150 volumes).
* [Download the data](#download-the-data)
* [Clean the data](#clean-the-data)
* [Run `melodic`](#run-melodic)
Firstly we will import the necessary packages for this notebook:
%% Cell type:code id: tags:
``` python
from nilearn import datasets
from nilearn import image
from nilearn import plotting
import nibabel as nb
import numpy as np
import os.path as op
import os
import glob
import matplotlib.pyplot as plt
```
%% Cell type:markdown id: tags:
This function will be used to plot ICs later:
%% Cell type:code id: tags:
``` python
def map_plot(d):
N = d.shape[-1]
fig, ax = plt.subplots(int(np.ceil((N/2))),2, figsize=(12, N))
for img, ax0 in zip(image.iter_img(d), ax.ravel()):
coord = plotting.find_xyz_cut_coords(img, activation_threshold=2.3)
plotting.plot_stat_map(img, cut_coords=coord, vmax=10, axes=ax0)
return fig
```
%% Cell type:markdown id: tags:
<a class="anchor" id="download-the-data"></a>
## Download the data
We use a method from [`nilearn`](https://nilearn.github.io/index.html) called `fetch_cobre` to download the fMRI data:
Create a directory in the users home directory to store the downloaded data:
`expanduser` will expand the `~` to the be users home directory:
%% Cell type:code id: tags:
``` python
# Create a directory in the users home directory to store the downloaded data
data_dir = op.expanduser('~/nilearn_data')
if not op.exists(data_dir):
os.makedirs(data_dir)
```
%% Cell type:markdown id: tags:
Download the data (if not already downloaded). We use a method from [`nilearn`](https://nilearn.github.io/index.html) called `fetch_cobre` to download the fMRI data:
%% Cell type:code id: tags:
# Download the data (if not already downloaded)
``` python
d = datasets.fetch_cobre(data_dir=data_dir)
```
%% Cell type:markdown id: tags:
<a class="anchor" id="clean-the-data"></a>
## Clean the data
We use methods from [`nilearn`](https://nilearn.github.io/index.html) to regress confounds from the data (```clean_img```) and to spatially smooth the data with a gaussian filter of 10mm FWHM (```smooth_img```):
%% Cell type:code id: tags:
``` python
# Create a list of filenames for cleaned and smoothed data
clean = [f.replace('.nii.gz', '_clean.nii.gz') for f in d.func]
smooth = [f.replace('.nii.gz', '_clean_smooth.nii.gz') for f in d.func]
# loop through each subject, regress confounds and smooth
for img, cleaned, smoothed, conf in zip(d.func, clean, smooth, d.confounds):
image.clean_img(img, confounds=conf).to_filename(cleaned)
image.smooth_img(img, 10).to_filename(smoothed)
```
%% Cell type:markdown id: tags:
To run ```melodic``` we will need a brain mask in MNI152 space at the same resolution as the fMRI. Here we use [`nilearn`](https://nilearn.github.io/index.html) methods to load the MNI152 mask (```load_mni152_brain_mask```), resample to the resolution of the fMRI (```resample_to_img```), and binarize (```math_img```):
%% Cell type:code id: tags:
``` python
# load a single fMRI dataset (func0)
func0 = nb.load(d.func[0].replace('.nii.gz', '_clean_smooth.nii.gz'))
# load MNI153 brainmask, resample to func0 resolution, binarize, and save to nifti
mask = datasets.load_mni152_brain_mask()
mask = image.resample_to_img(mask, func0)
mask = image.math_img('img > 0.5', img=mask)
mask.to_filename(op.join(data_dir, 'brain_mask.nii.gz'))
# plot brainmask to make sure it looks OK
disp = plotting.plot_anat(image.index_img(func0, 0))
disp.add_contours(mask, threshold=0.5)
```
%% Cell type:markdown id: tags:
<a class="anchor" id="run-melodic"></a>
### Run ```melodic```
Generate a command line string and run group ```melodic``` on the smoothed fMRI with a dimension of 10 components:
%% Cell type:code id: tags:
``` python
# generate melodic command line string
melodic_cmd = f"melodic -i {','.join(smooth)} --mask={op.join(data_dir, 'brain_mask.nii.gz')} -d 10 -v --Oall -o cobre.gica "
melodic_cmd = f"melodic -i {','.join(smooth)} --mask={op.join(data_dir, 'brain_mask.nii.gz')} -d 10 -v -o cobre.gica "
print(melodic_cmd)
# run melodic
! {melodic_cmd}
```
%% Cell type:markdown id: tags:
Now we can load and plot the group ICs generated by ```melodic```.
Hopefully you can see some familiar looking RSN spatial patterns:
%% Cell type:code id: tags:
``` python
ics = nb.load('cobre.gica/melodic_IC.nii.gz')
fig = map_plot(ics)
```
%% Cell type:code id: tags:
``` python
```
......
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Matlab MIGP\n",
"\n",
"This notebook will load the dimension reduced data from Matlab MIGP, run group ICA, and then plot the group ICs.\n",
"\n",
"* [Run `melodic`](#run-melodic)\n",
"* [Plot group ICs](#plot-group-ics)\n",
"\n",
"Firstly we will import the necessary packages for this notebook: "
]
},
{
"cell_type": "code",
"execution_count": null,
......@@ -10,7 +24,17 @@
"from nilearn import image\n",
"import nibabel as nb\n",
"import matplotlib.pyplot as plt\n",
"import numpy as np"
"import numpy as np\n",
"import os.path as op"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"It will be necessary to know the location where the data was stored so that we can load the brainmask. \n",
"\n",
"`expanduser` will expand the `~` to the be users home directory:"
]
},
{
......@@ -19,16 +43,28 @@
"metadata": {},
"outputs": [],
"source": [
"%%bash\n",
"data_dir = op.expanduser('~/nilearn_data')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a class=\"anchor\" id=\"run-melodic\"></a>\n",
"### Run ```melodic```\n",
"\n",
"melodic -i matMIGP.nii.gz \\\n",
"\t--mask=data/brain_mask.nii.gz \\\n",
"\t-d 20 \\\n",
"\t-v \\\n",
"\t--nobet \\\n",
"\t--disableMigp \\\n",
"\t--varnorm \\\n",
"\t-o matMIGP_dim20.ica"
"Generate a command line string and run group ```melodic``` on the Matlab MIGP dimension reduced data with a dimension of 10 components:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# generate melodic command line string\n",
"melodic_cmd = f\"melodic -i matMIGP.nii.gz --mask={op.join(data_dir, 'brain_mask.nii.gz')} -d 10 -v --nobet --disableMigp -o matmigp.gica\"\n",
"print(melodic_cmd)"
]
},
{
......@@ -37,15 +73,56 @@
"metadata": {},
"outputs": [],
"source": [
"ics = nb.load('matMIGP_dim20.ica/melodic_IC.nii.gz')\n",
"# run melodic\n",
"! {melodic_cmd}"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a class=\"anchor\" id=\"plot-group-ics\"></a>\n",
"### Plot group ICs\n",
"\n",
"Now we can load and plot the group ICs generated by ```melodic```.\n",
"\n",
"This function will be used to plot ICs:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def map_plot(d):\n",
"\n",
"N = ics.shape[-1]\n",
" N = d.shape[-1]\n",
"\n",
"fig, ax = plt.subplots(int(np.ceil((N/2))),2, figsize=(12, N))\n",
" fig, ax = plt.subplots(int(np.ceil((N/2))),2, figsize=(12, N))\n",
"\n",
"for img, ax0 in zip(image.iter_img(ics), ax.ravel()):\n",
" coord = plotting.find_xyz_cut_coords(img, activation_threshold=2.3)\n",
" plotting.plot_stat_map(img, cut_coords=coord, vmax=5, axes=ax0)"
" for img, ax0 in zip(image.iter_img(d), ax.ravel()):\n",
" coord = plotting.find_xyz_cut_coords(img, activation_threshold=3.5)\n",
" plotting.plot_stat_map(img, cut_coords=coord, vmax=10, axes=ax0)\n",
" \n",
" return fig"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Hopefully you can see some familiar looking RSN spatial patterns:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"ics = nb.load('matmigp.gica/melodic_IC.nii.gz')\n",
"fig = map_plot(ics)"
]
},
{
......
%% Cell type:markdown id: tags:
## Matlab MIGP
This notebook will load the dimension reduced data from Matlab MIGP, run group ICA, and then plot the group ICs.
* [Run `melodic`](#run-melodic)
* [Plot group ICs](#plot-group-ics)
Firstly we will import the necessary packages for this notebook:
%% Cell type:code id: tags:
``` python
from nilearn import plotting
from nilearn import image
import nibabel as nb
import matplotlib.pyplot as plt
import numpy as np
import os.path as op
```
%% Cell type:markdown id: tags:
It will be necessary to know the location where the data was stored so that we can load the brainmask.
`expanduser` will expand the `~` to the be users home directory:
%% Cell type:code id: tags:
``` python
data_dir = op.expanduser('~/nilearn_data')
```
%% Cell type:markdown id: tags:
<a class="anchor" id="run-melodic"></a>
### Run ```melodic```
Generate a command line string and run group ```melodic``` on the Matlab MIGP dimension reduced data with a dimension of 10 components:
%% Cell type:code id: tags:
``` python
%%bash
# generate melodic command line string
melodic_cmd = f"melodic -i matMIGP.nii.gz --mask={op.join(data_dir, 'brain_mask.nii.gz')} -d 10 -v --nobet --disableMigp -o matmigp.gica"
print(melodic_cmd)
```
%% Cell type:code id: tags:
melodic -i matMIGP.nii.gz \
--mask=data/brain_mask.nii.gz \
-d 20 \
-v \
--nobet \
--disableMigp \
--varnorm \
-o matMIGP_dim20.ica
``` python
# run melodic
! {melodic_cmd}
```
%% Cell type:markdown id: tags:
<a class="anchor" id="plot-group-ics"></a>
### Plot group ICs
Now we can load and plot the group ICs generated by ```melodic```.
This function will be used to plot ICs:
%% Cell type:code id: tags:
``` python
ics = nb.load('matMIGP_dim20.ica/melodic_IC.nii.gz')
def map_plot(d):
N = ics.shape[-1]
N = d.shape[-1]
fig, ax = plt.subplots(int(np.ceil((N/2))),2, figsize=(12, N))
fig, ax = plt.subplots(int(np.ceil((N/2))),2, figsize=(12, N))
for img, ax0 in zip(image.iter_img(ics), ax.ravel()):
coord = plotting.find_xyz_cut_coords(img, activation_threshold=2.3)
plotting.plot_stat_map(img, cut_coords=coord, vmax=5, axes=ax0)
for img, ax0 in zip(image.iter_img(d), ax.ravel()):
coord = plotting.find_xyz_cut_coords(img, activation_threshold=3.5)
plotting.plot_stat_map(img, cut_coords=coord, vmax=10, axes=ax0)
return fig
```
%% Cell type:markdown id: tags:
Hopefully you can see some familiar looking RSN spatial patterns:
%% Cell type:code id: tags:
``` python
ics = nb.load('matmigp.gica/melodic_IC.nii.gz')
fig = map_plot(ics)
```
%% Cell type:code id: tags:
``` python
```
......
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Python MIGP\n",
"\n",
"This notebook will perform *python* MIGP dimension reduction, run group ICA, and then plot the group ICs.\n",
"\n",
"* [Run python `MIGP`](#run-python-migp)\n",
"* [Run `melodic`](#run-melodic)\n",
"* [Plot group ICs](#plot-group-ics)\n",
"\n",
"Firstly we will import the necessary packages for this notebook: "
]
},
{
"cell_type": "code",
"execution_count": null,
......@@ -13,7 +28,17 @@
"from scipy.sparse.linalg import svds, eigs\n",
"import matplotlib.pyplot as plt\n",
"from nilearn import plotting\n",
"from nilearn import image"
"from nilearn import image\n",
"import os.path as op"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"It will be necessary to know the location where the data was stored so that we can load the brainmask. \n",
"\n",
"`expanduser` will expand the `~` to the be users home directory:"
]
},
{
......@@ -22,11 +47,17 @@
"metadata": {},
"outputs": [],
"source": [
"in_list = [nb.load(f) for f in glob.glob('data/cobre/fmri_*.nii.gz')]\n",
"in_mask = nb.load('data/brain_mask.nii.gz')\n",
"go = '3DMIGP'\n",
"dPCA_int = 550\n",
"dPCA_out = 100"
"data_dir = op.expanduser('~/nilearn_data')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a class=\"anchor\" id=\"run-python-migp\"></a>\n",
"### Run python `MIGP`\n",
"\n",
"Firstly we need to set the MIGP parameters:"
]
},
{
......@@ -35,7 +66,15 @@
"metadata": {},
"outputs": [],
"source": [
"random.shuffle(in_list)"
"# create lists of (nibabel) image objects\n",
"in_list = [nb.load(f) for f in glob.glob(f'{data_dir}/cobre/fmri_*_smooth.nii.gz')]\n",
"in_mask = nb.load(f'{data_dir}/brain_mask.nii.gz')\n",
"\n",
"# set user parameters (equivalent to melodic defaults)\n",
"GO = 'pyMIGP.nii.gz' # output filename\n",
"dPCA_int = 299 # internal number of components - typically 2-4 times number of timepoints in each run (if you have enough RAM for that)\n",
"dPCA_out = 299 # number of eigenvectors to output - should be less than dPCAint and more than the final ICA dimensionality\n",
"sep_vn = False # switch on separate variance nomalisation for each input dataset"
]
},
{
......@@ -44,26 +83,36 @@
"metadata": {},
"outputs": [],
"source": [
"# randomise the subject order\n",
"random.shuffle(in_list)\n",
"\n",
"# load and unravel brainmask\n",
"mask = in_mask.get_fdata().ravel()\n",
"\n",
"# function to demean the data\n",
"def demean(x):\n",
" return x - np.mean(x, axis=0)\n",
"\n",
"# loop through input files/subjects\n",
"for i, f in enumerate(in_list):\n",
" \n",
" print(i, end=' ')\n",
" \n",
" # read data\n",
" print(f'Reading data file {f.get_filename()}')\n",
" grot = f.get_fdata()\n",
" grot = np.reshape(grot, [-1, grot.shape[-1]])\n",
" grot = grot[mask!=0, :].T\n",
" \n",
" # demean\n",
" print(f'\\tRemoving mean image')\n",
" grot = demean(grot)\n",
" \n",
" # var-norm\n",
" [uu, ss, vt] = svds(grot, k=30)\n",
" vt[np.abs(vt) < (2.3 * np.std(vt))] = 0;\n",
" stddevs = np.maximum(np.std(grot - (uu @ np.diag(ss) @ vt), axis=0), 0.001)\n",
" grot = grot/stddevs\n",
" if sep_vn:\n",
" print(f'\\tNormalising by voxel-wise variance')\n",
" [uu, ss, vt] = svds(grot, k=30)\n",
" vt[np.abs(vt) < (2.3 * np.std(vt))] = 0;\n",
" stddevs = np.maximum(np.std(grot - (uu @ np.diag(ss) @ vt), axis=0), 0.001)\n",
" grot = grot/stddevs\n",
" \n",
" if i == 0:\n",
" W = demean(grot)\n",
......@@ -73,23 +122,28 @@
" \n",
" # reduce W to dPCA_int eigenvectors\n",
" if W.shape[0]-10 > dPCA_int:\n",
" print(f'\\tReducing data matrix to a {dPCA_int} dimensional subspace')\n",
" uu = eigs(W@W.T, dPCA_int)[1]\n",
" uu = np.real(uu)\n",
" W = uu.T @ W\n",
" \n",
" f.uncache()\n",
" \n"
"# reshape and save\n",
"grot = np.zeros([mask.shape[0], dPCA_out])\n",
"grot[mask!=0, :] = W[:dPCA_out, :].T\n",
"grot = np.reshape(grot, in_list[0].shape[:3] + (dPCA_out,))\n",
"\n",
"print(f'Save to {GO}')\n",
"nb.Nifti1Image(grot, affine=in_list[0].affine).to_filename(GO)\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"cell_type": "markdown",
"metadata": {},
"outputs": [],
"source": [
"grot = np.zeros([mask.shape[0], dPCA_out])\n",
"grot[mask!=0, :] = W[:dPCA_out, :].T\n",
"grot = np.reshape(grot, in_list[0].shape[:3] + (dPCA_out,))"
"<a class=\"anchor\" id=\"run-melodic\"></a>\n",
"### Run ```melodic```\n",
"\n",
"Generate a command line string and run group ```melodic``` on the Python MIGP dimension reduced data with a dimension of 10 components:"
]
},
{
......@@ -98,7 +152,9 @@
"metadata": {},
"outputs": [],
"source": [
"nb.Nifti1Image(grot, affine=in_list[0].affine, header=in_list[0].header).to_filename('pyMIGP.nii.gz')"
"# generate melodic command line string\n",
"melodic_cmd = f\"melodic -i pyMIGP.nii.gz --mask={op.join(data_dir, 'brain_mask.nii.gz')} -d 10 -v --nobet --disableMigp -o pymigp.gica\"\n",
"print(melodic_cmd)"
]
},
{
......@@ -107,16 +163,20 @@
"metadata": {},
"outputs": [],
"source": [
"%%bash\n",
"\n",
"melodic -i pyMIGP.nii.gz \\\n",
"\t--mask=data/brain_mask.nii.gz \\\n",
"\t-d 20 \\\n",
"\t-v \\\n",
"\t--nobet \\\n",
"\t--disableMigp \\\n",
"\t--varnorm \\\n",
"\t-o pymigp_dim20.ica"
"# run melodic\n",
"! {melodic_cmd}"
]