<jupyter_start><jupyter_text>Create and test the files/funcitons needed for the fMRIPrep of PPMI datasetVincent  1. Subject-session csv files for sdMRI/livingPark (get_subj_ses.py)ppmi_subject_session.csvppmi_subject_session_part2.csv<jupyter_code># libs and envs
import sys
from pathlib import Path
import pandas as pd
import numpy as np

#proj_dir='/data/pd/ppmi' # BIC
proj_dir='/scratch' # CC
codes_dir_str=proj_dir+'/mr_proc'
sys.path.append(codes_dir_str)

# main PATH
codes_dir    = Path(codes_dir_str)
fmriprep_dir = codes_dir / 'workflow'/ 'fMRIPrep' 

# output subject session list 
col_names = ['subject', 'session']
ppmi_subj_ses_file   = fmriprep_dir / 'ppmi_subject_session.csv'  # Information from download database.
# Read and check the existing subject session list 
subj_ses_exist_df = pd.read_csv(ppmi_subj_ses_file, sep=',', header=None, index_col=None, names= col_names)
print('sessions in dataset: ', 'ses-'+', ses-'.join([str(x) for x in list(subj_ses_exist_df.session.unique())]))
#subj_ses_exist_df<jupyter_output>sessions in dataset:  ses-1, ses-7, ses-21, ses-5, ses-11, ses-91, ses-90, ses-0, ses-9, ses-30<jupyter_text>2. Create bids-filter for all the sessions: mr_proc\workflow\fMRIPrep\anat_ses_\.json<jupyter_code># create json session files
save_file=0

import json
basic_json_filename = 'anat_ses-?.json'
session_json_dict={
    "t1w": {
        "datatype": "anat",
        "session": "0",
	"suffix": "T1w"
    }
}

for x in list(subj_ses_exist_df.session.unique()):
    ses_str=str(x)
    json_file_tmp = basic_json_filename.replace('?', ses_str)
    print(json_file_tmp)
    session_json_dict['t1w']['session']=ses_str
    json_string_tmp = json.dumps(session_json_dict, indent=4)
    #print(json_string_tmp)
    if save_file==1:
        with open((fmriprep_dir/json_file_tmp), 'w', encoding='utf-8') as outfile_tmp:
            outfile_tmp.write(json_string_tmp)<jupyter_output>anat_ses-1.json
anat_ses-7.json
anat_ses-21.json
anat_ses-5.json
anat_ses-11.json
anat_ses-91.json
anat_ses-90.json
anat_ses-0.json
anat_ses-9.json
anat_ses-30.json<jupyter_text>3. Subject-session csv files for livingparklivingpark_subject_session.csv<jupyter_code>save_file=0

livingpark_subj_metadata_file = codes_dir / 'tab_data' /'PPMI_livingpark_dcminfo.csv' # Information from download database.
livingpark_subj_session_file   = fmriprep_dir / 'livingpark_subject_session.csv'  

livingpark_col_names = ['subject', 'session']
livingpark_subj_ses_df = pd.read_csv(livingpark_subj_metadata_file, sep=',')
livingpark_subj_ses_df = livingpark_subj_ses_df[livingpark_col_names].drop_duplicates()
livingpark_subj_ses_df['Subject']=['sub-'+str(x) for x in livingpark_subj_ses_df['Subject']]
#print('sessions in dataset: ', 'ses-'+', ses-'.join([str(x) for x in list(subj_ses_exist_df.Visit.unique())]))
if save_file:
    livingpark_subj_ses_df.to_csv(livingpark_subj_session_file, index=False, header=False)
#livingpark_subj_ses_df<jupyter_output><empty_output><jupyter_text>4. Rerun sdMRI fmriprep failure subjectssdMRI_subject_session_rerun1.csv<jupyter_code>save_file=0

sdMRI_err_file = fmriprep_dir /'err_fmriprep_run1.log' # Information from download database.
sdMRI_rerun1_ses_file = fmriprep_dir / 'sdMRI_subject_session_rerun1.csv'  

subj_ses_col_names = ['subject', 'session']
sdMRI_rerun1_df = pd.read_csv(sdMRI_err_file, names=subj_ses_col_names, header=None, sep='_', usecols=[2,3])
sdMRI_rerun1_df = sdMRI_rerun1_df.drop_duplicates()

