<reponame>robbisg/mvpa_itab_wu
#!/usr/bin/env python
from bids import BIDSLayout
import os
import numpy as np
import pandas as pd
import click

@click.command()
@click.option("--path", default='.', help="Path where data is stored")
@click.option("--pipeline", default='afniglm')
@click.option("--task", default='rest')
@click.option("--hrf", default="SPM")
@click.option("--n_jobs", default=10)
@click.option("--polort", default=3)
@click.option("--add_confounds", is_flag=True)
@click.option("--event_label", default='trial_type')
@click.option("--extra_event", nargs=2, multiple=True)
def main(**args):
    
    path = "/mnt/DATA_4Tera/Dati_Sherlock/bids/"

    derivatives = os.path.join(path, "derivatives", "afniproc")
    print("mkdir -p {}".format(derivatives))
    os.system("mkdir -p {}".format(derivatives))

    subj_deriv = os.path.join(derivatives, 'sub-{subject}', "ses-{session}")

    layout = BIDSLayout(path)

    subjects = layout.get_subjects()
    sessions = layout.get_sessions()

    for session in sessions:
        for subj in subjects:

            deriv_dir = subj_deriv.format(session=session, subject=subj)
            print("mkdir -p {}".format(deriv_dir))
            os.system("mkdir -p {}".format(deriv_dir))
            
            # Create anat and func
            anat_dir = os.path.join(subj_deriv, "{datatype}").format(session=session, subject=subj, datatype='anat')
            func_dir = os.path.join(subj_deriv, "{datatype}").format(session=session, subject=subj, datatype='func')

            print("mkdir -p {}".format(anat_dir))
            os.system("mkdir -p {}".format(anat_dir))
            
            print("mkdir -p {}".format(func_dir))
            os.system("mkdir -p {}".format(func_dir))

            # Check and convert T1 to send Freesurfer segmentation 
            t1 = layout.get(subject=subj, session=session, suffix='T1w')[0]
            
            entities = t1.get_entities()
            pattern = os.path.join(subj_deriv, "{datatype}", "sub-{subject}[_ses-{session}][_desc-{desc}]_{suffix}.{extension}")
            entities['desc'] = 'fsprep'
            t1_fs = layout.build_path(entities, pattern, validate=False)
            
            entities['extension'] = 'txt'
            t1_log = layout.build_path(entities, pattern, validate=False)
            
            command = "check_dset_for_fs.py -input %s -fix_all -fix_out_prefix %s -fix_out_vox_dim 1 -verb > %s"
            command = command % (t1.path, t1_fs, t1_log)
            print(command)
            os.system(command)
            

            check_fs(layout, subj, session, subj_deriv)

            ### T1 to MNI space coreg
            entities = t1.get_entities()
            entities['desc'] = 'fsprep'
            entities['space'] = 'MNI152'

            pattern = os.path.join(subj_deriv, "{datatype}", "sub-{subject}[_ses-{session}][_space-{space}][_desc-{desc}]_{suffix}.{extension}")
            t1_mni = layout.build_path(entities, pattern, validate=False)
            
            command = "@auto_tlrc -base MNI152_2009_template.nii.gz -pad_base 35 -prefix {prefix} -input {input}"
            command = command.format(prefix=t1_mni, input=t1_fs)
            print(command)
            os.system(command)

            runs = layout.get_runs()
            ordered_bold = []
            for run in runs:
                fname = layout.get(session=session, subject=subj, run=run, suffix='bold')
                if len(fname) != 0:
                    ordered_bold.append(fname[0])

            # Slice time correction - motion correction - align EPI to Anat to MNI
            bold = layout.get(subject=subj, session=session, suffix='bold', extension='nii.gz')


            # 1D File
            slice_timing = np.array(bold[0].get_metadata()['SliceTiming'])
            slice_fname = os.path.join(path, "slice_timing.txt")
            np.savetxt(slice_fname, slice_timing, delimiter=' ', newline=' ', fmt='%.5f')

            epi = layout.get(subject=subj, session=session, suffix='bold', extension='nii.gz', run=1)[0].path
            child_epi = " ".join([b.path for b in ordered_bold])

            command = "align_epi_anat.py -anat {anat} -epi {epi} -child_epi {child_epi}"+\
                      " -epi_base 0 -tshift_opts -tpattern {tpattern} -epi2anat -giant_move"+\
                      " -tlrc_apar {tlrc_apar}"
            command = command.format(anat=t1_fs, epi=epi, child_epi=child_epi, tpattern=slice_fname, tlrc_apar=t1_mni)
            print(command)
            os.system(command)

            # Create mask
            automask = " ".join([b.filename[:-7]+'_tlrc_al+tlrc.HEAD' for b in bold])
            mean_mask_prefix = 'mean_mni.nii.gz'
            command = '3dTstat -prefix {prefix} {input}'.format(prefix=mean_mask_prefix, input=automask)
            print(command)
            os.system(command)

            entities = bold[0].get_entities()
            entities['suffix'] = 'mask'
            entities['space'] = 'MNI152'

            pattern = os.path.join(subj_deriv, "{datatype}", "sub-{subject}[_ses-{session}][_space-{space}][_desc-{desc}]_{suffix}.{extension}")
            mask_prefix = layout.build_path(entities, pattern, validate=False)
            command = '3dAutomask -prefix {prefix} {input}'.format(input=mean_mask_prefix, prefix=mask_prefix)
            print(command)
            os.system(command)

            # Clean files
            removed = " ".join([b.filename[:-7]+'_al+orig.*' for b in bold])
            command = "rm "+removed
            print(command)
            os.system(command)
            
            command = "rm __tt_*.*"
            print(command)
            os.system(command)

            command = "rm malldump.*"
            print(command)
            os.system(command)

            # Put files in BIDS
            header = ['trans_x', 'trans_y', 'trans_z', 'rot_x', 'rot_y', 'rot_z']
            motion_ordered = list()
            for bold in ordered_bold:

                motion_fname = bold.filename[:-7]+"_vr_motion.1D"
                motion = np.genfromtxt(motion_fname)

                motion_ordered.append(motion_fname)

                entities = bold.get_entities()
                entities['suffix'] = 'motion'
                entities['desc'] = 'volreg'
                entities['extension'] = 'tsv'
                
                pattern = os.path.join(subj_deriv, "{datatype}", "sub-{subject}_ses-{session}_task-{task}_run-{run:02d}_desc-{desc}_{suffix}.{extension}")
                motion_bids = pattern.format(**entities)

                motion = motion[:,[3,4,5,0,1,2]]
                np.savetxt(motion_bids, motion, fmt="%f", delimiter="\t", header="\t".join(header))

                print("rm "+motion_fname)
                os.system("rm "+motion_fname)


                for desc in ['mat', "reg_mat", "tlrc_mat"]:
                    affine_fname = bold.filename[:-7]+"_al_"+desc+".aff12.1D"
                    entities = bold.get_entities()
                    entities['suffix'] = 'affine'
                    entities['desc'] = desc.replace("_", "")
                    entities['extension'] = 'tsv'
                
                    pattern = os.path.join(subj_deriv, "{datatype}", "sub-{subject}_ses-{session}_task-{task}_run-{run:02d}_desc-{desc}_{suffix}.{extension}")
                    affine_bids = pattern.format(**entities)

                    command = "mv {0} {1}".format(affine_fname, affine_bids)
                    print(command)
                    os.system(command)
                    
               
                afni_bold = bold.filename[:-7]+"_tlrc_al+tlrc"
                entities = bold.get_entities()
                entities['desc'] = "afniproc"
                entities['extension'] = 'nii.gz'
                entities['space'] = 'MNI152'
                pattern = os.path.join(subj_deriv, "{datatype}", "sub-{subject}_ses-{session}_task-{task}_run-{run:02d}_space-{space}_desc-{desc}_{suffix}.{extension}")

                afni_bids = pattern.format(**entities)
                command = "3dcopy {0} {1}".format(afni_bold, afni_bids)
                print(command)
                os.system(command)

                print("rm {}*".format(afni_bold))
                os.system("rm {}*".format(afni_bold))


            # Create confound regressors
            # Motion
                 
            motion_files = list()
            for run in runs:
                f = layout.get(subject=subj, session=session, task=session, run=run, suffix='motion')
                if len(f) != 0:
                    motion_files.append(f[0])

            motion_df = [pd.read_csv(m.path, delimiter="\t") for m in motion_files]
            
            motion_demean = [m - m.mean(0) for m in motion_df]
            motion_demean = pd.concat(motion_demean)

            entities = motion_files[0].get_entities()
            entities['desc'] = 'demean'

            pattern = os.path.join(subj_deriv, "{datatype}", "sub-{subject}_ses-{session}_desc-{desc}_{suffix}.{extension}")
            demean_fname = pattern.format(**entities)
            motion_demean.to_csv(demean_fname, header=False, index=False, sep="\t")

            motion_deriv = [m.diff() for m in motion_df]
            motion_deriv = [m.fillna(0) for m in motion_deriv]
            motion_deriv = [m - m.mean(0) for m in motion_deriv]
            motion_deriv = pd.concat(motion_deriv)

            entities['desc'] = 'deriv'
            deriv_fname = pattern.format(**entities)
            motion_deriv.to_csv(deriv_fname, header=False, index=False, sep="\t")

            tr_counts = [m.shape[0] for m in motion_df]

            for j, t in enumerate(tr_counts):
                command = "1dBport -nodata {ntr} 1 -band 0.01 999 -invert -nozero > bpass.1D".format(ntr=t)
                print(command)
                os.system(command)

                command = "1d_tool.py -infile bpass.1D -pad_into_many_runs {run:1d} {n_runs}"+\
                          " -set_run_lengths {tr_counts} -write bpass.{run:02d}.1D"
                command = command.format(run=j+1, n_runs=str(len(tr_counts)), tr_counts=" ".join([str(t) for t in tr_counts]))
                print(command)
                os.system(command)

            entities['desc'] = 'bpass'
            entities['suffix'] = 'timeseries'
            entities['extension'] = '1D'
            bpass_fname = pattern.format(**entities)

            command = "1dcat bpass.*.1D > {}".format(bpass_fname)
            print(command)
            os.system(command)

            command = "rm bpass*"
            print(command)
            os.system(command) 


