In this notebook, I will be taking the (t/a)compcor and Cosine confounds output by fmriprep and attempt to non-aggressively denoise the fmriprep confounds using the melodic matrix and noise ICs output from ICA-AROMA.
EDIT: Since the non-aggressively denoised fmriprep confounds looked strange, I tried a simpler full regression, so that any variance explained by the noise components is removed from the other fmriprep confounds of interest.
%matplotlib inline
import matplotlib.pyplot as plt
from nipype.interfaces.fsl.utils import FilterRegressor
import numpy as np
import nibabel as nib
import pandas as pd
import csv
# name files
noise_idx_file = './.confound_regression/sub-controlGE154_ses-pre_task-flanker_bold_AROMAnoiseICs.csv'
melodic_mix_file = './.confound_regression/sub-controlGE154_ses-pre_task-flanker_bold_MELODICmix.tsv'
confounds_file = './.confound_regression/sub-controlGE154_ses-pre_task-flanker_bold_confounds.tsv'
# list the noise ics (correctly, e.g. subtract one to match 0 index)
with open(noise_idx_file, 'r') as f:
reader = csv.reader(f, delimiter=',')
noise_idx = list(reader)[0]
noise_idx = [int(x) - 1 for x in noise_idx]
melodic_mix = np.loadtxt(melodic_mix_file)
confounds = pd.read_csv(confounds_file, delimiter='\t')
# extract the compcor and Cosine columns
confound_columns = [col for col in confounds.columns
if 'CompCor' in col or 'Cosine' in col]
axis_labels = confound_columns + ['noise_ic'+str(x) for x in noise_idx]
confound_matrix = confounds.as_matrix(columns=confound_columns)
axis_labels
['tCompCor00', 'tCompCor01', 'tCompCor02', 'tCompCor03', 'tCompCor04', 'tCompCor05', 'aCompCor00', 'aCompCor01', 'aCompCor02', 'aCompCor03', 'aCompCor04', 'aCompCor05', 'Cosine00', 'Cosine01', 'Cosine02', 'Cosine03', 'Cosine04', 'Cosine05', 'Cosine06', 'Cosine07', 'noise_ic2', 'noise_ic3', 'noise_ic4', 'noise_ic6', 'noise_ic7', 'noise_ic8', 'noise_ic9', 'noise_ic11', 'noise_ic12', 'noise_ic13', 'noise_ic14', 'noise_ic15', 'noise_ic16', 'noise_ic19', 'noise_ic20', 'noise_ic21', 'noise_ic22', 'noise_ic24', 'noise_ic26', 'noise_ic28', 'noise_ic29', 'noise_ic30', 'noise_ic31', 'noise_ic33', 'noise_ic35', 'noise_ic36', 'noise_ic38', 'noise_ic40', 'noise_ic42', 'noise_ic44', 'noise_ic45', 'noise_ic47', 'noise_ic48', 'noise_ic53', 'noise_ic54', 'noise_ic55', 'noise_ic58', 'noise_ic62', 'noise_ic63']
# filter melodic mix to only see the noise ICs
noise_ics = melodic_mix.T[noise_idx]
print(confound_matrix.T.shape)
print(noise_ics.shape)
(20, 315) (39, 315)
before_correlation = np.corrcoef(confound_matrix.T, noise_ics, rowvar=True)
import seaborn as sns
with sns.axes_style("white"):
plt.figure(figsize=(18,15))
sns.heatmap(before_correlation,center=0,vmax=1, vmin=-1, xticklabels=axis_labels, yticklabels=axis_labels)
# need to add more zero voxels since fsl assumes there are more voxels than timepoints in the image
confound_img_data = np.expand_dims(np.expand_dims(np.hstack((np.hstack((confound_matrix, np.ones([315,1]))), np.zeros([315,314]))).T,0),0)
confound_img_data.shape
(1, 1, 335, 315)
nib.save(confound_img, './.confound_regression/confounds.nii.gz')
regfilt = FilterRegressor(in_file='./.confound_regression/confounds.nii.gz', design_file=melodic_mix_file, filter_columns=[x+1 for x in noise_idx], out_file='./.confound_regression/confounds_non-aggr_denoise.nii.gz')
regfilt.run()
<nipype.interfaces.base.InterfaceResult at 0x7f06fda0a438>
filt_confound_img = nib.load('./.confound_regression/confounds_non-aggr_denoise.nii.gz')
confound_filt_data = filt_confound_img.get_data()
confound_filt_matrix = confound_filt_data[0][0][:][:]
confound_filt_matrix.shape
(20, 315)
after_correlation = np.corrcoef(confound_filt_matrix, noise_ics, rowvar=True)
/home/james/bin/miniconda3/lib/python3.6/site-packages/numpy/lib/function_base.py:3162: RuntimeWarning: invalid value encountered in true_divide c /= stddev[:, None] /home/james/bin/miniconda3/lib/python3.6/site-packages/numpy/lib/function_base.py:3163: RuntimeWarning: invalid value encountered in true_divide c /= stddev[None, :]
with sns.axes_style("white"):
plt.figure(figsize=(18,15))
sns.heatmap(after_correlation,center=0,vmax=1, vmin=-1, xticklabels=axis_labels, yticklabels=axis_labels)
lets try just simply regressing the noise_ics from the confounds and returning the residuals
full_regression_matrix = np.dot(noise_ics.T, np.dot(confound_matrix.T, np.linalg.pinv(noise_ics)).T).T
confound_full_filt_matrix = confound_matrix.T - full_regression_matrix
full_regression_correlation = np.corrcoef(confound_full_filt_matrix, noise_ics, rowvar=True)
with sns.axes_style("white"):
plt.figure(figsize=(18,15))
sns.heatmap(full_regression_correlation,center=0,vmax=1, vmin=-1, xticklabels=axis_labels, yticklabels=axis_labels)
with sns.axes_style("white"):
plt.figure(figsize=(18,15))
sns.heatmap(full_regression_correlation-before_correlation,center=0,vmax=1, vmin=-1, xticklabels=axis_labels, yticklabels=axis_labels)