if save_file:
    sdMRI_rerun1_df.to_csv(sdMRI_rerun1_ses_file, index=False, header=False)
#sdMRI_rerun1_df<jupyter_output><empty_output><jupyter_text>5. simple QC scripts    1. Reference to the description of fMRIPrep docs [https://fmriprep.org/en/stable/outputs.html]    2. test code for: mr_proc/workflow/fMRIPrep/fmriprep_simple_qc.py<jupyter_code># Test code block for fmriprep_simple_qc.py
# libs and envs
import sys
from pathlib import Path
import pandas as pd
import numpy as np
import os
import glob
import shutil

# inputs
root_str = "/scratch/PPMI_ses-21_fmriprep_anat_20.2.7"
out_report_dir = '/scratch' # CC
out_csv_dir = out_report_dir+'/mr_proc/workflow/fMRIPrep'

root_dir = Path(root_str)
fmriprep_dir = root_dir / "fmriprep"
fs_dir = root_dir / "freesurfer-6.0.1"

## general fMRIPrep info
root_folder = root_str.split('/')[-1]
[dataset_name, ses_full_str, software_name, proc_name, software_version]=root_folder.split('_')
# Session info
ses_name=ses_full_str.split('-')[-1]
fs_ver = str(fs_dir).split('/')[-1].split('-')[-1]

# qc report folder
report_folder_name=root_folder+"_report"
qc_report_out_dir = Path(out_report_dir+'/'+report_folder_name)
if (not qc_report_out_dir.is_dir()): os.mkdir(qc_report_out_dir) 
else: shutil.rmtree(qc_report_out_dir); os.mkdir(qc_report_out_dir) 
# qc summary table file
qc_tab_file = Path(out_csv_dir+'/'+root_folder+"_report.csv")

# get subject list
sub_withFolder=[x for x in next(os.walk(fmriprep_dir))[1] if "sub" in x]
sub_withReport=[x.split('.')[0] for x in next(os.walk(fmriprep_dir))[2] if x.split('.')[-1]=="html"]
sub_fsFolder=[x for x in next(os.walk(fs_dir))[1] if "sub" in x]
sub_list=list(set(sub_withFolder) or set(sub_reported) or set(sub_fsFolder))
n_sub=len(sub_list)
print(str(n_sub)+" fMRIPreped (version "+software_version+") subjects found for Dataset -> "+dataset_name+" Session -> "+ses_name)
# create qc summary df
qc_df = pd.DataFrame({'dataset':[dataset_name]*n_sub,'session':[ses_name]*n_sub,'fmriprep_proc':[proc_name]*n_sub,'fMRIPrep_ver':[software_version]*n_sub}, index=sub_list)
for x in sub_withFolder: qc_df.loc[x,'sub_fMRIPrep']=1
for x in sub_withReport: qc_df.loc[x,'sub_report']=1
for x in sub_fsFolder: qc_df.loc[x,'sub_Freesurfer']=1
for x in sub_fsFolder: qc_df.loc[x,'Freesurfer_ver']=fs_ver

#QC image list
fMRIPrep_qc_dict = {'t1_image': 'desc-preproc_T1w.nii.gz',
                    't1_mask' : 'desc-brain_mask.nii.gz',
                    'aparcaseg' : 'desc-aparcaseg_dseg.nii.gz',
                    'aseg': 'desc-aseg_dseg.nii.gz',
                    'dseg': 'dseg.nii.gz',
                    'prob_csf': 'label-CSF_probseg.nii.gz',
                    'prob_gm' : 'label-GM_probseg.nii.gz', 
                    'prob_wm' : 'label-WM_probseg.nii.gz',
                    'mni_T1'  : 'space-MNI152NLin2009cAsym_res-2_desc-preproc_T1w.nii.gz',
                    'mni_mask': 'space-MNI152NLin2009cAsym_res-2_desc-brain_mask.nii.gz', 
                    'mni_dseg': 'space-MNI152NLin2009cAsym_res-2_dseg.nii.gz',
                    'mni_csf' : 'space-MNI152NLin2009cAsym_res-2_label-CSF_probseg.nii.gz',
                    'mni_gm'  : 'space-MNI152NLin2009cAsym_res-2_label-GM_probseg.nii.gz', 
                    'mni_wm'  : 'space-MNI152NLin2009cAsym_res-2_label-WM_probseg.nii.gz'
                   }
