diff --git a/lost+found/02_asymmetry-report.py b/lost+found/02_asymmetry-report.py
deleted file mode 100644
index e0614f9..0000000
--- a/lost+found/02_asymmetry-report.py
+++ /dev/null
@@ -1,272 +0,0 @@
-from xhtml2pdf import pisa
-import sys
-import pandas as pd
-import os
-import argparse
-
-# defined command line options
-CLI=argparse.ArgumentParser()
-CLI.add_argument("subject", nargs=1, type=str, default="")
-CLI.add_argument("session", nargs=1, type=str, default="")
-CLI.add_argument("out", nargs=1, type=str, default="")
-CLI.add_argument("demo", nargs=1, type=str, default="")
-CLI.add_argument("thr", nargs=1, type=float, default=1.96)
-CLI.add_argument(
- "--featList_ctx",
- nargs="*",
- type=str,
- default=['flair', 'qt1', 'adc', 'thickness'])
-CLI.add_argument(
- "--featList_sctx",
- nargs="*",
- type=str,
- default=[])
-CLI.add_argument(
- "--featList_hipp",
- nargs="*",
- type=str,
- default=[])
-
-# parse the command line
-args = CLI.parse_args()
-
-# access CLI options
-subject = args.subject[0]
-session = args.session[0]
-out = args.out[0]
-demo = args.demo[0]
-thr = args.thr[0]
-featList_ctx = ', '.join(args.featList_ctx)
-featList_sctx = ', '.join(args.featList_sctx)
-featList_hipp = ', '.join(args.featList_hipp)
-
-TBL = pd.read_excel(demo, engine='openpyxl', skiprows=[0,2]).dropna(how="all").dropna(axis=1, how="all")
-
-# Get participant's sex
-if TBL[TBL['ID'] == subject]['sex'].to_list()[0] == 'F':
- sex = 'Female'
-elif TBL[TBL['ID'] == subject]['sex'].to_list()[0] == 'M':
- sex = 'Male'
-else:
- sex = 'Not defined'
-
-# Get participant's age
-try:
- age = TBL[TBL['ID'] == subject]['age'].to_list()[0]
-except:
- age = 'Not defined'
-
-# Path to figures
-figPath = os.path.join(os.path.dirname(out), "analysis", "asymmetry", subject, session) + "/" + subject
-
-def report_block_template(figPath, subject='', session='', sex='', age='', featList_ctx='', featList_sctx='', featList_hipp=''):
- ctx_block = (' '
- '
' )
- sctx_block = (' '
- '
'
- '
')
- hipp_block = (
- ' '
- '
'
- '
')
-
- # FLAIR
- if 'flair' in featList_ctx:
- flair_ctx_block = (' '
- '
' )
- else:
- flair_ctx_block = ('No cortical FLAIR
')
-
- if 'flair' in featList_sctx:
- flair_sctx_block = (' '
- '
' )
- else:
- flair_sctx_block = ('No subcortical FLAIR
')
-
- if 'flair' in featList_hipp:
- flair_hipp_block = (' '
- '
' )
- else:
- flair_hipp_block = ('No hippocampal FLAIR
')
-
- # qT1
- if 'qt1' in featList_ctx:
- qt1_ctx_block = (' '
- '
' )
- else:
- qt1_ctx_block = ('No cortical qT1
')
-
- if 'qt1' in featList_sctx:
- qt1_sctx_block = (' '
- '
' )
- else:
- qt1_sctx_block = ('No subcortical qT1
')
-
- if 'qt1' in featList_hipp:
- qt1_hipp_block = (' '
- '
' )
- else:
- qt1_hipp_block = ('No hippocampal qT1
')
-
- # ADC
- if 'adc' in featList_ctx:
- adc_ctx_block = (' '
- '
' )
- else:
- adc_ctx_block = ('No cortical ADC
')
-
- if 'adc' in featList_sctx:
- adc_sctx_block = (' '
- '
' )
- else:
- adc_sctx_block = ('No subcortical ADC
')
-
- if 'adc' in featList_hipp:
- adc_hipp_block = (' '
- '
' )
- else:
- adc_hipp_block = ('No hippocampal ADC
')
-
- # THICKNESS
- if 'thickness' in featList_ctx:
- thickness_ctx_block = (' '
- '
' )
- else:
- thickness_ctx_block = ('No cortical atrophy
')
-
- if 'thickness' in featList_sctx:
- thickness_sctx_block = (' '
- '
' )
- else:
- thickness_sctx_block = ('No subcortical atrophy
')
-
- if 'thickness' in featList_hipp:
- thickness_hipp_block = (' '
- '
' )
- else:
- thickness_hipp_block = ('No hippocampal atrophy
')
-
- report_block = (
- # =============================================================================================================#
- # =============================================== P A G E # 1 ===============================================#
- # =============================================================================================================#
- # Title
- '
Clinical report | Asymmetry
'
-
- # Subject's ID - Session - Basic demographics
- ' Subject: {subject}, Session: {session}, Sex: {sex}, Age: {age}, '
- 'z-score threshold: {thr}
'
- ' Blue = Right > Left
'
- ' Red = Left > Right
' + '
' +
-
- # Multivariate cortical asymmetry
- ' Multivariate cortical asymmetry ({featList_ctx})
' + ctx_block +
- ''
-
- ' Multivariate subcortical ({featList_sctx})
' + sctx_block + ''
-
- ' Multivariate hippocampal asymmetry ({featList_hipp})
' +
- hipp_block + ''
-
- + '
'
-
- # =============================================================================================================#
- # =============================================== P A G E # 2 ===============================================#
- # ================================================= F L A I R =================================================#
- # =============================================================================================================#
- + '
'
- ' Univariate cortical FLAIR
' + flair_ctx_block + ''
-
- + ' Univariate subcortical FLAIR
'
- + flair_sctx_block + ''
-
- + ' Univariate hippocampal FLAIR
'
- + flair_hipp_block + '' + '
'
-
- # =============================================================================================================#
- # =============================================== P A G E # 3 ===============================================#
- # =================================================== q T 1 ===================================================#
- # =============================================================================================================#
- + '
'
- ' Univariate cortical qT1
' + qt1_ctx_block + ''
-
- + ' Univariate subcortical qT1
'
- + qt1_sctx_block + ''
-
- + ' Univariate hippocampal qT1
'
- + qt1_hipp_block + '' + '
'
-
- # =============================================================================================================#
- # =============================================== P A G E # 4 ===============================================#
- # =================================================== A D C ===================================================#
- # =============================================================================================================#
- + '
'
- ' Univariate cortical ADC
' + adc_ctx_block + ''
-
- + ' Univariate subcortical ADC
'
- + adc_sctx_block + ''
-
- + ' Univariate hippocampal ADC
'
- + adc_hipp_block + '' + '
'
-
- # =============================================================================================================#
- # =============================================== P A G E # 5 ===============================================#
- # ============================================= T H I C K N E S S =============================================#
- # =============================================================================================================#
- + '
'
- ' Univariate cortical atrophy
' + thickness_ctx_block + ''
-
- + ' Univariate subcortical atrophy
'
- + thickness_sctx_block + ''
-
- + ' Univariate hippocampal atrophy
'
- + thickness_hipp_block + '' + '
'
-
- )
-
- return report_block.format(figPath=figPath, subject=subject, session=session, sex=sex, age=age, thr=thr,
- featList_ctx=featList_ctx, featList_sctx=featList_sctx, featList_hipp=featList_hipp)
-
-# Create PDF report
-static_report = ''
-_static_block = report_block_template(figPath, subject=subject, session=session, sex=sex, age=age,
- featList_ctx=featList_ctx, featList_sctx=featList_sctx, featList_hipp=featList_hipp)
-static_report += _static_block
-
-# Utility function
-def convert_html_to_pdf(source_html, output_filename):
- # open output file for writing (truncated binary)
- result_file = open(output_filename, "w+b")
-
- # convert HTML to PDF
- pisa_status = pisa.CreatePDF(
- source_html, # the HTML to convert
- dest=result_file) # file handle to recieve result
-
- # close output file
- result_file.close() # close output file
-
- # return True on success and False on errors
- return pisa_status.err
-
-convert_html_to_pdf(static_report, figPath + '_asymmetryReport.pdf')
-
diff --git a/lost+found/02_asymmetry.py b/lost+found/02_asymmetry.py
deleted file mode 100644
index e1fdd82..0000000
--- a/lost+found/02_asymmetry.py
+++ /dev/null
@@ -1,462 +0,0 @@
-import os
-import sys
-import pandas as pd
-import numpy as np
-import nibabel as nib
-import glob
-from enigmatoolbox.utils.useful import zscore_matrix
-from enigmatoolbox.datasets import load_mask
-from enigmatoolbox.plotting import plot_cortical, plot_subcortical, plot_hippocampal
-from enigmatoolbox.datasets import getaffine
-import warnings
-
-warnings.simplefilter('ignore')
-
-import argparse
-# defined command line options
-# this also generates --help and error handling
-CLI=argparse.ArgumentParser()
-CLI.add_argument("subject", nargs=1, type=str, default="")
-CLI.add_argument("session", nargs=1, type=str, default="")
-CLI.add_argument("out", nargs=1, type=str, default="")
-CLI.add_argument("hipDir", nargs=1, type=str, default="")
-CLI.add_argument("demo", nargs=1, type=str, default="")
-CLI.add_argument("thr", nargs=1, type=float, default=1.96)
-CLI.add_argument(
- "--featList_ctx", # name on the CLI - drop the `--` for positional/required parameters
- nargs="*",
- type=str,
- default=['flair', 'qt1', 'adc', 'thickness'])
-CLI.add_argument(
- "--featList_sctx",
- nargs="*",
- type=str, # any type/callable can be used here
- default=[])
-CLI.add_argument(
- "--featList_hipp",
- nargs="*",
- type=str, # any type/callable can be used here
- default=[])
-
-# parse the command line
-args = CLI.parse_args()
-
-# access CLI options
-subject = args.subject[0]
-session = args.session[0]
-out = args.out[0]
-hipDir = args.hipDir[0]
-demo = args.demo[0]
-thr = args.thr[0]
-featList_ctx = args.featList_ctx
-featList_sctx = args.featList_sctx
-featList_hipp = args.featList_hipp
-
-def new_gifti_image(data, intent=0, datatype=16, metadata=None):
- """NiBabel wrapper to generate a gifti image with data array and metadata.
- Parameters
- ----------
- data : ndarray
- 1-D ndarray containing one hemisphere surface data.
- intent : int
- Intent code for Gifti File. Defaults to 0 (Intent = NONE).
- Available intent codes:
- NIFTI_INTENT_NONE - 0
- NIFTI_INTENT_CORREL - 2
- NIFTI_INTENT_TTEST - 3
- NIFTI_INTENT_ZSCORE - 5
- NIFTI_INTENT_PVAL - 22
- NIFTI_INTENT_LOGPVAL - 23
- NIFTI_INTENT_LOG10PVAL - 24
- NIFTI_INTENT_LABEL - 1002
- NIFTI_INTENT_POINTSET - 1008
- NIFTI_INTENT_TRIANGLE - 1009
- NIFTI_INTENT_TIME_SERIES - 2001
- NIFTI_INTENT_NODE_INDEX - 2002
- NIFTI_INTENT_SHAPE - 2005
- More intent codes can be found at: https://nifti.nimh.nih.gov/nifti-1/documentation/nifti1fields/nifti1fields_pages/group__NIFTI1__INTENT__CODES.html
- datatype : int
- Datatype for gifti image. Defaults to 16 (dtype = float32)
- Available datatypes:
- UINT8 - 2
- INT32 - 8
- FLOAT32 - 16
- metadata : nibabel.gifti.gifti.GiftiMetaData
- Metadata for gifti image.
- Returns
- -------
- nibabel.gifti.gifti.GiftiImage
- Gifti image with specified metadata and data array.
- """
- dtypes = {2: np.uint8, 8: np.int32, 16: np.float32}
- data = data.astype(dtypes[datatype])
- if metadata:
- metadata = nib.gifti.GiftiMetaData(metadata)
- gifti_data = nib.gifti.GiftiDataArray(data=data, intent=intent, datatype=datatype)
- gifti_img = nib.gifti.GiftiImage(meta=metadata, darrays=[gifti_data])
- return gifti_img
-
-
-def zscore(mtrx, TBL):
- grpk = TBL['grp']
- data_z = zscore_matrix(mtrx, grpk, 'HC')
- return data_z
-
-
-def matrix(area, path, filename, TBL):
- # Get number of vertices/structures
- if area == "ctx":
- nvert = 64984
- elif area == "sctx":
- nvert = 14
- elif area == "hipp":
- nvert = 14524
-
- # Subject IDs
- sub = TBL['ID']
- (nSub,) = sub.shape
-
- # Initialize matrix
- mtrx = np.empty((nSub, nvert))
- mtrx[:] = np.nan
-
- # Cortical feature matrix
- if area == "ctx":
- for i in range(nSub):
- try:
- dpath = out + "/" + sub[i] + "/" + session + "/" + path
- d = []
- for _, h in enumerate(['lh', 'rh']):
- d = np.append(d, nib.load(dpath + sub[i] + '_' + session +
- filename.format(h)).get_fdata().squeeze())
- mtrx[i, :] = d
- except:
- print(sub[i] + " : NO " + filename)
- pass
-
- # Subcortical feature matrix
- elif area == "sctx":
- for i in range(nSub):
- try:
- dpath = out + "/" + sub[i] + "/" + session + "/" + path
- mtrx[i] = np.loadtxt(dpath + sub[i] + "_" + session + filename,
- delimiter=",", skiprows=1, usecols=range(1,15))
- except:
- print(sub[i] + " : NO " + filename)
- pass
-
- # Hippocampal feature matrix
- elif area == "hipp" and "_thickness" in filename:
- for i in range(nSub):
- try:
- dpath = hipDir + "/" + sub[i] + "/" + path
- d = []
- for _, h in enumerate(['L', 'R']):
- d = np.append(d, nib.load(dpath + sub[i] +
- filename.format(h)).darrays[0].data)
- mtrx[i] = d
- except:
- print(sub[i] + " : NO " + filename)
- pass
-
-
- elif area == "hipp" and "_thickness" not in filename:
- for i in range(nSub):
- try:
- dpath = out + "/" + sub[i] + "/" + session + "/" + path
- d = []
- for _, h in enumerate(['lh', 'rh']):
- d = np.append(d, nib.load(dpath + sub[i] + '_' + session +
- filename.format(h)).darrays[0].data)
- mtrx[i] = d
- except:
- print(sub[i] + " : NO " + filename)
- pass
-
- # Compute asymmetry
- mtrx = (mtrx[:, 0:int(nvert / 2)] - mtrx[:, int(nvert / 2):nvert]) / (
- (mtrx[:, 0:int(nvert / 2)] + mtrx[:, int(nvert / 2):nvert]) / 2)
-
- # Create dataframe
- mtrx = pd.DataFrame(mtrx)
-
- # Compute z-score data
- zmtrx = zscore(mtrx, TBL)
- return zmtrx
-
-if __name__ == "__main__":
- # Load demographics
- TBL = pd.read_excel(demo, engine='openpyxl', skiprows=[0,2]).dropna(how="all").dropna(axis=1, how="all")
- #TBL = TBL.drop(TBL[TBL.incl == 0].index).reset_index()
-
- #=============================================================================
- # Load, concatenate, and save CORTICAL features (T2-FLAIR, qT1, AD, thickness)
- #=============================================================================
- mv_c = []
- mvdic_c = {}
- if "flair" in featList_ctx:
- t2z_c = matrix("ctx", '/anat/surfaces/flair/', '_space-conte69-32k_desc-{}_flair_10mm.mgh', TBL)
- mv_c.append(t2z_c.values)
- mvdic_c.update({'flair': t2z_c})
- if "qt1" in featList_ctx:
- qt1z_c = matrix("ctx", '/anat/surfaces/qt1/', '_space-conte69-32k_desc-{}_qt1_10mm.mgh', TBL)
- mv_c.append(qt1z_c.values)
- mvdic_c.update({'qt1': qt1z_c})
- if "adc" in featList_ctx:
- adcz_c = matrix("ctx", '/dwi/surfaces/', '_space-conte69-32k_desc-{}_model-DTI_map-ADC_10mm.mgh', TBL)
- mv_c.append(adcz_c.values)
- mvdic_c.update({'adc': adcz_c})
- if "thickness" in featList_ctx:
- ctz_c = matrix("ctx", '/anat/surfaces/morphology/', '_space-conte69-32k_desc-{}_thickness_10mm.mgh', TBL)
- ctz_c = ctz_c * -1
- mv_c.append(ctz_c.values)
- mvdic_c.update({'thickness': ctz_c})
-
- if mv_c:
- mv_c = np.nanmean(mv_c, axis=0)
-
- # Save unthrehsolded map, then filter out noise or non-significant findings and save
- m = load_mask(surface_name="conte69", join=True).astype(int)
- mv_c_orig = np.tile(mv_c, (1, 2)) * m
- mv_c_unthr = mv_c_orig
- np.savetxt(os.path.join(os.path.dirname(out), "analysis", "asymmetry", "allSubjects_unthr_ctx-z.csv"),
- mv_c_orig, delimiter=",")
-
- mv_c[(mv_c > -thr) & (mv_c < thr)] = 0
- mv_c = np.tile(mv_c, (1, 2)) * m
- np.savetxt(os.path.join(os.path.dirname(out), "analysis", "asymmetry", "allSubjects_ctx-z.csv"),
- mv_c, delimiter=",")
-
- # Save individual features (unthresholded and thresholded)
- for (_, mvName) in enumerate(mvdic_c):
- mv = mvdic_c[mvName]
- np.savetxt(os.path.join(os.path.dirname(out), "analysis", "asymmetry", "allSubjects_unthr_ctx-{}.csv".
- format(str(mvName))), mv, delimiter=",")
- mv[(mv > -thr) & (mv < thr)] = 0
- mv = np.tile(mv, (1, 2)) * m
- np.savetxt(os.path.join(os.path.dirname(out), "analysis", "asymmetry", "allSubjects_ctx-{}.csv".
- format(str(mvName))), mv, delimiter=",")
-
-
-
- #================================================================================
- # Load, concatenate, and save SUBCORTICAL features (T2-FLAIR, qT1, AD, thickness)
- #================================================================================
- mv_s = []
- mvdic_s = {}
- if "flair" in featList_sctx:
- t2z_s = matrix("sctx", '/anat/surfaces/flair/', '_space-flair_subcortical-intensities.csv', TBL)
- mv_s.append(t2z_s.values)
- mvdic_s.update({'flair': t2z_s})
- if "qt1" in featList_sctx:
- qt1z_s = matrix("sctx", '/anat/surfaces/qt1/', '_space-qt1_subcortical-intensities.csv', TBL)
- mv_s.append(qt1z_s.values)
- mvdic_s.update({'qt1': qt1z_s})
- if "adc" in featList_sctx:
- adcz_s = matrix("sctx", '/dwi/surfaces/', '_space-dwi_subcortical-ADC.csv', TBL)
- mv_s.append(adcz_s.values)
- mvdic_s.update({'adc': adcz_s})
- if "thickness" in featList_sctx:
- ctz_s = matrix("sctx", '/anat/surfaces/morphology/sctx_volume/', '_sctx_volume.csv', TBL)
- ctz_s = ctz_s * -1
- mv_s.append(ctz_s.values)
- mvdic_s.update({'thickness': ctz_s})
-
- if mv_s:
- mv_s = np.nanmean(mv_s, axis=0)
- mv_s_unthr = np.tile(mv_s, (1, 2))
-
- # Save unthrehsolded map, then filter out noise or non-significant findings and save
- np.savetxt(os.path.join(os.path.dirname(out), "analysis", "asymmetry", "allSubjects_unthr_sctx-z.csv"),
- mv_s, delimiter=",")
- mv_s[(mv_s > -thr) & (mv_s < thr)] = 0
- mv_s = np.tile(mv_s, (1, 2))
- np.savetxt(os.path.join(os.path.dirname(out), "analysis", "asymmetry", "allSubjects_sctx-z.csv"),
- mv_s, delimiter=",")
-
- # Save individual features (unthresholded and thresholded)
- for (_, mvName) in enumerate(mvdic_s):
- mv = mvdic_s[mvName]
- np.savetxt(os.path.join(os.path.dirname(out), "analysis", "asymmetry", "allSubjects_unthr_sctx-{}.csv".
- format(str(mvName))), mv, delimiter=",")
- mv[(mv > -thr) & (mv < thr)] = 0
- mv = np.tile(mv, (1, 2))
- np.savetxt(os.path.join(os.path.dirname(out), "analysis", "asymmetry", "allSubjects_sctx-{}.csv".
- format(str(mvName))), mv, delimiter=",")
-
-
-
- #=========================================================================================
- # Load, concatenate, and save HIPPOCAMPAL features (T2-FLAIR, qT1, AD, thickness)
- #=========================================================================================
- mv_h = []
- mvdic_h = {}
- if "flair" in featList_hipp:
- t2z_h = matrix("hipp", '/anat/surfaces/flair/', '_hemi-{}_space-flair_desc-flair_N4_den-0p5mm_label-hipp_midthickness_10mm.func.gii', TBL)
- mv_h.append(t2z_h.values)
- mvdic_h.update({'flair': t2z_h})
- if "qt1" in featList_hipp:
- qt1z_h = matrix("hipp", '/anat/surfaces/qt1/', '_hemi-{}_space-qt1_desc-qt1_den-0p5mm_label-hipp_midthickness_10mm.func.gii', TBL)
- mv_h.append(qt1z_h.values)
- mvdic_h.update({'qt1': qt1z_h})
- if "adc" in featList_hipp:
- adcz_h = matrix("hipp", '/dwi/surfaces/', '_hemi-{}_space-dwi_desc-dwi-ADC_den-0p5mm_label-hipp_midthickness_10mm.func.gii', TBL)
- mv_h.append(adcz_h.values)
- mvdic_h.update({'adc': adcz_h})
- if "thickness" in featList_hipp:
- ctz_h = matrix("hipp", '/surf/', '_hemi-{}_space-T1w_den-0p5mm_label-hipp_thickness.shape.gii', TBL)
- ctz_h = ctz_h * -1
- mv_h.append(ctz_h.values)
- mvdic_h.update({'thickness': ctz_h})
-
- if mv_h:
- mv_h = np.nanmean(mv_h, axis=0)
- mv_h_unthr = np.tile(mv_h, (1, 2))
-
- # Save unthrehsolded map, then filter out noise or non-significant findings and save
- np.savetxt(os.path.join(os.path.dirname(out), "analysis", "asymmetry", "allSubjects_unthr_hipp-z.csv"),
- mv_h, delimiter=",")
- mv_h[(mv_h > -thr) & (mv_h < thr)] = 0
- mv_h = np.tile(mv_h, (1, 2))
- np.savetxt(os.path.join(os.path.dirname(out), "analysis", "asymmetry", "allSubjects_hipp-z.csv"),
- mv_h, delimiter=",")
-
- # Save individual features (unthresholded and thresholded)
- for (_, mvName) in enumerate(mvdic_h):
- mv = mvdic_h[mvName]
- np.savetxt(os.path.join(os.path.dirname(out), "analysis", "asymmetry", "allSubjects_unthr_hipp-{}.csv".
- format(str(mvName))), mv, delimiter=",")
- mv[(mv > -thr) & (mv < thr)] = 0
- mv = np.tile(mv, (1, 2))
- np.savetxt(os.path.join(os.path.dirname(out), "analysis", "asymmetry", "allSubjects_hipp-{}.csv".
- format(str(mvName))), mv, delimiter=",")
-
-
-
- # ==================================================================
- # Plot and save multivariate cortical, subcortical, and hippocampal asymmetry
- # ==================================================================
- # Get row number of subject
- rn = (TBL.loc[TBL['ID'] == subject].index).tolist()
-
- # Save and plot multivariate z-score | cortical
- fname = os.path.join(os.path.dirname(out), "analysis", "asymmetry", subject, session, subject + "_ctx-mz-unthr_{}.gii")
- mv_c_half = mv_c_unthr[rn, :]
- mv_h_lh = new_gifti_image(mv_c_half[:len(mv_c_half)//2])
- mv_h_rh = new_gifti_image(mv_c_half[len(mv_c_half)//2:])
- nib.save(mv_h_lh, fname.format('lh'))
- nib.save(mv_h_rh, fname.format('rh'))
-
- fname = os.path.join(os.path.dirname(out), "analysis", "asymmetry", subject, session, subject + "_ctx-mz_{}.gii")
- mv_c_half = mv_c[rn, :]
- mv_h_lh = new_gifti_image(mv_c_half[:len(mv_c_half)//2])
- mv_h_rh = new_gifti_image(mv_c_half[len(mv_c_half)//2:])
- nib.save(mv_h_lh, fname.format('lh'))
- nib.save(mv_h_rh, fname.format('rh'))
-
- fname = os.path.join(os.path.dirname(out), "analysis", "asymmetry", subject, session, subject + "_ctx-mz.png")
- mv_c_half[:, :len(mv_c_half.flatten())//2] = np.nan
- plot_cortical(array_name=mv_c_half, surface_name='conte69', size=(800, 180), zoom=1.18, cmap='RdBu_r',
- color_bar=True, color_range=(-3, 3), screenshot=True, view=['lateral', 'medial', 'medial', 'lateral'],
- filename=fname)
-
- # Save and plot multivariate z-score | subcortical
- fname = os.path.join(os.path.dirname(out), "analysis", "asymmetry", subject, session, subject + "_sctx-mz-unthr_{}.gii")
- mv_s_half = mv_s_unthr[rn, :].flatten()
- mv_h_lh = new_gifti_image(mv_s_half[:len(mv_s_half)//2])
- mv_h_rh = new_gifti_image(mv_s_half[len(mv_s_half)//2:])
- nib.save(mv_h_lh, fname.format('lh'))
- nib.save(mv_h_rh, fname.format('rh'))
-
- fname = os.path.join(os.path.dirname(out), "analysis", "asymmetry", subject, session, subject + "_sctx-mz_{}.gii")
- mv_s_half = mv_s[rn, :].flatten()
- mv_h_lh = new_gifti_image(mv_s_half[:len(mv_s_half)//2])
- mv_h_rh = new_gifti_image(mv_s_half[len(mv_s_half)//2:])
- nib.save(mv_h_lh, fname.format('lh'))
- nib.save(mv_h_rh, fname.format('rh'))
-
- mv_s_half[:len(mv_s_half)//2] = np.nan
- fname = os.path.join(os.path.dirname(out), "analysis", "asymmetry", subject, session, subject + "_sctx-mz.png")
- plot_subcortical(array_name=mv_s_half, ventricles=False, size=(800, 180), zoom=1.18, cmap='RdBu_r',
- color_bar=True, color_range=(-3, 3), screenshot=True, view=['lateral', 'medial', 'medial', 'lateral'],
- filename=fname)
-
- # Save and plot multivariate z-score | hippocampal unthresholded
- fname = os.path.join(os.path.dirname(out), "analysis", "asymmetry", subject, session, subject + "_hipp-mz-unthr_{}.gii")
- mv_h_half = mv_h_unthr[rn, :].flatten()
- mv_h_lh = new_gifti_image(mv_h_half[:len(mv_h_half)//2])
- mv_h_rh = new_gifti_image(mv_h_half[len(mv_h_half)//2:])
- nib.save(mv_h_lh, fname.format('lh'))
- nib.save(mv_h_rh, fname.format('rh'))
-
- # Save and plot multivariate z-score | hippocampal thresholded
- fname = os.path.join(os.path.dirname(out), "analysis", "asymmetry", subject, session, subject + "_hipp-mz_{}.gii")
- mv_h_half = mv_h[rn, :].flatten()
- mv_h_lh = new_gifti_image(mv_h_half[:len(mv_h_half)//2])
- mv_h_rh = new_gifti_image(mv_h_half[len(mv_h_half)//2:])
- nib.save(mv_h_lh, fname.format('lh'))
- nib.save(mv_h_rh, fname.format('rh'))
-
- mv_h_half[:len(mv_h_half) // 2] = np.nan
- fname = os.path.join(os.path.dirname(out), "analysis", "asymmetry", subject, session, subject + "_hipp-mz.png")
- plot_hippocampal(array_name=mv_h_half, size=(800, 180), zoom=1.18, cmap='RdBu_r',
- color_bar=True, color_range=(-3, 3), screenshot=True,
- view=['lateral', 'medial', 'ventral', 'dorsal'], filename=fname)
-
-
-
- # ================================================================
- # Plot and save univariate cortical, subcortical, and hippocampal asymmetry
- # ================================================================
- # Plot univariate z-score | cortical
- for (_, mvName) in enumerate(mvdic_c):
- fname = os.path.join(os.path.dirname(out), "analysis", "asymmetry", subject, session, subject + "_ctx-{}-unthr_{}.gii")
- mv_tmp = np.loadtxt(os.path.join(os.path.dirname(out), "analysis", "asymmetry", "allSubjects_ctx-{}.csv".
- format(str(mvName))), delimiter=",")[rn, :]
-
- mv_h_lh = new_gifti_image(mv_tmp[:len(mv_tmp)//2])
- mv_h_rh = new_gifti_image(mv_tmp[len(mv_tmp)//2:])
- nib.save(mv_h_lh, fname.format(str(mvName), 'lh'))
- nib.save(mv_h_rh, fname.format(str(mvName), 'rh'))
-
- mv_tmp[:, :len(mv_c_half.flatten())//2] = np.nan
- plot_cortical(array_name=mv_tmp, surface_name='conte69', size=(800, 180), zoom=1.18, cmap='RdBu_r',
- color_bar=True, color_range=(-3, 3), screenshot=True, view=['lateral', 'medial', 'medial', 'lateral'],
- filename=os.path.join(os.path.dirname(out), "analysis", "asymmetry", subject, session, subject + "_ctx-{}.png".
- format(str(mvName))))
-
- # Plot univariate z-score | subcortical
- for (_, mvName) in enumerate(mvdic_s):
- fname = os.path.join(os.path.dirname(out), "analysis", "asymmetry", subject, session, subject + "_sctx-{}-unthr.gii")
- mv_tmp = np.loadtxt(os.path.join(os.path.dirname(out), "analysis", "asymmetry", "allSubjects_sctx-{}.csv".
- format(str(mvName))), delimiter=",")[rn, :].flatten()
-
- mv_h_lh = new_gifti_image(mv_tmp[:len(mv_tmp)//2])
- mv_h_rh = new_gifti_image(mv_tmp[len(mv_tmp)//2:])
- nib.save(mv_h_lh, fname.format(str(mvName), 'lh'))
- nib.save(mv_h_rh, fname.format(str(mvName), 'rh'))
-
- mv_tmp[:len(mv_tmp)//2] = np.nan
- plot_subcortical(array_name=mv_tmp, ventricles=False, size=(800, 180), zoom=1.18, cmap='RdBu_r',
- color_bar=True, color_range=(-3, 3), screenshot=True, view=['lateral', 'medial', 'medial', 'lateral'],
- filename=os.path.join(os.path.dirname(out), "analysis", "asymmetry", subject,
- session, subject + "_sctx-{}.png".format(str(mvName))))
-
- # Plot univariate z-score | hippocampal
- for (_, mvName) in enumerate(mvdic_h):
- fname = os.path.join(os.path.dirname(out), "analysis", "asymmetry", subject, session, subject + "_hipp-{}-unthr_{}.gii")
- mv_tmp = np.loadtxt(os.path.join(os.path.dirname(out), "analysis", "asymmetry", "allSubjects_hipp-{}.csv".
- format(str(mvName))), delimiter=",")[rn, :].flatten()
-
- mv_h_lh = new_gifti_image(mv_tmp[:len(mv_tmp)//2])
- mv_h_rh = new_gifti_image(mv_tmp[len(mv_tmp)//2:])
- nib.save(mv_h_lh, fname.format(str(mvName), 'lh'))
- nib.save(mv_h_rh, fname.format(str(mvName), 'rh'))
-
- mv_tmp[:len(mv_tmp) // 2] = np.nan
- plot_hippocampal(array_name=mv_tmp, size=(800, 180), zoom=1.18, cmap='RdBu_r',
- color_bar=True, color_range=(-3, 3), screenshot=True,
- view=['lateral', 'medial', 'ventral', 'dorsal'],
- filename=os.path.join(os.path.dirname(out), "analysis", "asymmetry", subject,
- session, subject + "_hipp-{}.png".format(str(mvName))))
-
\ No newline at end of file
diff --git a/lost+found/02_asymmetry.sh b/lost+found/02_asymmetry.sh
deleted file mode 100755
index a3f5364..0000000
--- a/lost+found/02_asymmetry.sh
+++ /dev/null
@@ -1,173 +0,0 @@
-#!/bin/bash
-#
-# Asymmetry analysis: Generates a multivariate featuer asymmety report
-#
-# This workflow makes use of many outputs and custom python scripts
-#
-# Atlas an templates are avaliable from:
-#
-# https://github.com/MICA-MNI/micapipe/tree/master/parcellations
-#
-# ARGUMENTS order:
-# $1 : BIDS directory
-# $2 : participant
-# $3 : Out Directory
-#
-umask 003
-BIDS=$1
-id=$2
-out=$3
-SES=$4
-nocleanup=$5
-threads=$6
-tmpDir=$7
-thr=$8
-demo=$9
-PROC=${10}
-export OMP_NUM_THREADS=$threads
-here=$(pwd)
-
-echo "PROC:"
-echo $PROC
-#------------------------------------------------------------------------------#
-# qsub configuration
-if [ "$PROC" = "qsub-MICA" ] || [ "$PROC" = "qsub-all.q" ];then
- export MICAPIPE=/data_/mica1/01_programs/micapipe
- source "${MICAPIPE}/functions/init.sh" "$threads"
-fi
-
-# source utilities
-source "$MICAPIPE/functions/utilities.sh"
-
-# Assigns variables names
-bids_variables "$BIDS" "$id" "$out" "$SES"
-
-#------------------------------------------------------------------------------#
-Title "asymmetry analysis\n\t\tmicapipe-z $Version, $PROC"
-micapipe_software
-bids_print.variables-post
-Info "wb_command will use $OMP_NUM_THREADS threads"
-Info "Saving temporal dir: $nocleanup"
-
-# Timer
-aloita=$(date +%s)
-Nsteps=0
-
-# Hippunfold directory
-hipDir="${out/micapipe/}/hippunfold_v1.0.0/hippunfold/"
-
-# Freesurfer SUBJECTs directory
-export SUBJECTS_DIR=${dir_surf}
-
-# Check CORTICAL input features (flair, qt1, thickness, adc); if missing a feature, skip this module
-featList_ctx=()
-if [ -f "${proc_struct}/surfaces/flair/${idBIDS}_space-conte69-32k_desc-rh_flair_10mm.mgh" ]; then
- featList_ctx+=("flair"); fi
-if [ -f "${proc_struct}/surfaces/qt1/${idBIDS}_space-conte69-32k_desc-rh_qt1_10mm.mgh" ]; then
- featList_ctx+=("qt1"); fi
-if [ -f "${proc_dwi}/surfaces/${idBIDS}_space-conte69-32k_desc-rh_model-DTI_map-ADC_10mm.mgh" ]; then
- featList_ctx+=("adc"); fi
-if [ -f "${proc_struct}/surfaces/morphology/${idBIDS}_space-conte69-32k_desc-rh_thickness_10mm.mgh" ]; then
- featList_ctx+=("thickness"); fi
-
-# Check SUBCORTICAL input features (flair, qt1, thickness, adc); if missing a feature, skip this module
-featList_sctx=()
-if [ -f "${proc_struct}/surfaces/flair/${idBIDS}_space-flair_subcortical-intensities.csv" ]; then
- featList_sctx+=("flair"); fi
-if [ -f "${proc_struct}/surfaces/qt1/${idBIDS}_space-qt1_subcortical-intensities.csv" ]; then
- featList_sctx+=("qt1"); fi
-if [ -f "${proc_dwi}/surfaces/${idBIDS}_space-dwi_subcortical-ADC.csv" ]; then
- featList_sctx+=("adc"); fi
-if [ -f "${proc_struct}/surfaces/morphology/sctx_volume/${idBIDS}_sctx_volume.csv" ]; then
- featList_sctx+=("thickness"); fi
-
-# Check HIPPOCAMPAL input features (flair, qt1, thickness, adc); if missing a feature, skip this module
-featList_hipp=()
-if [ -f "${proc_struct}/surfaces/flair/${idBIDS}_hemi-rh_space-flair_desc-flair_N4_den-0p5mm_label-hipp_midthickness_2mm.func.gii" ]; then
- featList_hipp+=("flair"); fi
-if [ -f "${proc_struct}/surfaces/qt1/${idBIDS}_hemi-rh_space-qt1_desc-qt1_den-0p5mm_label-hipp_midthickness_2mm.func.gii" ]; then
- featList_hipp+=("qt1"); fi
-if [ -f "${proc_dwi}/surfaces/${idBIDS}_hemi-rh_space-dwi_desc-dwi-ADC_den-0p5mm_label-hipp_midthickness_2mm.func.gii" ]; then
- featList_hipp+=("adc"); fi
-if [ -f "${hipDir}/sub-${id}/surf/sub-${id}_hemi-R_space-T1w_den-0p5mm_label-hipp_thickness.shape.gii" ]; then
- featList_hipp+=("thickness"); fi
-
-# Check input feature
-Info "Inputs:"
-if [[ "$thr" == DEFAULT ]]; then Note "default z-score threshold at |1.96|"; else
-Note "thr :" "$thr"; fi
-if [[ "$demo" == '' ]]; then Note "No demographic file specified"; else
-Note "demo :" "$demo"; fi
-
-# Create script specific temp directory
-tmp="${tmpDir}/${RANDOM}_micapipe-z_asymmetry_${idBIDS}"
-Do_cmd mkdir -p "$tmp"
-
-# TRAP in case the script fails
-trap 'cleanup $tmp $nocleanup $here' SIGINT SIGTERM
-
-# Make output directory
-outDir="${out/micapipe/}/analysis/asymmetry/sub-${id}/${SES}/"
-[[ ! -d "$outDir" ]] && Do_cmd mkdir -p "$outDir"
-
-
-#------------------------------------------------------------------------------#
-### asymmetry analysis ###
-Do_cmd python "$ZBRAINS"/functions/02_asymmetry.py "sub-$id" "$SES" "$out" "$hipDir" \
- "${demo}" "$thr" \
- --featList_ctx "${featList_ctx[@]}" --featList_sctx "${featList_sctx[@]}" \
- --featList_hipp "${featList_hipp[@]}"
-
-# Generate pdf report: asymmetry
-Do_cmd python "$ZBRAINS"/functions/02_asymmetry-report.py "sub-$id" "$SES" "$out" \
- "${demo}" "${thr}" \
- --featList_ctx "${featList_ctx[@]}" --featList_sctx "${featList_sctx[@]}" \
- --featList_hipp "${featList_hipp[@]}"
-
-if [[ -f "${outDir}/sub-${id}_asymmetryReport.pdf" ]]; then ((Nsteps++)); fi
-
-
-#------------------------------------------------------------------------------#
-# Surface register back to fsspace
-for hemi in lh rh; do
- [[ "$hemi" == lh ]] && hemisphere=l || hemisphere=r
- HEMICAP=$(echo $hemisphere | tr [:lower:] [:upper:])
-
- Do_cmd mri_convert ${outDir}/sub-${id}_ctx-mz_${hemi}.mgh ${outDir}/sub-${id}_ctx-mz_${hemi}.func.gii
-
- Do_cmd wb_command -metric-resample \
- ${outDir}/sub-${id}_ctx-mz_${hemi}.func.gii \
- "${util_surface}/fs_LR-deformed_to-fsaverage.${HEMICAP}.sphere.32k_fs_LR.surf.gii" \
- "${dir_conte69}/${idBIDS}_${hemi}_sphereReg.surf.gii" \
- ADAP_BARY_AREA \
- "${tmp}/${idBIDS}_space-fsspace-desc-${hemi}_ctx-mz.func.gii" \
- -area-surfs \
- "${dir_conte69}/${idBIDS}_space-conte69-32k_desc-${hemi}_midthickness.surf.gii" \
- "${dir_freesurfer}/surf/${hemi}.midthickness.surf.gii"
-
- Do_cmd mri_convert "${tmp}/${idBIDS}_space-fsspace-desc-${hemi}_ctx-mz.func.gii" \
- "${tmp}/${idBIDS}_space-fsspace-desc-${hemi}_ctx-mz.mgh"
-done
-
-# Map surface to volume
-Do_cmd mri_surf2vol --o "${outDir}/sub-${id}_ctx-mz.nii.gz" \
- --subject "$idBIDS" \
- --so ${dir_freesurfer}/surf/lh.white "${tmp}/${idBIDS}_space-fsspace-desc-lh_ctx-mz.mgh" \
- --so ${dir_freesurfer}/surf/rh.white "${tmp}/${idBIDS}_space-fsspace-desc-rh_ctx-mz.mgh"
-
-if [[ -f "${outDir}/sub-${id}_ctx-mz.nii.gz" ]]; then ((Nsteps++)); fi
-
-#------------------------------------------------------------------------------#
-# QC notification of completion
-lopuu=$(date +%s)
-eri=$(echo "$lopuu - $aloita" | bc)
-eri=$(echo print "$eri"/60 | perl)
-
-# Notification of completion
-if [ "$Nsteps" -eq 2 ]; then status="COMPLETED"; else status="ERROR asymmetry is missing a processing step"; fi
-Title "asymmetry processing ended in \033[38;5;220m $(printf "%0.3f\n" "$eri") minutes \033[38;5;141m.
-\tSteps completed : $(printf "%02d" "$Nsteps")/02
-\tStatus : ${status}
-\tCheck logs : $(ls "${dir_logs}"/asymmetry_*.txt)"
-echo "${id}, ${SES/ses-/}, asymmetry, $status N=$(printf "%02d" "$Nsteps")/02, $(whoami), $(uname -n), $(date), $(printf "%0.3f\n" "$eri"), ${PROC}, ${Version}" >> "${out}/micapipez_processed_sub.csv"
-cleanup "$tmp" "$nocleanup" "$here"
diff --git a/lost+found/02_regional_analysis.py b/lost+found/02_regional_analysis.py
deleted file mode 100644
index 611e701..0000000
--- a/lost+found/02_regional_analysis.py
+++ /dev/null
@@ -1,404 +0,0 @@
-import os
-import sys
-import pandas as pd
-import numpy as np
-import nibabel as nib
-import glob
-from enigmatoolbox.utils.useful import zscore_matrix
-from enigmatoolbox.datasets import load_mask
-from enigmatoolbox.plotting import plot_cortical, plot_subcortical, plot_hippocampal
-from enigmatoolbox.datasets import getaffine
-import warnings
-
-warnings.simplefilter('ignore')
-
-import argparse
-# defined command line options
-# this also generates --help and error handling
-CLI=argparse.ArgumentParser()
-CLI.add_argument("subject", nargs=1, type=str, default="")
-CLI.add_argument("session", nargs=1, type=str, default="")
-CLI.add_argument("out", nargs=1, type=str, default="")
-CLI.add_argument("hipDir", nargs=1, type=str, default="")
-CLI.add_argument("demo", nargs=1, type=str, default="")
-CLI.add_argument("thr", nargs=1, type=float, default=1.96)
-CLI.add_argument(
- "--featList_ctx", # name on the CLI - drop the `--` for positional/required parameters
- nargs="*",
- type=str,
- default=['flair', 'qt1', 'adc', 'thickness'])
-CLI.add_argument(
- "--featList_sctx",
- nargs="*",
- type=str, # any type/callable can be used here
- default=[])
-CLI.add_argument(
- "--featList_hipp",
- nargs="*",
- type=str, # any type/callable can be used here
- default=[])
-
-# parse the command line
-args = CLI.parse_args()
-
-# access CLI options
-subject = args.subject[0]
-session = args.session[0]
-out = args.out[0]
-hipDir = args.hipDir[0]
-demo = args.demo[0]
-thr = args.thr[0]
-featList_ctx = args.featList_ctx
-featList_sctx = args.featList_sctx
-featList_hipp = args.featList_hipp
-
-def zscore(mtrx, TBL):
- grpk = TBL['grp']
- data_z = zscore_matrix(mtrx, grpk, 'HC')
- return data_z
-
-def matrix(area, path, filename, TBL):
- # Get number of vertices/structures
- if area == "ctx":
- nvert = 64984
- elif area == "sctx":
- nvert = 14
- elif area == "hipp":
- nvert = 14524
-
- # Subject IDs
- sub = TBL['ID']
- (nSub,) = sub.shape
-
- # Initialize matrix
- mtrx = np.empty((nSub, nvert))
- mtrx[:] = np.nan
-
- # Cortical feature matrix
- if area == "ctx":
- for i in range(nSub):
- try:
- dpath = out + "/" + sub[i] + "/" + session + "/" + path
- d = []
- for _, h in enumerate(['lh', 'rh']):
- d = np.append(d, nib.load(dpath + sub[i] + '_' + session +
- filename.format(h)).get_fdata().squeeze())
- mtrx[i, :] = d
- except:
- print(sub[i] + " : NO " + filename)
- pass
-
- # Subcortical feature matrix
- elif area == "sctx":
- for i in range(nSub):
- try:
- dpath = out + "/" + sub[i] + "/" + session + "/" + path
- mtrx[i] = np.loadtxt(dpath + sub[i] + "_" + session + filename,
- delimiter=",", skiprows=1, usecols=range(1,15))
- except:
- print(sub[i] + " : NO " + filename)
- pass
-
- # Hippocampal feature matrix
- elif area == "hipp" and "_thickness" in filename:
- for i in range(nSub):
- try:
- dpath = hipDir + "/" + sub[i] + "/" + path
- d = []
- for _, h in enumerate(['L', 'R']):
- d = np.append(d, nib.load(dpath + sub[i] +
- filename.format(h)).darrays[0].data)
- mtrx[i] = d
- except:
- print(sub[i] + " : NO " + filename)
- pass
-
-
- elif area == "hipp" and "_thickness" not in filename:
- for i in range(nSub):
- try:
- dpath = out + "/" + sub[i] + "/" + session + "/" + path
- d = []
- for _, h in enumerate(['lh', 'rh']):
- d = np.append(d, nib.load(dpath + sub[i] + '_' + session +
- filename.format(h)).darrays[0].data)
- mtrx[i] = d
- except:
- print(sub[i] + " : NO " + filename)
- pass
-
- # Get vertexwise z-scores
- mtrx_full = mtrx
- mtrx_full = pd.DataFrame(mtrx_full)
- zmtrx_full = zscore(mtrx_full, TBL)
-
- return zmtrx_full
-
-if __name__ == "__main__":
- # Load demographics
- TBL = pd.read_excel(demo, engine='openpyxl', skiprows=[0,2]).dropna(how="all").dropna(axis=1, how="all")
- #TBL = TBL.drop(TBL[TBL.incl == 0].index).reset_index()
-
- #=============================================================================
- # Load, concatenate, and save CORTICAL features (T2-FLAIR, qT1, AD, thickness)
- #=============================================================================
- mvFull_c = []
- mvFullDic_c = {}
- if "flair" in featList_ctx:
- t2zFull_c = matrix("ctx", '/anat/surfaces/flair/', '_space-conte69-32k_desc-{}_flair_10mm.mgh', TBL)
- mvFull_c.append(t2zFull_c.values)
- mvFullDic_c.update({'flair': t2zFull_c})
- if "qt1" in featList_ctx:
- qt1zFull_c = matrix("ctx", '/anat/surfaces/qt1/', '_space-conte69-32k_desc-{}_qt1_10mm.mgh', TBL)
- mvFull_c.append(qt1zFull_c.values)
- mvFullDic_c.update({'qt1': qt1zFull_c})
- if "adc" in featList_ctx:
- adczFull_c = matrix("ctx", '/dwi/surfaces/', '_space-conte69-32k_desc-{}_model-DTI_map-ADC_10mm.mgh', TBL)
- mvFull_c.append(adczFull_c.values)
- mvFullDic_c.update({'adc': adczFull_c})
- if "thickness" in featList_ctx:
- ctzFull_c = matrix("ctx", '/anat/surfaces/morphology/', '_space-conte69-32k_desc-{}_thickness_10mm.mgh', TBL)
- ctzFull_c = ctzFull_c * -1
- mvFull_c.append(ctzFull_c.values)
- mvFullDic_c.update({'thickness': ctzFull_c})
-
-
- if mvFull_c:
- mvFull_c = np.nanmean(mvFull_c, axis=0)
-
- # Save unthrehsolded map, then filter out noise or non-significant findings and save
- m = load_mask(surface_name="conte69", join=True).astype(int)
- mv_c_orig = mvFull_c * m
- mvFull_c_unthr = mv_c_orig
- np.savetxt(os.path.join(os.path.dirname(out), "analysis", "regional", "allSubjects_unthr_ctx-z.csv"),
- mv_c_orig, delimiter=",")
- mvFull_c[(mvFull_c > -thr) & (mvFull_c < thr)] = 0
- np.savetxt(os.path.join(os.path.dirname(out), "analysis", "regional", "allSubjects_ctx-z.csv"),
- mvFull_c, delimiter=",")
-
- # Save individual features (unthresholded and thresholded)
- for (_, mvName) in enumerate(mvFullDic_c):
- mv = mvFullDic_c[mvName]
- np.savetxt(os.path.join(os.path.dirname(out), "analysis", "regional", "allSubjects_unthr_ctx-{}.csv".
- format(str(mvName))), mv, delimiter=",")
- mv[(mv > -thr) & (mv < thr)] = 0
- np.savetxt(os.path.join(os.path.dirname(out), "analysis", "regional", "allSubjects_ctx-{}.csv".
- format(str(mvName))), mv, delimiter=",")
-
-
-
-
- #================================================================================
- # Load, concatenate, and save SUBCORTICAL features (T2-FLAIR, qT1, AD, thickness)
- #================================================================================
- mvFull_s = []
- mvFullDic_s = {}
- if "flair" in featList_sctx:
- t2zFull_s = matrix("sctx", '/anat/surfaces/flair/', '_space-flair_subcortical-intensities.csv', TBL)
- mvFull_s.append(t2zFull_s.values)
- mvFullDic_s.update({'flair': t2zFull_s})
- if "qt1" in featList_sctx:
- qt1zFull_s = matrix("sctx", '/anat/surfaces/qt1/', '_space-qt1_subcortical-intensities.csv', TBL)
- mvFull_s.append(qt1zFull_s.values)
- mvFullDic_s.update({'qt1': qt1zFull_s})
- if "adc" in featList_sctx:
- adczFull_s = matrix("sctx", '/dwi/surfaces/', '_space-dwi_subcortical-ADC.csv', TBL)
- mvFull_s.append(adczFull_s.values)
- mvFullDic_s.update({'adc': adczFull_s})
- if "thickness" in featList_sctx:
- ctzFull_s = matrix("sctx", '/anat/surfaces/morphology/sctx_volume/', '_sctx_volume.csv', TBL)
- ctzFull_s = ctzFull_s * -1
- mvFull_s.append(ctzFull_s.values)
- mvFullDic_s.update({'thickness': ctzFull_s})
-
-
- if mvFull_s:
- mvFull_s = np.nanmean(mvFull_s, axis=0)
- mvFull_s_unthr = mvFull_s
-
- # Save unthrehsolded map, then filter out noise or non-significant findings and save
- np.savetxt(os.path.join(os.path.dirname(out), "analysis", "regional", "allSubjects_unthr_sctx-z.csv"),
- mvFull_s_unthr, delimiter=",")
- mvFull_s[(mvFull_s > -thr) & (mvFull_s < thr)] = 0
- np.savetxt(os.path.join(os.path.dirname(out), "analysis", "regional", "allSubjects_sctx-z.csv"),
- mvFull_s, delimiter=",")
-
- # Save individual features (unthresholded and thresholded)
- for (_, mvName) in enumerate(mvFullDic_s):
- mv = mvFullDic_s[mvName]
- np.savetxt(os.path.join(os.path.dirname(out), "analysis", "regional", "allSubjects_unthr_sctx-{}.csv".
- format(str(mvName))), mv, delimiter=",")
- mv[(mv > -thr) & (mv < thr)] = 0
- np.savetxt(os.path.join(os.path.dirname(out), "analysis", "regional", "allSubjects_sctx-{}.csv".
- format(str(mvName))), mv, delimiter=",")
-
-
-
- #=========================================================================================
- # Load, concatenate, and save HIPPOCAMPAL features (T2-FLAIR, qT1, AD, thickness)
- #=========================================================================================
- mvFull_h = []
- mvFullDic_h = {}
- if "flair" in featList_hipp:
- t2zFull_h = matrix("hipp", '/anat/surfaces/flair/', '_hemi-{}_space-flair_desc-flair_N4_den-0p5mm_label-hipp_midthickness_10mm.func.gii', TBL)
- mvFull_h.append(t2zFull_h.values)
- mvFullDic_h.update({'flair': t2zFull_h})
- if "qt1" in featList_hipp:
- qt1zFull_h = matrix("hipp", '/anat/surfaces/qt1/', '_hemi-{}_space-qt1_desc-qt1_den-0p5mm_label-hipp_midthickness_10mm.func.gii', TBL)
- mvFull_h.append(qt1zFull_h.values)
- mvFullDic_h.update({'qt1': qt1zFull_h})
- if "adc" in featList_hipp:
- adczFull_h = matrix("hipp", '/dwi/surfaces/', '_hemi-{}_space-dwi_desc-dwi-ADC_den-0p5mm_label-hipp_midthickness_10mm.func.gii', TBL)
- mvFull_h.append(adczFull_h.values)
- mvFullDic_h.update({'adc': adczFull_h})
- if "thickness" in featList_hipp:
- ctzFull_h = matrix("hipp", '/surf/', '_hemi-{}_space-T1w_den-0p5mm_label-hipp_thickness.shape.gii', TBL)
- ctzFull_h = ctzFull_h * -1
- mvFull_h.append(ctzFull_h.values)
- mvFullDic_h.update({'thickness': ctzFull_h})
-
- if mvFull_h:
- mvFull_h = np.nanmean(mvFull_h, axis=0)
- mvFull_h_unthr = mvFull_h
-
- # Save unthrehsolded map, then filter out noise or non-significant findings and save
- np.savetxt(os.path.join(os.path.dirname(out), "analysis", "regional", "allSubjects_unthr_hipp-z.csv"),
- mvFull_h_unthr, delimiter=",")
- mvFull_h[(mvFull_h > -thr) & (mvFull_h < thr)] = 0
- np.savetxt(os.path.join(os.path.dirname(out), "analysis", "regional", "allSubjects_hipp-z.csv"),
- mvFull_h, delimiter=",")
-
- # Save individual features (unthresholded and thresholded)
- for (_, mvName) in enumerate(mvFullDic_h):
- mv = mvFullDic_h[mvName]
- np.savetxt(os.path.join(os.path.dirname(out), "analysis", "regional", "allSubjects_unthr_hipp-{}.csv".
- format(str(mvName))), mv, delimiter=",")
- mv[(mv > -thr) & (mv < thr)] = 0
- np.savetxt(os.path.join(os.path.dirname(out), "analysis", "regional", "allSubjects_hipp-{}.csv".
- format(str(mvName))), mv, delimiter=",")
-
-
-
- # ========================================================================
- # Plot multivariate cortical, subcortical, and hippocampal features (full)
- # ========================================================================
- # Get row number of subject
- rn = (TBL.loc[TBL['ID'] == subject].index).tolist()
-
- # Save and plot multivariate z-score | cortical
- mv_c_full = mvFull_c_unthr[rn, :]
- fname = os.path.join(os.path.dirname(out), "analysis", "regional", subject, session, subject + "_ctx-mz-unthr_{}.mgh")
- nib.freesurfer.mghformat.MGHImage(np.float32(mv_c_full[:, :mv_c_full.shape[1]//2].flatten()),
- getaffine('conte69', 'lh')).to_filename(fname.format('lh'))
- nib.freesurfer.mghformat.MGHImage(np.float32(mv_c_full[:, mv_c_full.shape[1]//2:].flatten()),
- getaffine('conte69', 'rh')).to_filename(fname.format('rh'))
-
- mv_c_full = mvFull_c[rn, :]
- fname = os.path.join(os.path.dirname(out), "analysis", "regional", subject, session, subject + "_ctx-mz_{}.mgh")
- nib.freesurfer.mghformat.MGHImage(np.float32(mv_c_full[:, :mv_c_full.shape[1]//2].flatten()),
- getaffine('conte69', 'lh')).to_filename(fname.format('lh'))
- nib.freesurfer.mghformat.MGHImage(np.float32(mv_c_full[:, mv_c_full.shape[1]//2:].flatten()),
- getaffine('conte69', 'rh')).to_filename(fname.format('rh'))
-
- fname = os.path.join(os.path.dirname(out), "analysis", "regional", subject, session, subject + "_ctx-mz.png")
- plot_cortical(array_name=mv_c_full, surface_name='conte69', size=(800, 180), zoom=1.18, cmap='RdBu_r',
- color_bar=True, color_range=(-3, 3), screenshot=True, view=['lateral', 'medial', 'medial', 'lateral'],
- filename=fname)
-
- # Save and plot multivariate z-score | subcortical
- mv_tmp = np.loadtxt(os.path.join(os.path.dirname(out), "analysis", "regional", "allSubjects_unthr_sctx-z.csv"),
- delimiter=",")[rn, :].flatten()
- np.savetxt(os.path.join(os.path.dirname(out), "analysis", "regional", subject, session, subject + "_sctx-mz-unthr.txt"),
- mv_tmp, delimiter=",")
-
- mv_s_full = mvFull_s[rn, :].flatten()
- np.savetxt(os.path.join(os.path.dirname(out), "analysis", "regional", subject, session, subject + "_sctx-mz.txt"),
- mv_s_full, delimiter=",")
-
- fname = os.path.join(os.path.dirname(out), "analysis", "regional", subject, session, subject + "_sctx-mz.png")
- plot_subcortical(array_name=mv_s_full, ventricles=False, size=(800, 180), zoom=1.18, cmap='RdBu_r',
- color_bar=True, color_range=(-3, 3), screenshot=True, view=['lateral', 'medial', 'medial', 'lateral'],
- filename=fname)
-
- # Save and plot multivariate z-score | hippocampal
- mv_h_full = np.loadtxt(os.path.join(os.path.dirname(out), "analysis", "regional", "allSubjects_unthr_hipp-z.csv"),
- delimiter=",")[rn, :].flatten()
- mv_h_lh = nib.gifti.gifti.GiftiImage()
- mv_h_lh.add_gifti_data_array(nib.gifti.gifti.GiftiDataArray(data=mv_h_full[:len(mv_h_full)//2]))
- mv_h_rh = nib.gifti.gifti.GiftiImage()
- mv_h_rh.add_gifti_data_array(nib.gifti.gifti.GiftiDataArray(data=mv_h_full[len(mv_h_full)//2:]))
- fname = os.path.join(os.path.dirname(out), "analysis", "regional", subject, session, subject + "_hipp-mz-unthr_{}.gii")
- nib.save(mv_h_lh, fname.format('lh'))
- nib.save(mv_h_rh, fname.format('rh'))
-
- mv_h_full = mvFull_h[rn, :].flatten()
- mv_h_lh = nib.gifti.gifti.GiftiImage()
- mv_h_lh.add_gifti_data_array(nib.gifti.gifti.GiftiDataArray(data=mv_h_full[:len(mv_h_full)//2]))
- mv_h_rh = nib.gifti.gifti.GiftiImage()
- mv_h_rh.add_gifti_data_array(nib.gifti.gifti.GiftiDataArray(data=mv_h_full[len(mv_h_full)//2:]))
- fname = os.path.join(os.path.dirname(out), "analysis", "regional", subject, session, subject + "_hipp-mz_{}.gii")
- nib.save(mv_h_lh, fname.format('lh'))
- nib.save(mv_h_rh, fname.format('rh'))
-
- fname = os.path.join(os.path.dirname(out), "analysis", "regional", subject, session, subject + "_hipp-mz.png")
- plot_hippocampal(array_name=mv_h_full, size=(800, 180), zoom=1.18, cmap='RdBu_r',
- color_bar=True, color_range=(-3, 3), screenshot=True,
- view=['lateral', 'ventral', 'dorsal', 'medial'], filename=fname)
-
-
-
- # ======================================================================
- # Plot univariate cortical, subcortical, and hippocampal features (full)
- # ======================================================================
- # Plot univariate z-score | cortical
- for (_, mvName) in enumerate(mvFullDic_c):
- mv_tmp = np.loadtxt(os.path.join(os.path.dirname(out), "analysis", "regional", "allSubjects_ctx-{}.csv".
- format(str(mvName))), delimiter=",")[rn, :]
- plot_cortical(array_name=mv_tmp, surface_name='conte69', size=(800, 180), zoom=1.18, cmap='RdBu_r',
- color_bar=True, color_range=(-3, 3), screenshot=True, view=['lateral', 'medial', 'medial', 'lateral'],
- filename=os.path.join(os.path.dirname(out), "analysis", "regional", subject, session, subject + "_ctx-{}.png".
- format(str(mvName))))
-
- mv_tmp = np.loadtxt(os.path.join(os.path.dirname(out), "analysis", "regional", "allSubjects_unthr_ctx-{}.csv".
- format(str(mvName))), delimiter=",")[rn, :]
- fname = os.path.join(os.path.dirname(out), "analysis", "regional", subject, session, subject + "_ctx-{}-unthr_{}.mgh")
- nib.freesurfer.mghformat.MGHImage(np.float32(mv_tmp[:, :mv_tmp.shape[1]//2].flatten()),
- getaffine('conte69', 'lh')).to_filename(fname.format(str(mvName), 'lh'))
- nib.freesurfer.mghformat.MGHImage(np.float32(mv_tmp[:, mv_tmp.shape[1]//2:].flatten()),
- getaffine('conte69', 'rh')).to_filename(fname.format(str(mvName), 'rh'))
-
-
- # Plot univariate z-score | subcortical
- for (_, mvName) in enumerate(mvFullDic_s):
- mv_tmp = np.loadtxt(os.path.join(os.path.dirname(out), "analysis", "regional", "allSubjects_sctx-{}.csv".
- format(str(mvName))), delimiter=",")[rn, :].flatten()
- plot_subcortical(array_name=mv_tmp, ventricles=False, size=(800, 180), zoom=1.18, cmap='RdBu_r',
- color_bar=True, color_range=(-3, 3), screenshot=True, view=['lateral', 'medial', 'medial', 'lateral'],
- filename=os.path.join(os.path.dirname(out), "analysis", "regional", subject,
- session, subject + "_sctx-{}.png".format(str(mvName))))
-
- mv_tmp = np.loadtxt(os.path.join(os.path.dirname(out), "analysis", "regional", "allSubjects_unthr_sctx-{}.csv").
- format(str(mvName)), delimiter=",")[rn, :].flatten()
- np.savetxt(os.path.join(os.path.dirname(out), "analysis", "regional", subject, session, subject + "_sctx-{}-unthr.txt").
- format(str(mvName)), mv_tmp, delimiter=",")
-
-
- # Plot univariate z-score | hippocampal
- for (_, mvName) in enumerate(mvFullDic_h):
- mv_tmp = np.loadtxt(os.path.join(os.path.dirname(out), "analysis", "regional", "allSubjects_hipp-{}.csv".
- format(str(mvName))), delimiter=",")[rn, :].flatten()
- plot_hippocampal(array_name=mv_tmp, size=(800, 180), zoom=1.18, cmap='RdBu_r',
- color_bar=True, color_range=(-3, 3), screenshot=True,
- view=['lateral', 'ventral', 'dorsal', 'medial'],
- filename=os.path.join(os.path.dirname(out), "analysis", "regional", subject,
- session, subject + "_hipp-{}.png".format(str(mvName))))
-
- mv_tmp = np.loadtxt(os.path.join(os.path.dirname(out), "analysis", "regional", "allSubjects_unthr_hipp-{}.csv".
- format(str(mvName))), delimiter=",")[rn, :].flatten()
- mv_h_lh = nib.gifti.gifti.GiftiImage()
- mv_h_lh.add_gifti_data_array(nib.gifti.gifti.GiftiDataArray(data=mv_tmp[:len(mv_tmp)//2]))
- mv_h_rh = nib.gifti.gifti.GiftiImage()
- mv_h_rh.add_gifti_data_array(nib.gifti.gifti.GiftiDataArray(data=mv_tmp[len(mv_tmp)//2:]))
- fname = os.path.join(os.path.dirname(out), "analysis", "regional", subject, session, subject + "_hipp-{}-unthr_{}.gii")
- nib.save(mv_h_lh, fname.format(str(mvName), 'lh'))
- nib.save(mv_h_rh, fname.format(str(mvName), 'rh'))
diff --git a/lost+found/02_regional_analysis.sh b/lost+found/02_regional_analysis.sh
deleted file mode 100755
index a63e12a..0000000
--- a/lost+found/02_regional_analysis.sh
+++ /dev/null
@@ -1,176 +0,0 @@
-#!/bin/bash
-#
-# Regional analysis: Generates a multivariate feature regional changes report
-#
-# This workflow makes use of many outputs and custom python scripts
-#
-# Atlases and templates are avaliable from:
-#
-# https://github.com/MICA-MNI/micapipe/tree/master/parcellations
-#
-# ARGUMENTS order:
-# $1 : BIDS directory
-# $2 : participant
-# $3 : Out Directory
-#
-umask 003
-BIDS=$1
-id=$2
-out=$3
-SES=$4
-nocleanup=$5
-threads=$6
-tmpDir=$7
-thr=$8
-demo=$9
-PROC=${10}
-export OMP_NUM_THREADS=$threads
-here=$(pwd)
-
-echo "PROC:"
-echo $PROC
-#------------------------------------------------------------------------------#
-# qsub configuration
-if [ "$PROC" = "qsub-MICA" ] || [ "$PROC" = "qsub-all.q" ];then
- export MICAPIPE=/data_/mica1/01_programs/micapipe
- source "${MICAPIPE}/functions/init.sh" "$threads"
-fi
-
-# source utilities
-source "$MICAPIPE/functions/utilities.sh"
-
-# Assigns variables names
-bids_variables "$BIDS" "$id" "$out" "$SES"
-
-#------------------------------------------------------------------------------#
-Title "regional analysis\n\t\tmicapipe-z $Version, $PROC"
-micapipe_software
-bids_print.variables-post
-Info "wb_command will use $OMP_NUM_THREADS threads"
-Info "Saving temporal dir: $nocleanup"
-
-# Timer
-aloita=$(date +%s)
-Nsteps=0
-
-# Hippunfold directory
-hipDir="${out/micapipe/}/hippunfold_v1.0.0/hippunfold/"
-
-# Freesurfer SUBJECTs directory
-export SUBJECTS_DIR=${dir_surf}
-
-# Check CORTICAL input features (flair, qt1, thickness, adc); if missing a feature, skip this module
-featList_ctx=()
-if [ -f "${proc_struct}/surfaces/flair/${idBIDS}_space-conte69-32k_desc-rh_flair_10mm.mgh" ]; then
- featList_ctx+=("flair"); fi
-if [ -f "${proc_struct}/surfaces/qt1/${idBIDS}_space-conte69-32k_desc-rh_qt1_10mm.mgh" ]; then
- featList_ctx+=("qt1"); fi
-if [ -f "${proc_dwi}/surfaces/${idBIDS}_space-conte69-32k_desc-rh_model-DTI_map-ADC_10mm.mgh" ]; then
- featList_ctx+=("adc"); fi
-if [ -f "${proc_struct}/surfaces/morphology/${idBIDS}_space-conte69-32k_desc-rh_thickness_10mm.mgh" ]; then
- featList_ctx+=("thickness"); fi
-
-# Check SUBCORTICAL input features (flair, qt1, thickness, adc); if missing a feature, skip this module
-featList_sctx=()
-if [ -f "${proc_struct}/surfaces/flair/${idBIDS}_space-flair_subcortical-intensities.csv" ]; then
- featList_sctx+=("flair"); fi
-if [ -f "${proc_struct}/surfaces/qt1/${idBIDS}_space-qt1_subcortical-intensities.csv" ]; then
- featList_sctx+=("qt1"); fi
-if [ -f "${proc_dwi}/surfaces/${idBIDS}_space-dwi_subcortical-ADC.csv" ]; then
- featList_sctx+=("adc"); fi
-if [ -f "${proc_struct}/surfaces/morphology/sctx_volume/${idBIDS}_sctx_volume.csv" ]; then
- featList_sctx+=("thickness"); fi
-
-# Check HIPPOCAMPAL input features (flair, qt1, thickness, adc); if missing a feature, skip this module
-featList_hipp=()
-if [ -f "${proc_struct}/surfaces/flair/${idBIDS}_hemi-rh_space-flair_desc-flair_N4_den-0p5mm_label-hipp_midthickness_10mm.func.gii" ]; then
- featList_hipp+=("flair"); fi
-if [ -f "${proc_struct}/surfaces/qt1/${idBIDS}_hemi-rh_space-qt1_desc-qt1_den-0p5mm_label-hipp_midthickness_10mm.func.gii" ]; then
- featList_hipp+=("qt1"); fi
-if [ -f "${proc_dwi}/surfaces/${idBIDS}_hemi-rh_space-dwi_desc-dwi-ADC_den-0p5mm_label-hipp_midthickness_10mm.func.gii" ]; then
- featList_hipp+=("adc"); fi
-if [ -f "${hipDir}/sub-${id}/surf/sub-${id}_hemi-R_space-T1w_den-0p5mm_label-hipp_thickness.shape.gii" ]; then
- featList_hipp+=("thickness"); fi
-
-# Check input feature
-Info "Inputs:"
-if [[ "$thr" == DEFAULT ]]; then Note "default z-score threshold at |1.96|"; else
-Note "thr :" "$thr"; fi
-if [[ "$demo" == '' ]]; then Note "No demographic file specified"; else
-Note "demo :" "$demo"; fi
-
-# Create script specific temp directory
-tmp="${tmpDir}/${RANDOM}_micapipe-z_asymmetry_${idBIDS}"
-Do_cmd mkdir -p "$tmp"
-
-# TRAP in case the script fails
-trap 'cleanup $tmp $nocleanup $here' SIGINT SIGTERM
-
-# Make output directory
-outDir="${out/micapipe/}/analysis/regional/sub-${id}/${SES}/"
-[[ ! -d "$outDir" ]] && Do_cmd mkdir -p "$outDir"
-
-
-#------------------------------------------------------------------------------#
-### regional analysis ###
-Do_cmd python "$ZBRAINS"/functions/02_regional_analysis.py "sub-$id" "$SES" "$out" "$hipDir" \
- "${demo}" "$thr" \
- --featList_ctx "${featList_ctx[@]}" --featList_sctx "${featList_sctx[@]}" \
- --featList_hipp "${featList_hipp[@]}"
-
-# Generate pdf report
-Do_cmd python "$ZBRAINS"/functions/02_regional_analysis_report.py "sub-$id" "$SES" "$out" \
- "${demo}" "${thr}" \
- --featList_ctx "${featList_ctx[@]}" --featList_sctx "${featList_sctx[@]}" \
- --featList_hipp "${featList_hipp[@]}"
-
-if [[ -f "${outDir}/sub-${id}_regional_analysis.pdf" ]]; then ((Nsteps++)); fi
-
-
-################################################
-# continue here
-################################################
-
-# Surface register back to fsspace
-for hemi in lh rh; do
- [[ "$hemi" == lh ]] && hemisphere=l || hemisphere=r
- HEMICAP=$(echo $hemisphere | tr [:lower:] [:upper:])
-
- Do_cmd mri_convert ${outDir}/sub-${id}_ctx-mz_${hemi}.mgh ${outDir}/sub-${id}_ctx-mz_${hemi}.func.gii
-
- Do_cmd wb_command -metric-resample \
- ${outDir}/sub-${id}_ctx-mz_${hemi}.func.gii \
- "${util_surface}/fs_LR-deformed_to-fsaverage.${HEMICAP}.sphere.32k_fs_LR.surf.gii" \
- "${dir_conte69}/${idBIDS}_${hemi}_sphereReg.surf.gii" \
- ADAP_BARY_AREA \
- "${tmp}/${idBIDS}_space-fsspace-desc-${hemi}_ctx-mz.func.gii" \
- -area-surfs \
- "${dir_conte69}/${idBIDS}_space-conte69-32k_desc-${hemi}_midthickness.surf.gii" \
- "${dir_freesurfer}/surf/${hemi}.midthickness.surf.gii"
-
- Do_cmd mri_convert "${tmp}/${idBIDS}_space-fsspace-desc-${hemi}_ctx-mz.func.gii" \
- "${tmp}/${idBIDS}_space-fsspace-desc-${hemi}_ctx-mz.mgh"
-done
-
-# Map surface to volume
-Do_cmd mri_surf2vol --o "${outDir}/sub-${id}_ctx-mz.nii.gz" \
- --subject "$idBIDS" \
- --so ${dir_freesurfer}/surf/lh.white "${tmp}/${idBIDS}_space-fsspace-desc-lh_ctx-mz.mgh" \
- --so ${dir_freesurfer}/surf/rh.white "${tmp}/${idBIDS}_space-fsspace-desc-rh_ctx-mz.mgh"
-
-if [[ -f "${outDir}/sub-${id}_ctx-mz.nii.gz" ]]; then ((Nsteps++)); fi
-
-#------------------------------------------------------------------------------#
-# QC notification of completion
-lopuu=$(date +%s)
-eri=$(echo "$lopuu - $aloita" | bc)
-eri=$(echo print "$eri"/60 | perl)
-
-# Notification of completion
-if [ "$Nsteps" -eq 2 ]; then status="COMPLETED"; else status="ERROR regional analysis is missing a processing step"; fi
-Title "asymmetry processing ended in \033[38;5;220m $(printf "%0.3f\n" "$eri") minutes \033[38;5;141m.
-\tSteps completed : $(printf "%02d" "$Nsteps")/02
-\tStatus : ${status}
-\tCheck logs : $(ls "${dir_logs}"/regional_z_*.txt)"
-echo "${id}, ${SES/ses-/}, regional_z, $status N=$(printf "%02d" "$Nsteps")/02, $(whoami), $(uname -n), $(date), $(printf "%0.3f\n" "$eri"), ${PROC}, ${Version}" >> "${out}/micapipez_processed_sub.csv"
-cleanup "$tmp" "$nocleanup" "$here"
diff --git a/lost+found/02_regional_analysis_report.py b/lost+found/02_regional_analysis_report.py
deleted file mode 100644
index 015b39b..0000000
--- a/lost+found/02_regional_analysis_report.py
+++ /dev/null
@@ -1,272 +0,0 @@
-from xhtml2pdf import pisa
-import sys
-import pandas as pd
-import os
-import argparse
-
-# defined command line options
-CLI=argparse.ArgumentParser()
-CLI.add_argument("subject", nargs=1, type=str, default="")
-CLI.add_argument("session", nargs=1, type=str, default="")
-CLI.add_argument("out", nargs=1, type=str, default="")
-CLI.add_argument("demo", nargs=1, type=str, default="")
-CLI.add_argument("thr", nargs=1, type=float, default=1.96)
-CLI.add_argument(
- "--featList_ctx",
- nargs="*",
- type=str,
- default=['flair', 'qt1', 'adc', 'thickness'])
-CLI.add_argument(
- "--featList_sctx",
- nargs="*",
- type=str,
- default=[])
-CLI.add_argument(
- "--featList_hipp",
- nargs="*",
- type=str,
- default=[])
-
-# parse the command line
-args = CLI.parse_args()
-
-# access CLI options
-subject = args.subject[0]
-session = args.session[0]
-out = args.out[0]
-demo = args.demo[0]
-thr = args.thr[0]
-featList_ctx = ', '.join(args.featList_ctx)
-featList_sctx = ', '.join(args.featList_sctx)
-featList_hipp = ', '.join(args.featList_hipp)
-
-TBL = pd.read_excel(demo, engine='openpyxl', skiprows=[0,2]).dropna(how="all").dropna(axis=1, how="all")
-
-# Get participant's sex
-if TBL[TBL['ID'] == subject]['sex'].to_list()[0] == 'F':
- sex = 'Female'
-elif TBL[TBL['ID'] == subject]['sex'].to_list()[0] == 'M':
- sex = 'Male'
-else:
- sex = 'Not defined'
-
-# Get participant's age
-try:
- age = TBL[TBL['ID'] == subject]['age'].to_list()[0]
-except:
- age = 'Not defined'
-
-# Path to figures
-figPath = os.path.join(os.path.dirname(out), "analysis", "regional", subject, session) + "/" + subject
-
-def report_block_template(figPath, subject='', session='', sex='', age='', featList_ctx='', featList_sctx='', featList_hipp=''):
- ctx_block = (' '
- '
' )
- sctx_block = (' '
- ''
- '
')
- hipp_block = (
- ' '
- ''
- '
')
-
- # FLAIR
- if 'flair' in featList_ctx:
- flair_ctx_block = (' '
- '
' )
- else:
- flair_ctx_block = ('No cortical FLAIR
')
-
- if 'flair' in featList_sctx:
- flair_sctx_block = (' '
- '
' )
- else:
- flair_sctx_block = ('No subcortical FLAIR
')
-
- if 'flair' in featList_hipp:
- flair_hipp_block = (' '
- '
' )
- else:
- flair_hipp_block = ('No hippocampal FLAIR
')
-
- # qT1
- if 'qt1' in featList_ctx:
- qt1_ctx_block = (' '
- '
' )
- else:
- qt1_ctx_block = ('No cortical qT1
')
-
- if 'qt1' in featList_sctx:
- qt1_sctx_block = (' '
- '
' )
- else:
- qt1_sctx_block = ('No subcortical qT1
')
-
- if 'qt1' in featList_hipp:
- qt1_hipp_block = (' '
- '
' )
- else:
- qt1_hipp_block = ('No hippocampal qT1
')
-
- # ADC
- if 'adc' in featList_ctx:
- adc_ctx_block = (' '
- '
' )
- else:
- adc_ctx_block = ('No cortical ADC
')
-
- if 'adc' in featList_sctx:
- adc_sctx_block = (' '
- '
' )
- else:
- adc_sctx_block = ('No subcortical ADC
')
-
- if 'adc' in featList_hipp:
- adc_hipp_block = (' '
- '
' )
- else:
- adc_hipp_block = ('No hippocampal ADC
')
-
- # THICKNESS
- if 'thickness' in featList_ctx:
- thickness_ctx_block = (' '
- '
' )
- else:
- thickness_ctx_block = ('No cortical atrophy
')
-
- if 'thickness' in featList_sctx:
- thickness_sctx_block = (' '
- '
' )
- else:
- thickness_sctx_block = ('No subcortical atrophy
')
-
- if 'thickness' in featList_hipp:
- thickness_hipp_block = (' '
- '
' )
- else:
- thickness_hipp_block = ('No hippocampal atrophy
')
-
- report_block = (
- # =============================================================================================================#
- # =============================================== P A G E # 1 ===============================================#
- # =============================================================================================================#
- # Title
- '
Clinical report | Regional changes
'
-
- # Subject's ID - Session - Basic demographics
- ' Subject: {subject}, Session: {session}, Sex: {sex}, Age: {age}, '
- 'z-score threshold: {thr}
'
- ' Blue = patient < controls (note: thickness is flipped)
'
- ' Red = patient > controls (note: thickness is flipped)
' + '
' +
-
- # Multivariate cortical changes
- ' Multivariate cortical changes ({featList_ctx})
' + ctx_block +
- ''
-
- ' Multivariate subcortical changes ({featList_sctx})
' + sctx_block + ''
-
- ' Multivariate hippocampal changes ({featList_hipp})
' +
- hipp_block + ''
-
- + '
'
-
- # =============================================================================================================#
- # =============================================== P A G E # 2 ===============================================#
- # ================================================= F L A I R =================================================#
- # =============================================================================================================#
- + '
'
- ' Univariate cortical FLAIR
' + flair_ctx_block + ''
-
- + ' Univariate subcortical FLAIR
'
- + flair_sctx_block + ''
-
- + ' Univariate hippocampal FLAIR
'
- + flair_hipp_block + '' + '
'
-
- # =============================================================================================================#
- # =============================================== P A G E # 3 ===============================================#
- # =================================================== q T 1 ===================================================#
- # =============================================================================================================#
- + '
'
- ' Univariate cortical qT1
' + qt1_ctx_block + ''
-
- + ' Univariate subcortical qT1
'
- + qt1_sctx_block + ''
-
- + ' Univariate hippocampal qT1
'
- + qt1_hipp_block + '' + '
'
-
- # =============================================================================================================#
- # =============================================== P A G E # 4 ===============================================#
- # =================================================== A D C ===================================================#
- # =============================================================================================================#
- + '
'
- ' Univariate cortical ADC
' + adc_ctx_block + ''
-
- + ' Univariate subcortical ADC
'
- + adc_sctx_block + ''
-
- + ' Univariate hippocampal ADC
'
- + adc_hipp_block + '' + '
'
-
- # =============================================================================================================#
- # =============================================== P A G E # 5 ===============================================#
- # ============================================= T H I C K N E S S =============================================#
- # =============================================================================================================#
- + '
'
- ' Univariate cortical atrophy
' + thickness_ctx_block + ''
-
- + ' Univariate subcortical atrophy
'
- + thickness_sctx_block + ''
-
- + ' Univariate hippocampal atrophy
'
- + thickness_hipp_block + '' + '
'
-
- )
-
- return report_block.format(figPath=figPath, subject=subject, session=session, sex=sex, age=age, thr=thr,
- featList_ctx=featList_ctx, featList_sctx=featList_sctx, featList_hipp=featList_hipp)
-
-# Create PDF report
-static_report = ''
-_static_block = report_block_template(figPath, subject=subject, session=session, sex=sex, age=age,
- featList_ctx=featList_ctx, featList_sctx=featList_sctx, featList_hipp=featList_hipp)
-static_report += _static_block
-
-# Utility function
-def convert_html_to_pdf(source_html, output_filename):
- # open output file for writing (truncated binary)
- result_file = open(output_filename, "w+b")
-
- # convert HTML to PDF
- pisa_status = pisa.CreatePDF(
- source_html, # the HTML to convert
- dest=result_file) # file handle to recieve result
-
- # close output file
- result_file.close() # close output file
-
- # return True on success and False on errors
- return pisa_status.err
-
-convert_html_to_pdf(static_report, figPath + '_regionZ_Report.pdf')
-
diff --git a/lost+found/config.cfg b/lost+found/config.cfg
deleted file mode 100644
index dcd06cd..0000000
--- a/lost+found/config.cfg
+++ /dev/null
@@ -1,37 +0,0 @@
-VERSION="v0.0.2 'reborn'"
-
-# Default csv file for controls
-PATH_CSV_CONTROLS='/data_/mica1/03_projects/jessica/hackathon2023/lists/control.csv'
-
-LIST_FEATURES=("thickness" "flair" "FA" "ADC" "qT1")
-LIST_STRUCTURES=("cortex" "subcortex" "hippocampus")
-LIST_TASKS=("proc" "analysis" "report")
-LIST_APPROACHES=("zscore" "norm")
-LIST_RESOLUTIONS=("low" "high")
-LIST_LABELS_CTX=("white" "midthickness" "pial" "swm[0-9]+")
-LIST_LABELS_HIP=("midthickness")
-
-# Subject dir folder names
-FOLDER_LOGS="logs"
-FOLDER_MAPS="maps"
-FOLDER_NORM_Z="norm-z"
-FOLDER_NORM_MODEL="norm-normative"
-
-# Structure folders
-FOLDER_SCTX="subcortex"
-FOLDER_CTX="cortex"
-FOLDER_HIP="hippocampus"
-
-
-# Default smoothing - in mm
-DEFAULT_SMOOTH_CTX=5
-DEFAULT_SMOOTH_HIP=2
-
-# Default threshold
-DEFAULT_THRESHOLD=1.96
-
-# Resolution
-LOW_RESOLUTION_CTX="5k"
-HIGH_RESOLUTION_CTX="32k"
-LOW_RESOLUTION_HIP="0p5mm"
-HIGH_RESOLUTION_HIP="2mm"
diff --git a/lost+found/micapipe-z b/lost+found/micapipe-z
deleted file mode 100755
index 733279c..0000000
--- a/lost+found/micapipe-z
+++ /dev/null
@@ -1,579 +0,0 @@
-#!/bin/bash
-#
-# MICAPIPE-Z BIDS - Normalization and visualization of structural imaging features
-#
-version() {
- echo -e "\z-brains April 2023 (Version v.0.0.2 'reborn')\n"
-}
-export this_version="v0.0.2 'reborn'"
-#---------------- FUNCTION: HELP ----------------#
-help() {
-echo -e "
-\033[38;5;141mCOMMAND:\033[0m
- $(basename $0)
-
-
-\033[38;5;141mARGUMENTS:\033[0m
-\t\033[38;5;197m-sub\033[0m : idBIDS identification
-\t\033[38;5;197m-out\033[0m : Output directory for the processed files .
-\t\033[38;5;197m-bids\033[0m : Path to BIDS directory
-\t\033[38;5;120m-ses\033[0m : OPTIONAL flag that indicates the session name (if omitted will manage as SINGLE session)
-
- Flags for cortical, hippocampal, and subcortical feature post-processing:
-\t\033[38;5;197m-proc_cortex\033[0m : run post-processing of cortical features for report generation. Relies on micapipe derivatives
-\t\t\033[38;5;120m-featureStr\033[0m : Specify which cortical features should undergo post-processing
-\t\t\t Default is '-featureStr all', meaning all features found under /derivatives/micapipe/subject/maps
-\t\t\t To only process one feature, specify using string: e.g. '-featureStr flair'
-\t\t\t Multiple inputs example = '-featureStr flair,FA,ADC'
-\t\t\t Currently supported features: thickness,flair,FA,ADC,qT1
-\t\t\033[38;5;120m-smoothCtx\033[0m : Specify size of gaussian smoothing kernel (Default: 5mm). If set to 0, no smoothing is applied.
-\t\033[38;5;197m-proc_hippocampus\033[0m : run post-processing of hippocampal features for report generation. Relies on hippunfold derivatives
-\t\t\033[38;5;120m-featureStr\033[0m : Specify which hippocampal features should undergo post-processing
-\t\t\t Default is '-featureStr all', meaning all available features based on BIDS raw data
-\t\t\t To only process one feature, specify using string: e.g. '-featureStr flair'
-\t\t\t Multiple inputs example = '-featureStr flair,FA,ADC'
-\t\t\t Currently supported features: thickness,flair,FA,ADC,qT1
-\t\t\033[38;5;120m-smoothHipp\033[0m : Specify size of gaussian smoothing kernel (Default: 2mm). If set to 0, no smoothing is applied.
-\t\033[38;5;197m-proc_subcortex\033[0m : run post-processing of subcortical structures for report generation. Relies on Freesurfer and micapipe
-\t\t\033[38;5;120m-featureStr\033[0m : Specify which subcortical features should undergo post-processing
-\t\t\t Default is '-featureStr all', meaning all available features based on BIDS outputs
-\t\t\t To only process one feature, specify using string: e.g. '-featureStr flair'
-\t\t\t Multiple inputs example = '-featureStr flair,FA,ADC'
-\t\t\t Currently supported features: volume,flair,FA,ADC,qT1
-
- Flags for case-control analysis and clinical report generation:
-\t\033[38;5;197m-regional\033[0m : Multivariate regional analysis
-\t\t\033[38;5;120m-csvControls\033[0m : A list of control subjects
-\t\t\t Default='/data_/mica1/03_projects/jessica/hackathon2023/lists/control.csv'
-\t\t\033[38;5;120m-demo\033[0m : Path to demographic file
-\t\t\t Default=''
-\t\t\033[38;5;120m-featureStr\033[0m : Specify which cortical/subcortical/hippocampal features should be analyzed
-\t\t\t Default=all
-\t\t\t Multiple inputs example = '-featureStr thickness,FA,ADC'
-\t\t\033[38;5;120m-smoothCtx\033[0m : Specify size of gaussian smoothing kernel used to generate cortical input features
-\t\t\t Default=5mm
-\t\t\033[38;5;120m-smoothHipp\033[0m : Specify size of gaussian smoothing kernel used to generate hippocampal input features
-\t\t\t Default=2mm
-\t\033[38;5;197m-asymmetry\033[0m : Multivariate feature asymmetry analysis
-\t\t\033[38;5;120m-thr\033[0m : Z-score threshold for asymmetry maps
-\t\t\t Default=|1.96|
-\t\t\033[38;5;120m-demo\033[0m : Path to demographic file
-\t\t\t Default=''
-\t\t\033[38;5;120m-featureStr\033[0m : Specify which cortical/subcortical/hippocampal features should be analyzed
-\t\t\t Default=all
-\t\t\t Multiple inputs example = '-featureStr thickness,FA,ADC'
-\t\t\033[38;5;120m-smoothCtx\033[0m : Specify size of gaussian smoothing kernel used to generate cortical input features
-\t\t\t Default=5mm
-\t\t\033[38;5;120m-smoothHipp\033[0m : Specify size of gaussian smoothing kernel used to generate hippocampal input features
-\t\t\t Default=2mm
-
-
-\033[38;5;141mOPTIONS:\033[0m
-\t\033[38;5;197m-h|-help\033[0m : Print help
-\t\033[38;5;197m-v|-version\033[0m : Print software version
-\t\033[38;5;197m-force\033[0m : WARNING this will overwrite the idBIDS directory
-\t\033[38;5;197m-quiet\033[0m : Do not print comments
-\t\033[38;5;197m-nocleanup\033[0m : Do not delete temporary directory at script completion
-\t\033[38;5;197m-threads\033[0m : Number of threads (Default is 6)
-\t\033[38;5;197m-tmpDir\033[0m : Specify location of temporary directory (Default is /tmp)
-\t\033[38;5;197m-slim\033[0m : This option will keep only the main outputs and erase all the intermediate files
-\t\033[38;5;197m-mica\033[0m : Only for MICA local processing
-\t\033[38;5;197m-qsub\033[0m : Only for MICA network processing (SGE mica.q)
-\t\033[38;5;197m-qall\033[0m : Only for MICA network processing (SGE all.q)
-
-
-\033[38;5;141mUSAGE:\033[0m
- \033[38;5;141m$(basename $0)\033[0m \033[38;5;197m-sub\033[0m \033[38;5;197m-out\033[0m \033[38;5;197m-bids\033[0m \033[38;5;197m-sctx_vol\033[0m
-
-
-\033[38;5;141mDEPENDENCIES:\033[0m
- > workbench 1.4.2 (https://www.humanconnectome.org/software/workbench-command)
- > ANTS 2.3.3 (https://github.com/ANTsX/ANTs)
- > python 3.7.6 (https://www.python.org)
-
-
-McGill University, MNI, MICA lab, April 2023
-https://github.com/MICA-MNI/micapipe
-https://github.com/MICA-MNI/z-brains
-http://mica-mni.github.io/
-"
-}
-
-# Source utilities functions from MICAPIPE
-if [ -z ${ZBRAINS} ]; then echo "ZBRAINS NOT DEFINED ð¤"; exit 0; fi
-if [ -z ${MICAPIPE} ]; then echo "MICAPIPE NOT DEFINED ð¤"; exit 0; fi
-source "${MICAPIPE}/functions/utilities.sh"
-
-
-# -----------------------------------------------------------------------------------------------#
-# ARGUMENTS
-# Create VARIABLES
-for arg in "$@"
-do
- case "$arg" in
- -h|-help)
- help
- exit 1
- ;;
- -v|-version)
- version
- exit 1
- ;;
- -sub)
- id=$2
- shift;shift
- ;;
- -out)
- out=$2
- shift;shift
- ;;
- -bids)
- BIDS=$2
- shift;shift
- ;;
- -ses)
- SES=$2
- shift;shift
- ;;
- -proc_cortex)
- ctxPROC=TRUE
- shift
- ;;
- -featureStr)
- feat=$2
- shift;shift
- ;;
- -smoothCtx)
- smoothCtx=$2
- shift;shift
- ;;
- -proc_hippocampus)
- hippoPROC=TRUE
- shift
- ;;
- -featureStr)
- feat=$2
- shift;shift
- ;;
- -smoothHipp)
- smoothHipp=$2
- shift;shift
- ;;
- -proc_subcortex)
- sctxPROC=TRUE
- shift
- ;;
- -featureStr)
- feat=$2
- shift;shift
- ;;
- -regional)
- REGZ=TRUE
- shift
- ;;
- -csvControls)
- csvControls=$2
- shift;shift
- ;;
- -demo)
- demo=$2
- shift;shift
- ;;
- -smoothCtx)
- smoothCtx=$2
- shift;shift
- ;;
- -smoothHipp)
- smoothHipp=$2
- shift;shift
- ;;
- -asymmetry)
- ASYMMETRY=TRUE
- shift
- ;;
- -thr)
- thr=$2
- shift;shift
- ;;
- -demo)
- demo=$2
- shift;shift
- ;;
- -featureStr)
- feat=$2
- shift;shift
- ;;
- -smoothCtx)
- smoothCtx=$2
- shift;shift
- ;;
- -smoothHipp)
- smoothHipp=$2
- shift;shift
- ;;
- -mica)
- mica=TRUE
- shift
- ;;
- -tmpDir)
- tmpDir=$2
- shift;shift;
- ;;
- -qsub)
- micaq=TRUE
- shift
- ;;
- -qall)
- qall=TRUE
- shift
- ;;
- -nocleanup)
- nocleanup=TRUE
- shift
- ;;
- -threads)
- threads=$2
- shift;shift
- ;;
- -force)
- force=TRUE
- shift
- ;;
- -quiet)
- quiet=TRUE
- shift
- ;;
- -*)
- Error "Unknown option ${2}"
- help
- exit 1
- ;;
- esac
-done
-
-# argument check out & WARNINGS
-arg=($id $out $BIDS)
-if [ ${#arg[@]} -lt 3 ]; then
-Error "One or more mandatory arguments are missing:
- -sub : $id
- -out : $out
- -bids : $BIDS"
-help; exit 1; fi
-runs=($ALL $ctxPROC $hippoPROC $sctxPROC $REGZ $ASYMMETRY)
-if [ "${#runs[@]}" -lt 1 ]; then
-Error "A processing flag is missing:
- -proc_cortex
- -proc_hippocampus
- -proc_subcortex
- -regional
- -asymmetry"
-help; exit 1; fi
-
-# Get the real path of the Inputs
-outDeriv=$(realpath $out)
-outz=$(realpath $out)/zbrains
-outm=$(realpath $out)/micapipe_v0.2.0
-BIDS=$(realpath $BIDS)
-id=${id/sub-/}
-here=$(pwd)
-
-# Number of session (Default is "ses-pre")
-if [ -z ${SES} ]; then SES="SINGLE"; else SES="ses-${SES/ses-/}"; fi
-
-# Assigns variables names
-bids_variables "$BIDS" "$id" "$outm" "$SES"
-
-# Check BIDS Directory
-if [ ! -d "${subject_bids}" ]; then Error "$id was not found on the BIDS directory\n\t Check ls ${subject_bids}"; exit 1; fi
-
-# Erase temporary files by default
-if [ -z ${nocleanup} ]; then nocleanup=FALSE; fi
-
-# No print Do_cmd
-if [ "$quiet" = "TRUE" ]; then export quiet=TRUE; fi
-
-# Temporary directory
-if [ -z ${tmpDir} ]; then export tmpDir=/tmp; else tmpDir=$(realpath $tmpDir); fi
-
-# Setting Surface Directory
-Nrecon=($(ls "${dir_QC}/${idBIDS}_module-proc_surf-"*.json 2>/dev/null | wc -l))
-if [[ "$Nrecon" -lt 1 ]]; then
- Error "Subject $id doesn't have a module-proc_surf: run -proc_surf"; exit 1
-elif [[ "$Nrecon" -eq 1 ]]; then
- module_qc=$(ls "${dir_QC}/${idBIDS}_module-proc_surf-"*.json 2>/dev/null)
- recon="$(echo ${module_qc/.json/} | awk -F 'proc_surf-' '{print $2}')"
-elif [[ "$Nrecon" -gt 1 ]]; then
- Warning "${idBIDS} has been processed with freesurfer and fastsurfer."
- if [[ "$FastSurfer" == "TRUE" ]]; then
- Note "fastsurfer will run: $FastSurfer\n"; recon="fastsurfer";
- else
- Note "freesurfer is the default"; recon="freesurfer"
- fi
-fi
-# Set surface directory
-set_surface_directory "${recon}"
-if [ -z ${dir_subjsurf} ]; then FSdir="FALSE"; else FSdir=$(realpath $dir_subjsurf); fi
-
-
-# -----------------------------------------------------------------------------------------------#
-# Set optional arguments
-
-# Feature selection and smoothing
-if [ -z ${feat} ]; then featStr='all'; else featStr="$feat"; fi
-if [ -z ${smoothCtx} ]; then smoothCtx='5mm'; else smoothingCTX="$smoothCTX"; fi
-if [ -z ${smoothHipp} ]; then smoothHipp='2mm'; else smoothingHIPP="$smoothHIPP"; fi
-
-# arguments for regional and asymmetry analysis
-if [ -z "$csvControls" ]; then csvControls="/data_/mica1/03_projects/jessica/hackathon2023/lists/control.csv"; fi
-if [ -z "$demo" ]; then demo=''; fi
-
-# -----------------------------------------------------------------------------------------------#
-Title "z-brains pipeline - (Version $this_version) \n\t\tidBIDS: $id Session: $SES"
-
-# -----------------------------------------------------------------------------------------------#
-# Launch the init file for local processing at MICA lab
-if [ "$micaq" = "TRUE" ] || [ "$qall" = "TRUE" ] || [ "$mica" = "TRUE" ]; then
- source "${MICAPIPE}/functions/init.sh"
-else
-
- # -----------------------------------------------------------------------------------------------#
- # CHECK PACKAGES DEPENDENCIES
- # tree display
- if [[ -z $(which tree) ]]; then Warn "tree function was not found"; fi
- # freesurfer
- if [[ -z $(which recon-all) ]]; then Error "FreeSurfer was not found"; exit 1; fi
- # FSL
- if [[ -z $(which flirt) ]]; then Error "FSL was not found"; exit 1; fi
- # ANTSx
- if [[ -z $(which antsRegistration) ]]; then Error "ANTs was not found"; exit 1; fi
- # workbench
- if [[ -z $(which wb_command) ]]; then Error "WorkBench was not found"; exit 1; fi
-fi
-
-# Processing
-if [[ -z $PROC ]]; then export PROC="LOCAL"; fi
-
-# Number of THREADS used by ANTs and workbench, default is 6 if not defined by -threads
-if [[ -z $threads ]]; then export threads=6; fi
-Info "z-brains will use $threads threads for multicore processing"
-
-# Directories check
-if [[ ${force} == TRUE ]]; then
- Warning "$id processing directory will be overwritten"
- rm -rf $outz/${idBIDS};
-fi
-
-# -----------------------------------------------------------------------------------------------#
-# Timer & Beginning
-aloita=$(date +%s)
-
-# Create tmp dir
-if [ ! -d ${tmpDir} ]; then Do_cmd mkdir -p $tmpDir; fi
-
-# Create output directories
-if [ ! -d ${outz} ]; then Do_cmd mkdir -p ${outz}; fi
-
-# Handle Single Session
-if [ "$SES" == "SINGLE" ]; then
- export subject_dirz=$outz/${subject}
- ses=""
-else
- export subject_dirz=$outz/${subject}/${SES}
- ses="_${SES}"
-fi
-
-# subject_dir
-
-# Creates subject directory
-Do_cmd mkdir -m 770 -p "${subject_dirz}"/{logs,maps,norm-z,norm-propensity,norm-normative}
-Do_cmd mkdir -m 770 -p "${subject_dirz}"/maps/{subcortex,cortex,hippocampus}
-Do_cmd mkdir -m 770 -p "${subject_dirz}"/norm-z/{subcortex,cortex,hippocampus}
-Do_cmd mkdir -m 770 -p "${subject_dirz}"/norm-propensity/{subcortex,cortex,hippocampus}
-Do_cmd mkdir -m 770 -p "${subject_dirz}"/norm-normative/{subcortex,cortex,hippocampus}
-chmod g+xX -R --silent "${subject_dirz}"
-
-# -----------------------------------------------------------------------------------------------#
-# subject_dir
-if [ ! -d "${subject_dir}" ]; then
- Error "Subject ${id} directory doesn't exist, run micapipe v0.2.2 before running z-brains"
-exit 1; fi
-
-# freesurfer Directory
-if [ ! -d "$FSdir" ]; then
- Error "Subject ${id} freesurfer directory doesn't exist, run micapipe v0.2.0 before running z-brains"
-exit 1; fi
-
-# Log directory
-outLogs="${subject_dirz}/logs/"
-
-# Pipeline description json
-micapipe_json
-
-
-# -----------------------------------------------------------------------------------------------#
-# # # Structural processing: Cortex
-# -----------------------------------------------------------------------------------------------#
-if [ "$ALL" = "TRUE" ] || [ "$ctxPROC" = "TRUE" ]; then
- log_file_str=$outLogs/ctx_proc_$(date +'%d-%m-%Y')
-
- # Define UTILITIES directory
- export scriptDir=${ZBRAINS}/functions
-
- COMMAND="${scriptDir}/01_ctx_proc.sh $BIDS $id $out $SES $FSdir $nocleanup $threads $tmpDir $featStr $smoothCtx"
- # mica.q - Structural processing: Cortex
- if [[ $micaq == "TRUE" ]]; then
- Info "MICA qsub - Structural processing: Cortex"
- Warning "This option only works on the MICA workstations"
- export PROC="qsub-MICA"
- # Info "Setting a warper to run on mica.q"
- qsub -q mica.q -pe smp $threads -l h_vmem=6G -N q${id}_ctx_proc -e ${log_file_str}.e -o ${log_file_str}.txt $COMMAND $PROC
- # all.q - Structural processing: Cortex
- elif [[ $qall == "TRUE" ]]; then
- Info "all.q qsub - Structural processing: Cortex"
- Warning "This option only works on the McGill BIC network"
- export PROC="qsub-all.q"
- # Info "Setting a warper to run on all.q"
- qsub -q all.q -pe all.pe $threads -l h_vmem=6G -N q${id}_ctx_proc -e ${log_file_str}.e -o ${log_file_str}.txt $COMMAND $PROC
- else
- # LOCAL - Structural processing: Cortex
- $COMMAND $PROC 2>&1 | tee -a ${log_file_str}.txt
- fi
-fi
-
-
-# -----------------------------------------------------------------------------------------------#
-# # # Structural processing: Hippocampus
-# -----------------------------------------------------------------------------------------------#
-if [ "$ALL" = "TRUE" ] || [ "$hippoPROC" = "TRUE" ]; then
- log_file_str=$outLogs/hippo_proc_$(date +'%d-%m-%Y')
-
- # Define UTILITIES directory
- export scriptDir=${ZBRAINS}/functions
-
- COMMAND="${scriptDir}/01_hippo_proc.sh $BIDS $id $out $SES $FSdir $nocleanup $threads $tmpDir $featStr $smoothHipp"
- # mica.q - Structural processing: Hippocampus
- if [[ $micaq == "TRUE" ]]; then
- Info "MICA qsub - Structural processing: Hippocampus"
- Warning "This option only works on the MICA workstations"
- export PROC="qsub-MICA"
- # Info "Setting a warper to run on mica.q"
- qsub -q mica.q -pe smp $threads -l h_vmem=6G -N q${id}_hippo_proc -e ${log_file_str}.e -o ${log_file_str}.txt $COMMAND $PROC
- # all.q - Structural processing: Hippocampus
- elif [[ $qall == "TRUE" ]]; then
- Info "all.q qsub - Structural processing: Hippocampus"
- Warning "This option only works on the McGill BIC network"
- export PROC="qsub-all.q"
- # Info "Setting a warper to run on all.q"
- qsub -q all.q -pe all.pe $threads -l h_vmem=6G -N q${id}_hippo_proc -e ${log_file_str}.e -o ${log_file_str}.txt $COMMAND $PROC
- else
- # LOCAL - Structural processing: Hippocampus
- $COMMAND $PROC 2>&1 | tee -a ${log_file_str}.txt
- fi
-fi
-
-
-# -----------------------------------------------------------------------------------------------#
-# # # Structural processing: Subcortex
-# -----------------------------------------------------------------------------------------------#
-if [ "$ALL" = "TRUE" ] || [ "$sctxPROC" = "TRUE" ]; then
- log_file_str=$outLogs/sctx_proc_$(date +'%d-%m-%Y')
-
- # Define UTILITIES directory
- export scriptDir=${ZBRAINS}/functions
-
- COMMAND="${scriptDir}/01_sctx_proc.sh $BIDS $id $outDeriv $SES $FSdir $nocleanup $threads $tmpDir $featStr"
- # mica.q - Structural processing: Subcortex
- if [[ $micaq == "TRUE" ]]; then
- Info "MICA qsub - Structural processing: Subcortex"
- Warning "This option only works on the MICA workstations"
- export PROC="qsub-MICA"
- # Info "Setting a warper to run on mica.q"
- qsub -q mica.q -pe smp $threads -l h_vmem=6G -N q${id}_sctx_proc -e ${log_file_str}.e -o ${log_file_str}.txt $COMMAND $PROC
- # all.q - Structural processing: Subcortex
- elif [[ $qall == "TRUE" ]]; then
- Info "all.q qsub - Structural processing: Subcortex"
- Warning "This option only works on the McGill BIC network"
- export PROC="qsub-all.q"
- # Info "Setting a warper to run on all.q"
- qsub -q all.q -pe all.pe $threads -l h_vmem=6G -N q${id}_sctx_proc -e ${log_file_str}.e -o ${log_file_str}.txt $COMMAND $PROC
- else
- # LOCAL - Structural processing: Subcortex
- $COMMAND $PROC 2>&1 | tee -a ${log_file_str}.txt
- fi
-fi
-
-
-# -----------------------------------------------------------------------------------------------#
-# # # Multivariate analysis: REGIONAL Z
-# -----------------------------------------------------------------------------------------------#
-if [ "$ALL" = "TRUE" ] || [ "$REGZ" = "TRUE" ]; then
- log_file_str=$outLogs/regional_$(date +'%d-%m-%Y')
-
- export scriptDir=${ZBRAINS}/functions
-
- COMMAND="python ${scriptDir}/02_zscores.py $csvControls $id $featStr $smoothCtx $smoothHipp $outz $SES"
- echo $COMMAND
- # mica.q - Multivariate analysis: REGIONAL Z
- if [[ $micaq == "TRUE" ]]; then
- Info "MICA qsub - Multivariate analysis: REGIONAL"
- Warning "This option only works on the MICA workstations"
- export PROC="qsub-MICA"
- # Info "Setting a warper to run on mica.q"
- qsub -q mica.q -pe smp $threads -l h_vmem=6G -N q${id}_regionZ -e ${log_file_str}.e -o ${log_file_str}.txt $COMMAND $PROC
- # all.q - Multivariate analysis: REGIONAL Z
- elif [[ $qall == "TRUE" ]]; then
- Info "all.q qsub - Multivariate analysis: REGIONAL"
- Warning "This option only works on the McGill BIC network"
- export PROC="qsub-all.q"
- # Info "Setting a warper to run on all.q"
- qsub -q all.q -pe all.pe $threads -l h_vmem=6G -N q${id}_regionZ -e ${log_file_str}.e -o ${log_file_str}.txt $COMMAND $PROC
- else
- # LOCAL - Multivariate analysis: REGIONAL Z
- $COMMAND $PROC 2>&1 | tee -a ${log_file_str}.txt
- fi
-fi
-
-# -----------------------------------------------------------------------------------------------#
-# # # Multivariate analysis: ASYMMETRY
-# -----------------------------------------------------------------------------------------------#
-if [ "$ALL" = "TRUE" ] || [ "$ASYMMETRY" = "TRUE" ]; then
- log_file_str=$outLogs/asymmetry_$(date +'%d-%m-%Y')
-
- export scriptDir=${ZBRAINS}/functions
-
- COMMAND="${scriptDir}/02_asymmetry.sh $BIDS $id $out $SES $nocleanup $threads $tmpDir $thr $demo $feat $smoothingCTX $smoothingHIPP"
- # mica.q - Multivariate analysis: ASYMMETRY
- if [[ $micaq == "TRUE" ]]; then
- Info "MICA qsub - Multivariate analysis: ASYMMETRY"
- Warning "This option only works on the MICA workstations"
- export PROC="qsub-MICA"
- # Info "Setting a warper to run on mica.q"
- qsub -q mica.q -pe smp $threads -l h_vmem=6G -N q${id}_asymmetry -e ${log_file_str}.e -o ${log_file_str}.txt $COMMAND $PROC
- # all.q - Multivariate analysis: ASYMMETRY
- elif [[ $qall == "TRUE" ]]; then
- Info "all.q qsub - Multivariate analysis: ASYMMETRY"
- Warning "This option only works on the McGill BIC network"
- export PROC="qsub-all.q"
- # Info "Setting a warper to run on all.q"
- qsub -q all.q -pe all.pe $threads -l h_vmem=6G -N q${id}_asymmetry -e ${log_file_str}.e -o ${log_file_str}.txt $COMMAND $PROC
- else
- # LOCAL - Multivariate analysis: ASYMMETRY
- $COMMAND $PROC 2>&1 | tee -a ${log_file_str}.txt
- fi
-fi
-
-: <<'END'
-END
-
-
-# -----------------------------------------------------------------------------------------------#
-# Total Running Time
-lopuu=$(date +%s)
-eri=$(echo "$lopuu - $aloita" | bc)
-eri=$(echo print $eri/60 | perl)
-
-# Cleanup if processing was local
-if [ $PROC == "LOCAL" ] || [ "$mica" = "TRUE" ]; then
- cleanup "$tmpDir/*_micapipe*_${id}" "$nocleanup" "$here"
-fi
-Title "GLOBAL z-brains running time with $PROC processing:\033[38;5;220m $(printf "%0.3f\n" ${eri}) minutes \033[38;5;141m"
diff --git a/lost+found/propensity_scores.py b/lost+found/propensity_scores.py
deleted file mode 100644
index b6d9f56..0000000
--- a/lost+found/propensity_scores.py
+++ /dev/null
@@ -1,87 +0,0 @@
-"""
-Propensity scores estimation and matching.
-"""
-
-import numpy as np
-import pandas as pd
-
-from sklearn.preprocessing import StandardScaler, OneHotEncoder
-from sklearn.compose import make_column_transformer
-from sklearn.calibration import CalibratedClassifierCV
-from sklearn.linear_model import LogisticRegression
-from sklearn.pipeline import Pipeline
-
-
-def estimate_ps(df_cov: pd.DataFrame, y: np.ndarray, cat: list[str] | None = None) -> np.ndarray:
- """Estimate propensity scores from covariates and class label.
-
- Parameters
- ----------
- df_cov: pd.DataFrame of shape=(n_subjects, n_covariates)
- DataFrame holding subjects covariates.
- y: np.ndarray[int], shape=(n_subjects,)
- Binary class label (e.g., diagnosis)
- cat: list of str, default=None
- Categorical covariates in `df_cov` (e.g., site). If None, assumes no categorical
- columns in dataframe.
-
- Returns
- -------
- ps: np.ndarray of shape (n_subjects,)
- Estimated propensity scores.
- """
-
- cat = [] if cat is None else cat
-
- scale_keys = np.setdiff1d(df_cov.columns, cat)
- ct = make_column_transformer((StandardScaler(), scale_keys),
- (OneHotEncoder(drop='first'), cat))
-
- logit = LogisticRegression(penalty=None, random_state=0)
- clf = Pipeline([('ct', ct), ('clf', CalibratedClassifierCV(logit))])
- ps = clf.fit(df_cov, y).predict_proba(df_cov)[:, 1]
-
- return ps
-
-
-def get_matches(ps: float, ps_cn: np.ndarray[float], caliper: float | None = .2, n_min: int | None = None,
- n_max: int | None = None):
- """ Get the indices of the closest propensity scores in `pc_cn` to `ps`.
-
- Parameters
- ----------
- ps: float
- Propensity score of target subject
- ps_cn: np.ndarray, shape = (n_subjects,)
- Propensity scores of control subjects
- caliper: float, default=0.2
- Caliper to use for imperfect matches. If None, no caliper is used.
- n_min: int, default=None
- Minimum number of matches.
- n_max: int default=None
- Maximum number of matches to consider.
-
- Returns
- -------
- matches: np.ndarray, shape=(n_matches,)
- Indices of the closest matches in `pc_cn`.
- """
- d = np.abs(ps_cn - ps)
-
- # Indices of closest matches in controls
- idx = np.argsort(d)
-
- if caliper is not None:
- thresh = caliper * d.std()
- n = np.count_nonzero(d < thresh)
- if n_min is not None:
- n_min = max(n, n_min)
- else:
- n_min = n
- elif n_min is None:
- n_min = d.size
-
- if n_max is not None:
- n_min = min(n_min, n_max)
-
- return idx[:n_min]