From 3494efa280c24023c9be74727f5257df46487658 Mon Sep 17 00:00:00 2001 From: rcruces Date: Fri, 8 Mar 2024 16:57:34 -0500 Subject: [PATCH] removed Lost+found updated gitignore --- lost+found/02_asymmetry-report.py | 272 ---------- lost+found/02_asymmetry.py | 462 ----------------- lost+found/02_asymmetry.sh | 173 ------- lost+found/02_regional_analysis.py | 404 --------------- lost+found/02_regional_analysis.sh | 176 ------- lost+found/02_regional_analysis_report.py | 272 ---------- lost+found/config.cfg | 37 -- lost+found/micapipe-z | 579 ---------------------- lost+found/propensity_scores.py | 87 ---- 9 files changed, 2462 deletions(-) delete mode 100644 lost+found/02_asymmetry-report.py delete mode 100644 lost+found/02_asymmetry.py delete mode 100755 lost+found/02_asymmetry.sh delete mode 100644 lost+found/02_regional_analysis.py delete mode 100755 lost+found/02_regional_analysis.sh delete mode 100644 lost+found/02_regional_analysis_report.py delete mode 100644 lost+found/config.cfg delete mode 100755 lost+found/micapipe-z delete mode 100644 lost+found/propensity_scores.py 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]