Freesurfer_qc_dict = {'fs_aseg'  : 'aseg.stats',
                      'fs_wmparc':'wmparc.stats',
                      # left hemisphere
                      'lh.a2009s' : 'lh.aparc.a2009s.stats',
                      'lh.DKT'    : 'lh.aparc.DKTatlas.stats',
                      'lh.pial'   : 'lh.aparc.pial.stats',
                      'lh.aparc'  : 'lh.aparc.stats',
                      'lh.BA_exvivo': 'lh.BA_exvivo.stats',
                      'lh.BA_exvivo_th' : 'lh.BA_exvivo.thresh.stats', 
                      'lh.curv' : 'lh.curv.stats',
                      'lh.w-g.pct' : 'lh.w-g.pct.stats',
                      # right hemisphere
                      'rh.a2009s': 'rh.aparc.a2009s.stats', 
                      'rh.DKT'   : 'rh.aparc.DKTatlas.stats',
                      'rh.pial'  : 'rh.aparc.pial.stats',
                      'rh.aparc' : 'rh.aparc.stats',
                      'rh.BA_exvivo' : 'rh.BA_exvivo.stats',
                      'rh.BA_exvivo_th' : 'rh.BA_exvivo.thresh.stats', 
                      'rh.curv' : 'rh.curv.stats',
                      'rh.w-g.pct' : 'rh.w-g.pct.stats'
                     }
# set run
current_run_str = '_run-1_';
for sub_id_str in sub_list:
    # copy out fMRIPrep report
        if (fmriprep_dir/sub_id_str).is_dir():
            qc_df.loc[sub_id_str, 'sub_fMRIPrep']=int(1)
            if (fmriprep_dir/sub_id_str/"figures").is_dir():
                shutil.copytree(fmriprep_dir/sub_id_str/"figures", qc_report_out_dir/sub_id_str/"figures", dirs_exist_ok=True)
            if (fmriprep_dir/sub_id_str/"log").is_dir():
                shutil.copytree(fmriprep_dir/sub_id_str/"log", qc_report_out_dir/sub_id_str/"log", dirs_exist_ok=True)
        else:
            qc_df.loc[sub_id_str, 'sub_fMRIPrep']=int(0)
        # html file
        if (fmriprep_dir/(sub_id_str+'.html')).is_file():
            qc_df.loc[sub_id_str, 'sub_report']=int(1)
            shutil.copy2(fmriprep_dir/(sub_id_str+'.html'), qc_report_out_dir/(sub_id_str+'.html'))
        else:
            qc_df.loc[sub_id_str, 'sub_report']=int(0)
    ## qc fmriprep
    sub_img_dir = fmriprep_dir/sub_id_str/ses_full_str/'anat'
    for k, v in fMRIPrep_qc_dict.items():
        current_fmriprep_file_name = sub_id_str+'_'+ses_full_str+current_run_str+v
        current_fmriprep_file = sub_img_dir/current_fmriprep_file_name
        if current_fmriprep_file.is_file():
            qc_df.loc[sub_id_str, k]=int(1)
        else:
            qc_df.loc[sub_id_str, k]=int(0)
            print(sub_id_str+' missing fMRIPrep results: '+current_fmriprep_file_name)
    ## qc freesurfer
    sub_fs_dir = fs_dir/sub_id_str/'stats'
    for k, v in Freesurfer_qc_dict.items():
        #print(str(sub_img_dir/(sub_id_str+'_'+ses_full_str+current_run_str+v)))
        current_fmriprep_file = sub_fs_dir/v
        if current_fmriprep_file.is_file():
            qc_df.loc[sub_id_str, k]=int(1)
        else:
            qc_df.loc[sub_id_str, k]=int(0)
            print(sub_id_str+' missing Freesurfer results: '+v)
            