def check_fs(layout, subj, session, subj_deriv):

    # Check and convert T1 to send Freesurfer segmentation 
    t1 = layout.get(subject=subj, session=session, suffix='T1w')[0]
    
    entities = t1.get_entities()
    pattern = os.path.join(subj_deriv, "{datatype}", "sub-{subject}[_ses-{session}][_desc-{desc}]_{suffix}.{extension}")
    entities['desc'] = 'fsprep'
    t1_fs = layout.build_path(entities, pattern, validate=False)
    
    entities['extension'] = 'txt'
    t1_log = layout.build_path(entities, pattern, validate=False)
    
    command = "check_dset_for_fs.py -input %s -fix_all -fix_out_prefix %s -fix_out_vox_dim 1 -verb > %s"
    command = command % (t1.path, t1_fs, t1_log)
    print(command)
    os.system(command)




if __name__ == '__main__':
    main()
#!/usr/bin/env python
from bids import BIDSLayout
import os
import numpy as np
import pandas as pd
import click

@click.command()
@click.option("--path", default='.', help="Path where data is stored")
@click.option("--pipeline", default='afniglm')
@click.option("--task", default='rest')
@click.option("--hrf", default="SPM")
@click.option("--n_jobs", default=10)
@click.option("--polort", default=3)
@click.option("--add_confounds", is_flag=True)
@click.option("--event_label", default='trial_type')
@click.option("--extra_event", nargs=2, multiple=True)
def main(**args):
    
    path = "/mnt/DATA_4Tera/Dati_Sherlock/bids/"

    derivatives = os.path.join(path, "derivatives", "afniproc")
    print("mkdir -p {}".format(derivatives))
    os.system("mkdir -p {}".format(derivatives))

    subj_deriv = os.path.join(derivatives, 'sub-{subject}', "ses-{session}")

    layout = BIDSLayout(path)

    subjects = layout.get_subjects()
    sessions = layout.get_sessions()

    for session in sessions:
        for subj in subjects:

            deriv_dir = subj_deriv.format(session=session, subject=subj)
            print("mkdir -p {}".format(deriv_dir))
            os.system("mkdir -p {}".format(deriv_dir))
            
            # Create anat and func
            anat_dir = os.path.join(subj_deriv, "{datatype}").format(session=session, subject=subj, datatype='anat')
            func_dir = os.path.join(subj_deriv, "{datatype}").format(session=session, subject=subj, datatype='func')

            print("mkdir -p {}".format(anat_dir))
            os.system("mkdir -p {}".format(anat_dir))
            
            print("mkdir -p {}".format(func_dir))
            os.system("mkdir -p {}".format(func_dir))

            # Check and convert T1 to send Freesurfer segmentation 
            t1 = layout.get(subject=subj, session=session, suffix='T1w')[0]
            
            entities = t1.get_entities()
            pattern = os.path.join(subj_deriv, "{datatype}", "sub-{subject}[_ses-{session}][_desc-{desc}]_{suffix}.{extension}")
            entities['desc'] = 'fsprep'
            t1_fs = layout.build_path(entities, pattern, validate=False)
            
            entities['extension'] = 'txt'
            t1_log = layout.build_path(entities, pattern, validate=False)
            
            command = "check_dset_for_fs.py -input %s -fix_all -fix_out_prefix %s -fix_out_vox_dim 1 -verb > %s"
            command = command % (t1.path, t1_fs, t1_log)
            print(command)
            os.system(command)
            

            check_fs(layout, subj, session, subj_deriv)

            ### T1 to MNI space coreg
            entities = t1.get_entities()
            entities['desc'] = 'fsprep'
            entities['space'] = 'MNI152'

            pattern = os.path.join(subj_deriv, "{datatype}", "sub-{subject}[_ses-{session}][_space-{space}][_desc-{desc}]_{suffix}.{extension}")
            t1_mni = layout.build_path(entities, pattern, validate=False)
            
            command = "@auto_tlrc -base MNI152_2009_template.nii.gz -pad_base 35 -prefix {prefix} -input {input}"
            command = command.format(prefix=t1_mni, input=t1_fs)
            print(command)
            os.system(command)

            runs = layout.get_runs()
            ordered_bold = []
            for run in runs:
                fname = layout.get(session=session, subject=subj, run=run, suffix='bold')
                if len(fname) != 0:
                    ordered_bold.append(fname[0])

            # Slice time correction - motion correction - align EPI to Anat to MNI
            bold = layout.get(subject=subj, session=session, suffix='bold', extension='nii.gz')


            # 1D File
            slice_timing = np.array(bold[0].get_metadata()['SliceTiming'])
            slice_fname = os.path.join(path, "slice_timing.txt")
            np.savetxt(slice_fname, slice_timing, delimiter=' ', newline=' ', fmt='%.5f')

            epi = layout.get(subject=subj, session=session, suffix='bold', extension='nii.gz', run=1)[0].path
            child_epi = " ".join([b.path for b in ordered_bold])

            command = "align_epi_anat.py -anat {anat} -epi {epi} -child_epi {child_epi}"+\
                      " -epi_base 0 -tshift_opts -tpattern {tpattern} -epi2anat -giant_move"+\
                      " -tlrc_apar {tlrc_apar}"
            command = command.format(anat=t1_fs, epi=epi, child_epi=child_epi, tpattern=slice_fname, tlrc_apar=t1_mni)
            print(command)
            os.system(command)

            # Create mask
            automask = " ".join([b.filename[:-7]+'_tlrc_al+tlrc.HEAD' for b in bold])
            mean_mask_prefix = 'mean_mni.nii.gz'
            command = '3dTstat -prefix {prefix} {input}'.format(prefix=mean_mask_prefix, input=automask)
            print(command)
            os.system(command)

            entities = bold[0].get_entities()
            entities['suffix'] = 'mask'
            entities['space'] = 'MNI152'

            pattern = os.path.join(subj_deriv, "{datatype}", "sub-{subject}[_ses-{session}][_space-{space}][_desc-{desc}]_{suffix}.{extension}")
            mask_prefix = layout.build_path(entities, pattern, validate=False)
            command = '3dAutomask -prefix {prefix} {input}'.format(input=mean_mask_prefix, prefix=mask_prefix)
            print(command)
            os.system(command)

            # Clean files
            removed = " ".join([b.filename[:-7]+'_al+orig.*' for b in bold])
            command = "rm "+removed
            print(command)
            os.system(command)
            
            command = "rm __tt_*.*"
            print(command)
            os.system(command)

            command = "rm malldump.*"
            print(command)
            os.system(command)

            # Put files in BIDS
            header = ['trans_x', 'trans_y', 'trans_z', 'rot_x', 'rot_y', 'rot_z']
            motion_ordered = list()
            for bold in ordered_bold:

                motion_fname = bold.filename[:-7]+"_vr_motion.1D"
                motion = np.genfromtxt(motion_fname)

                motion_ordered.append(motion_fname)

                entities = bold.get_entities()
                entities['suffix'] = 'motion'
                entities['desc'] = 'volreg'
                entities['extension'] = 'tsv'
                
                pattern = os.path.join(subj_deriv, "{datatype}", "sub-{subject}_ses-{session}_task-{task}_run-{run:02d}_desc-{desc}_{suffix}.{extension}")
                motion_bids = pattern.format(**entities)

                motion = motion[:,[3,4,5,0,1,2]]
                np.savetxt(motion_bids, motion, fmt="%f", delimiter="\t", header="\t".join(header))

                print("rm "+motion_fname)
                os.system("rm "+motion_fname)


                for desc in ['mat', "reg_mat", "tlrc_mat"]:
                    affine_fname = bold.filename[:-7]+"_al_"+desc+".aff12.1D"
                    entities = bold.get_entities()
                    entities['suffix'] = 'affine'
                    entities['desc'] = desc.replace("_", "")
                    entities['extension'] = 'tsv'
                
                    pattern = os.path.join(subj_deriv, "{datatype}", "sub-{subject}_ses-{session}_task-{task}_run-{run:02d}_desc-{desc}_{suffix}.{extension}")
                    affine_bids = pattern.format(**entities)

                    command = "mv {0} {1}".format(affine_fname, affine_bids)
                    print(command)
                    os.system(command)
                    
               
                afni_bold = bold.filename[:-7]+"_tlrc_al+tlrc"
                entities = bold.get_entities()
                entities['desc'] = "afniproc"
                entities['extension'] = 'nii.gz'
                entities['space'] = 'MNI152'
                pattern = os.path.join(subj_deriv, "{datatype}", "sub-{subject}_ses-{session}_task-{task}_run-{run:02d}_space-{space}_desc-{desc}_{suffix}.{extension}")

                afni_bids = pattern.format(**entities)
                command = "3dcopy {0} {1}".format(afni_bold, afni_bids)
                print(command)
                os.system(command)

                print("rm {}*".format(afni_bold))
                os.system("rm {}*".format(afni_bold))


            # Create confound regressors
            # Motion
                 
            motion_files = list()
            for run in runs:
                f = layout.get(subject=subj, session=session, task=session, run=run, suffix='motion')
                if len(f) != 0:
                    motion_files.append(f[0])

            motion_df = [pd.read_csv(m.path, delimiter="\t") for m in motion_files]
            
            motion_demean = [m - m.mean(0) for m in motion_df]
            motion_demean = pd.concat(motion_demean)

            entities = motion_files[0].get_entities()
            entities['desc'] = 'demean'

            pattern = os.path.join(subj_deriv, "{datatype}", "sub-{subject}_ses-{session}_desc-{desc}_{suffix}.{extension}")
            demean_fname = pattern.format(**entities)
            motion_demean.to_csv(demean_fname, header=False, index=False, sep="\t")

            motion_deriv = [m.diff() for m in motion_df]
            motion_deriv = [m.fillna(0) for m in motion_deriv]
            motion_deriv = [m - m.mean(0) for m in motion_deriv]
            motion_deriv = pd.concat(motion_deriv)

            entities['desc'] = 'deriv'
            deriv_fname = pattern.format(**entities)
            motion_deriv.to_csv(deriv_fname, header=False, index=False, sep="\t")

            tr_counts = [m.shape[0] for m in motion_df]

            for j, t in enumerate(tr_counts):
                command = "1dBport -nodata {ntr} 1 -band 0.01 999 -invert -nozero > bpass.1D".format(ntr=t)
                print(command)
                os.system(command)

                command = "1d_tool.py -infile bpass.1D -pad_into_many_runs {run:1d} {n_runs}"+\
                          " -set_run_lengths {tr_counts} -write bpass.{run:02d}.1D"
                command = command.format(run=j+1, n_runs=str(len(tr_counts)), tr_counts=" ".join([str(t) for t in tr_counts]))
                print(command)
                os.system(command)

            entities['desc'] = 'bpass'
            entities['suffix'] = 'timeseries'
            entities['extension'] = '1D'
            bpass_fname = pattern.format(**entities)

            command = "1dcat bpass.*.1D > {}".format(bpass_fname)
            print(command)
            os.system(command)

            command = "rm bpass*"
            print(command)
            os.system(command) 


def check_fs(layout, subj, session, subj_deriv):

    # Check and convert T1 to send Freesurfer segmentation 
    t1 = layout.get(subject=subj, session=session, suffix='T1w')[0]
    
    entities = t1.get_entities()
    pattern = os.path.join(subj_deriv, "{datatype}", "sub-{subject}[_ses-{session}][_desc-{desc}]_{suffix}.{extension}")
    entities['desc'] = 'fsprep'
    t1_fs = layout.build_path(entities, pattern, validate=False)
    
    entities['extension'] = 'txt'
    t1_log = layout.build_path(entities, pattern, validate=False)
    
    command = "check_dset_for_fs.py -input %s -fix_all -fix_out_prefix %s -fix_out_vox_dim 1 -verb > %s"
    command = command % (t1.path, t1_fs, t1_log)
    print(command)
    os.system(command)




if __name__ == '__main__':
    main()
#!/usr/bin/env python # Create anat and func # Check and convert T1 to send Freesurfer segmentation ### T1 to MNI space coreg # Slice time correction - motion correction - align EPI to Anat to MNI # 1D File # Create mask # Clean files # Put files in BIDS # Create confound regressors # Motion # Check and convert T1 to send Freesurfer segmentation