## retrieve qc reports
# zip fMRIPrep QC report
shutil.make_archive(out_report_dir+'/'+report_folder_name, 'zip', str(qc_report_out_dir))
# save QC report
qc_df.loc[:, 'fmriprep_rerun']=qc_df.loc[:,fMRIPrep_qc_dict.keys()].apply(lambda x: int(not all(x)), axis=1)
qc_df.loc[:, 'freesurfer_rerun']=qc_df.loc[:,Freesurfer_qc_dict.keys()].apply(lambda x: int(not all(x)), axis=1)
ordered_columns = ['dataset', 'session', 'fmriprep_proc', 'fMRIPrep_ver', 'sub_fMRIPrep',
                   'sub_report', 'sub_Freesurfer', 'Freesurfer_ver', 'fmriprep_rerun', 'freesurfer_rerun', 
                   't1_image', 't1_mask','aparcaseg', 'aseg', 'dseg', 'prob_csf', 'prob_gm', 'prob_wm', 
                   'mni_T1', 'mni_mask', 'mni_dseg', 'mni_csf', 'mni_gm', 'mni_wm', 
                   'fs_aseg', 'fs_wmparc',
                   'lh.a2009s', 'lh.DKT', 'lh.pial', 'lh.aparc', 'lh.BA_exvivo', 'lh.BA_exvivo_th', 'lh.curv', 'lh.w-g.pct', 
                   'rh.a2009s', 'rh.DKT', 'rh.pial', 'rh.aparc', 'rh.BA_exvivo', 'rh.BA_exvivo_th', 'rh.curv', 'rh.w-g.pct']
qc_df=qc_df[ordered_columns]
rerun_list=qc_df[(qc_df['fmriprep_rerun']==1)|(qc_df['freesurfer_rerun']==1)].index
if len(rerun_list)==0:
    print(root_folder, " passed simple QC!")
else:
    print('Subjects need to rerun for',root_folder, ':',rerun_list)
qc_df.to_csv(qc_tab_file)
qc_df<jupyter_output>10 fMRIPreped (version 20.2.7) subjects found for Dataset -> PPMI Session -> 21<jupyter_text>6. Create rerun sub-session list after simple QC scripts    1. Input:  The simple fMRIPrep QC results for all sessions;    2. Output: The sub-ses csv rerun_list;<jupyter_code># Test code block for fmriprep_simple_qc.py
# libs and envs
save_file=1
import sys
import pandas as pd

# inputs
out_csv_dir = '/scratch/mr_proc/workflow/fMRIPrep'

# output
sub_ses_rerun2_file = '/scratch/mr_proc/workflow/fMRIPrep/sdMRI_subject_session_rerun2.csv'

ses_list = ['0', '1', '5', '7', '9', '11', '21', '30', '90', '91'];
ses_df = [pd.read_csv(out_csv_dir+'/PPMI_ses-'+x+'_fmriprep_anat_20.2.7_report.csv')for x in ses_list]
rerun2_df = pd.concat(ses_df)
rerun2_df=rerun2_df[(rerun2_df['fmriprep_rerun']==1)|(rerun2_df['Freesurfer_ver']==1)].copy()
sdMRI_rerun2_df=rerun2_df[['Unnamed: 0','session']]
sdMRI_rerun2_df.columns=['subject', 'session']

if save_file:
    sdMRI_rerun2_df.to_csv(sub_ses_rerun2_file, index=False, header=False)
sdMRI_rerun2_df<jupyter_output><empty_output><jupyter_text>Code below for testing:<jupyter_code># read bids layout and output all sessions
from bids import BIDSLayout
ppmi_layout=BIDSLayout(bids_dir)
print(ppmi_layout.get_sessions())
# select all the available T1w images and create the subject session list 
save_file=0

suffix    = 'T1w'
extension = 'nii.gz'
ppmi_file_list=ppmi_layout.get(suffix=suffix, extension=extension, return_type='file')
ppmi_file_names=[x.split('/')[-1] for x in ppmi_file_list]

subj_ses_df = pd.DataFrame({col_names[0]:[x.split('_')[0] for x in ppmi_file_names], col_names[1]: [x.split('_')[1].split('-')[-1] for x in ppmi_file_names]})
print(subj_ses_df)

# Generate subject,session file for fMRIPrep preporocessing 
if save_file:
    subj_ses_df.to_csv(ppmi_subj_ses_file, header=False, index=False)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1385790679908327, 'fraction_numerical': 0.017723453017570664, 'mean_word_length': 3.9177310293012773, 'pattern_counts': {'":': 4, '<': 26, '<?xml version=': 0, '>': 30, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '389381', 'license': 'MIT', 'n_tokens_mistral': 5545, 'n_tokens_neox': 5234, 'n_words': 794, 'path': 'workflow/fMRIPrep/fMRIPrep_help.ipynb', 'repo_name': 'Vincent-wq/mr_proc', 'size': '13090'}