diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..7e99e36 --- /dev/null +++ b/.gitignore @@ -0,0 +1 @@ +*.pyc \ No newline at end of file diff --git a/ICA_AROMA.py b/ICA_AROMA.py index 76e31b7..a769fac 100644 --- a/ICA_AROMA.py +++ b/ICA_AROMA.py @@ -10,6 +10,7 @@ import subprocess import ICA_AROMA_functions as aromafunc import shutil +import classification_plots # Change to script directory cwd = os.path.realpath(os.path.curdir) @@ -22,26 +23,26 @@ # Required options reqoptions = parser.add_argument_group('Required arguments') -reqoptions.add_argument('-o', '-out', dest="outDir", required=True, help='Output directory name' ) +reqoptions.add_argument('-o', '-out', dest="outDir", required=True, help='Output directory name') # Required options in non-Feat mode nonfeatoptions = parser.add_argument_group('Required arguments - generic mode') -nonfeatoptions.add_argument('-i', '-in',dest="inFile", required=False, help='Input file name of fMRI data (.nii.gz)') +nonfeatoptions.add_argument('-i', '-in', dest="inFile", required=False, help='Input file name of fMRI data (.nii.gz)') nonfeatoptions.add_argument('-mc', dest="mc", required=False, help='File name of the motion parameters obtained after motion realingment (e.g., FSL mcflirt). Note that the order of parameters does not matter, should your file not originate from FSL mcflirt. (e.g., /home/user/PROJECT/SUBJECT.feat/mc/prefiltered_func_data_mcf.par') -nonfeatoptions.add_argument('-a','-affmat', dest="affmat", default="", help='File name of the mat-file describing the affine registration (e.g., FSL FLIRT) of the functional data to structural space (.mat file). (e.g., /home/user/PROJECT/SUBJECT.feat/reg/example_func2highres.mat') -nonfeatoptions.add_argument('-w','-warp', dest="warp", default="", help='File name of the warp-file describing the non-linear registration (e.g., FSL FNIRT) of the structural data to MNI152 space (.nii.gz). (e.g., /home/user/PROJECT/SUBJECT.feat/reg/highres2standard_warp.nii.gz') -nonfeatoptions.add_argument('-m','-mask', dest="mask", default="", help='File name of the mask to be used for MELODIC (denoising will be performed on the original/non-masked input data)') +nonfeatoptions.add_argument('-a', '-affmat', dest="affmat", default="", help='File name of the mat-file describing the affine registration (e.g., FSL FLIRT) of the functional data to structural space (.mat file). (e.g., /home/user/PROJECT/SUBJECT.feat/reg/example_func2highres.mat') +nonfeatoptions.add_argument('-w', '-warp', dest="warp", default="", help='File name of the warp-file describing the non-linear registration (e.g., FSL FNIRT) of the structural data to MNI152 space (.nii.gz). (e.g., /home/user/PROJECT/SUBJECT.feat/reg/highres2standard_warp.nii.gz') +nonfeatoptions.add_argument('-m', '-mask', dest="mask", default="", help='File name of the mask to be used for MELODIC (denoising will be performed on the original/non-masked input data)') # Required options in Feat mode featoptions = parser.add_argument_group('Required arguments - FEAT mode') -featoptions.add_argument('-f', '-feat',dest="inFeat", required=False, help='Feat directory name (Feat should have been run without temporal filtering and including registration to MNI152)') +featoptions.add_argument('-f', '-feat', dest="inFeat", required=False, help='Feat directory name (Feat should have been run without temporal filtering and including registration to MNI152)') # Optional options optoptions = parser.add_argument_group('Optional arguments') -optoptions.add_argument('-tr', dest="TR", help='TR in seconds',type=float) +optoptions.add_argument('-tr', dest="TR", help='TR in seconds', type=float) optoptions.add_argument('-den', dest="denType", default="nonaggr", help='Type of denoising strategy: \'no\': only classification, no denoising; \'nonaggr\': non-aggresssive denoising (default); \'aggr\': aggressive denoising; \'both\': both aggressive and non-aggressive denoising (seperately)') -optoptions.add_argument('-md','-meldir', dest="melDir", default="",help='MELODIC directory name, in case MELODIC has been run previously.') -optoptions.add_argument('-dim', dest="dim", default=0,help='Dimensionality reduction into #num dimensions when running MELODIC (default: automatic estimation; i.e. -dim 0)',type=int) +optoptions.add_argument('-md', '-meldir', dest="melDir", default="",help='MELODIC directory name, in case MELODIC has been run previously.') +optoptions.add_argument('-dim', dest="dim", default=0, help='Dimensionality reduction into #num dimensions when running MELODIC (default: automatic estimation; i.e. -dim 0)', type=int) print('\n------------------------------- RUNNING ICA-AROMA ------------------------------- ') print('--------------- \'ICA-based Automatic Removal Of Motion Artifacts\' --------------- \n') @@ -51,69 +52,70 @@ args = parser.parse_args() # Define variables based on the type of input (i.e. Feat directory or specific input arguments), and check whether the specified files exist. -cancel=False +cancel = False + if args.inFeat: - inFeat = args.inFeat - - # Check whether the Feat directory exists - if not os.path.isdir(inFeat): - print('The specified Feat directory does not exist.') - print('\n----------------------------- ICA-AROMA IS CANCELED -----------------------------\n') - exit() - - # Define the variables which should be located in the Feat directory - inFile = os.path.join(args.inFeat,'filtered_func_data.nii.gz') - mc = os.path.join(args.inFeat,'mc','prefiltered_func_data_mcf.par') - affmat = os.path.join(args.inFeat,'reg','example_func2highres.mat') - warp = os.path.join(args.inFeat,'reg','highres2standard_warp.nii.gz') - - # Check whether these files actually exist - if not os.path.isfile(inFile): - print('Missing filtered_func_data.nii.gz in Feat directory.') - cancel=True - if not os.path.isfile(mc): - print('Missing mc/prefiltered_func_data_mcf.mat in Feat directory.') - cancel=True - if not os.path.isfile(affmat): - print('Missing reg/example_func2highres.mat in Feat directory.') - cancel=True - if not os.path.isfile(warp): - print('Missing reg/highres2standard_warp.nii.gz in Feat directory.') - cancel=True - - # Check whether a melodic.ica directory exists - if os.path.isdir(os.path.join(args.inFeat,'filtered_func_data.ica')): - melDir = os.path.join(args.inFeat,'filtered_func_data.ica') - else: - melDir=args.melDir + inFeat = args.inFeat + + # Check whether the Feat directory exists + if not os.path.isdir(inFeat): + print('The specified Feat directory does not exist.') + print('\n----------------------------- ICA-AROMA IS CANCELED -----------------------------\n') + exit() + + # Define the variables which should be located in the Feat directory + inFile = os.path.join(args.inFeat, 'filtered_func_data.nii.gz') + mc = os.path.join(args.inFeat, 'mc', 'prefiltered_func_data_mcf.par') + affmat = os.path.join(args.inFeat, 'reg', 'example_func2highres.mat') + warp = os.path.join(args.inFeat, 'reg', 'highres2standard_warp.nii.gz') + + # Check whether these files actually exist + if not os.path.isfile(inFile): + print('Missing filtered_func_data.nii.gz in Feat directory.') + cancel = True + if not os.path.isfile(mc): + print('Missing mc/prefiltered_func_data_mcf.mat in Feat directory.') + cancel = True + if not os.path.isfile(affmat): + print('Missing reg/example_func2highres.mat in Feat directory.') + cancel = True + if not os.path.isfile(warp): + print('Missing reg/highres2standard_warp.nii.gz in Feat directory.') + cancel = True + + # Check whether a melodic.ica directory exists + if os.path.isdir(os.path.join(args.inFeat, 'filtered_func_data.ica')): + melDir = os.path.join(args.inFeat, 'filtered_func_data.ica') + else: + melDir = args.melDir else: - inFile = args.inFile - mc = args.mc - affmat = args.affmat - warp = args.warp - melDir = args.melDir - - # Check whether the files exist - if not inFile: - print('No input file specified.') - else: - if not os.path.isfile(inFile): - print('The specified input file does not exist.') - cancel=True - if not mc: - print('No mc file specified.') - else: - if not os.path.isfile(mc): - print('The specified mc file does does not exist.') - cancel=True - if affmat: - if not os.path.isfile(affmat): - print('The specified affmat file does not exist.') - cancel=True - if warp: - if not os.path.isfile(warp): - print('The specified warp file does not exist.') - cancel=True + inFile = args.inFile + mc = args.mc + affmat = args.affmat + warp = args.warp + melDir = args.melDir + + # Check whether the files exist + if not inFile: + print('No input file specified.') + else: + if not os.path.isfile(inFile): + print('The specified input file does not exist.') + cancel = True + if not mc: + print('No mc file specified.') + else: + if not os.path.isfile(mc): + print('The specified mc file does does not exist.') + cancel = True + if affmat: + if not os.path.isfile(affmat): + print('The specified affmat file does not exist.') + cancel = True + if warp: + if not os.path.isfile(warp): + print('The specified warp file does not exist.') + cancel = True # Parse the arguments which do not depend on whether a Feat directory has been specified outDir = args.outDir @@ -122,67 +124,67 @@ # Check if the mask exists, when specified. if args.mask: - if not os.path.isfile(args.mask): - print('The specified mask does not exist.') - cancel=True + if not os.path.isfile(args.mask): + print('The specified mask does not exist.') + cancel = True # Check if the type of denoising is correctly specified, when specified if not (denType == 'nonaggr') and not (denType == 'aggr') and not (denType == 'both') and not (denType == 'no'): - print('Type of denoising was not correctly specified. Non-aggressive denoising will be run.') - denType='nonaggr' + print('Type of denoising was not correctly specified. Non-aggressive denoising will be run.') + denType = 'nonaggr' # If the criteria for file/directory specifications have not been met. Cancel ICA-AROMA. if cancel: - print('\n----------------------------- ICA-AROMA IS CANCELED -----------------------------\n') - exit() + print('\n----------------------------- ICA-AROMA IS CANCELED -----------------------------\n') + exit() #------------------------------------------- PREPARE -------------------------------------------# # Define the FSL-bin directory -fslDir = os.path.join(os.environ["FSLDIR"],'bin','') +fslDir = os.path.join(os.environ["FSLDIR"], 'bin', '') # Create output directory if needed if not os.path.isdir(outDir): - os.makedirs(outDir) + os.makedirs(outDir) # Get TR of the fMRI data, if not specified if args.TR: - TR = args.TR + TR = args.TR else: - cmd = ' '.join([os.path.join(fslDir,'fslinfo'), - inFile, - '| grep pixdim4 | awk \'{print $2}\'']) - TR=float(subprocess.getoutput(cmd)) + cmd = ' '.join([os.path.join(fslDir, 'fslinfo'), + inFile, + '| grep pixdim4 | awk \'{print $2}\'']) + TR = float(subprocess.getoutput(cmd)) # Check TR if TR == 1: - print('Warning! Please check whether the determined TR (of ' + str(TR) + 's) is correct!\n') + print('Warning! Please check whether the determined TR (of ' + str(TR) + 's) is correct!\n') elif TR == 0: - print('TR is zero. ICA-AROMA requires a valid TR and will therefore exit. Please check the header, or define the TR as an additional argument.\n----------------------------- ICA-AROMA IS CANCELED -----------------------------\n') + print('TR is zero. ICA-AROMA requires a valid TR and will therefore exit. Please check the header, or define the TR as an additional argument.\n----------------------------- ICA-AROMA IS CANCELED -----------------------------\n') # Define/create mask. Either by making a copy of the specified mask, or by creating a new one. -mask = os.path.join(outDir,'mask.nii.gz') +mask = os.path.join(outDir, 'mask.nii.gz') if args.mask: - shutil.copyfile(args.mask, mask) + shutil.copyfile(args.mask, mask) else: - # If a Feat directory is specified, and an example_func is present use example_func to create a mask - if args.inFeat and os.path.isfile(os.path.join(inFeat,'example_func.nii.gz')): - os.system(' '.join([os.path.join(fslDir,'bet'), - os.path.join(inFeat,'example_func.nii.gz'), - os.path.join(outDir,'bet'), - '-f 0.3 -n -m -R'])) - os.system(' '.join(['mv', - os.path.join(outDir,'bet_mask.nii.gz'), - mask])) - if os.path.isfile(os.path.join(outDir,'bet.nii.gz')): - os.remove(os.path.join(outDir,'bet.nii.gz')) - else: - if args.inFeat: - print(' - No example_func was found in the Feat directory. A mask will be created including all voxels with varying intensity over time in the fMRI data. Please check!\n') - os.system(' '.join([os.path.join(fslDir,'fslmaths'), - inFile, - '-Tstd -bin', - mask])) + # If a Feat directory is specified, and an example_func is present use example_func to create a mask + if args.inFeat and os.path.isfile(os.path.join(inFeat, 'example_func.nii.gz')): + os.system(' '.join([os.path.join(fslDir, 'bet'), + os.path.join(inFeat, 'example_func.nii.gz'), + os.path.join(outDir, 'bet'), + '-f 0.3 -n -m -R'])) + os.system(' '.join(['mv', + os.path.join(outDir, 'bet_mask.nii.gz'), + mask])) + if os.path.isfile(os.path.join(outDir, 'bet.nii.gz')): + os.remove(os.path.join(outDir, 'bet.nii.gz')) + else: + if args.inFeat: + print(' - No example_func was found in the Feat directory. A mask will be created including all voxels with varying intensity over time in the fMRI data. Please check!\n') + os.system(' '.join([os.path.join(fslDir, 'fslmaths'), + inFile, + '-Tstd -bin', + mask])) #---------------------------------------- Run ICA-AROMA ----------------------------------------# @@ -192,27 +194,30 @@ print('Step 2) Automatic classification of the components') print(' - registering the spatial maps to MNI') -melIC = os.path.join(outDir,'melodic_IC_thr.nii.gz') -melIC_MNI = os.path.join(outDir,'melodic_IC_thr_MNI2mm.nii.gz') +melIC = os.path.join(outDir, 'melodic_IC_thr.nii.gz') +melIC_MNI = os.path.join(outDir, 'melodic_IC_thr_MNI2mm.nii.gz') aromafunc.register2MNI(fslDir, melIC, melIC_MNI, affmat, warp) print(' - extracting the CSF & Edge fraction features') edgeFract, csfFract = aromafunc.feature_spatial(fslDir, outDir, scriptDir, melIC_MNI) print(' - extracting the Maximum RP correlation feature') -melmix = os.path.join(outDir,'melodic.ica','melodic_mix') +melmix = os.path.join(outDir, 'melodic.ica', 'melodic_mix') maxRPcorr = aromafunc.feature_time_series(melmix, mc) print(' - extracting the High-frequency content feature') -melFTmix = os.path.join(outDir,'melodic.ica','melodic_FTmix') +melFTmix = os.path.join(outDir, 'melodic.ica', 'melodic_FTmix') HFC = aromafunc.feature_frequency(melFTmix, TR) print(' - classification') motionICs = aromafunc.classification(outDir, maxRPcorr, edgeFract, HFC, csfFract) +classification_plots.classification_plot(os.path.join(outDir, 'classification_overview.txt'), + outDir) + if (denType != 'no'): - print('Step 3) Data denoising') - aromafunc.denoising(fslDir, inFile, outDir, melmix, denType, motionICs) + print('Step 3) Data denoising') + aromafunc.denoising(fslDir, inFile, outDir, melmix, denType, motionICs) # Remove thresholded melodic_IC file os.remove(melIC) diff --git a/ICA_AROMA_functions.py b/ICA_AROMA_functions.py index 8584966..8ba20fe 100644 --- a/ICA_AROMA_functions.py +++ b/ICA_AROMA_functions.py @@ -58,7 +58,8 @@ def runICA(fslDir, inFile, outDir, melDirIn, mask, dim, TR): # Create symbolic links to the items in the specified melodic directory os.makedirs(melDir) for item in os.listdir(melDirIn): - os.symlink(os.path.join(melDirIn, item), os.path.join(melDir, item)) + os.symlink(os.path.join(melDirIn, item), + os.path.join(melDir, item)) # Run mixture modeling os.system(' '.join([os.path.join(fslDir, 'melodic'), @@ -508,14 +509,25 @@ def classification(outDir, maxRPcorr, edgeFract, HFC, csfFract): # Create a summary overview of the classification txt = open(os.path.join(outDir, 'classification_overview.txt'), 'w') - txt.write('IC' + '\t' + 'Motion/noise' + '\t' + 'maximum RP correlation' + '\t' + 'Edge-fraction' + '\t\t' + 'High-frequency content' + '\t' + 'CSF-fraction') + txt.write('\t'.join(['IC', + 'Motion/noise', + 'maximum RP correlation', + 'Edge-fraction', + 'High-frequency content', + 'CSF-fraction'])) txt.write('\n') for i in range(0, len(csfFract)): if (proj[i] > 0) or (csfFract[i] > thr_csf) or (HFC[i] > thr_HFC): classif = "True" else: classif = "False" - txt.write('%.0f\t%s\t\t%.2f\t\t\t%.2f\t\t\t%.2f\t\t\t%.2f\n' % (i + 1, classif, maxRPcorr[i], edgeFract[i], HFC[i], csfFract[i])) + txt.write('\t'.join([i + 1, + classif, + maxRPcorr[i], + edgeFract[i], + HFC[i], + csfFract[i]])) + txt.write('\n') txt.close() return motionICs diff --git a/classification_plots.py b/classification_plots.py new file mode 100644 index 0000000..61e2c4c --- /dev/null +++ b/classification_plots.py @@ -0,0 +1,243 @@ +from __future__ import print_function + + +def classification_plot(myinput, outDir): + + import pandas as pd + import numpy as np + import seaborn as sns + import matplotlib.pyplot as plt + from matplotlib import gridspec + import glob + import os + + ###---Start---### + # find files + myfiles = glob.glob(myinput) + print('Found', len(myfiles), 'file(s)') + + # load in data from files + count = 0 + for m in myfiles: + + res = [] + + tmp = open(m, 'r').read().split('\n') + + for t in tmp[1:-1]: + vals = t.split('\t') + res.append([vals[1], + float(vals[2]), + float(vals[3]), + float(vals[4]), + float(vals[5])]) + + if count == 0: + df = pd.DataFrame.from_records(res) + else: + df2 = pd.DataFrame.from_records(res) + df = df.append(df2, ignore_index=True) + + count += 1 + + # get counts + ncomp = len(df) + nmot = len(df.loc[df[0] == "True"]) + print('Found', nmot, 'head motion-related components in a total of', ncomp, 'components.') + + # add dummy components if needed, this is just for making the plots look nice + tmp = df.loc[df[0] == "True"] + if len(tmp) == 0: + df3 = pd.DataFrame.from_records([["True", 1., 1., 0., 0.], + ["True", 1., 1., 0., 0.], + ["True", 1., 1., 0., 0.]]) + df = df.append(df3, ignore_index=True) + tmp = df.loc[df[0] == "False"] + if len(tmp) == 0: + df3 = pd.DataFrame.from_records([["False", 0., 0., 0., 0.], + ["False", 0., 0., 0., 0.], + ["False", 0., 0., 0., 0.]]) + df = df.append(df3, ignore_index=True) + + # rename columns + df = df.rename(index=str, columns={0: 'Motion', + 1: 'RP', + 2: 'Edge', + 3: 'Freq', + 4: 'CSF'}) + + # Make pretty figure + # styling + sns.set_style('white') + colortrue = "#FFBF17" + colorfalse = "#69A00A" + + # create figure + fig = plt.figure(figsize=[12, 4]) + + # define grids + gs = gridspec.GridSpec(4, 7, wspace=1) + gs00 = gridspec.GridSpecFromSubplotSpec(4, 4, subplot_spec=gs[:, 0:3]) + gs01 = gridspec.GridSpecFromSubplotSpec(4, 1, subplot_spec=gs[:, 3:5]) + gs02 = gridspec.GridSpecFromSubplotSpec(4, 1, subplot_spec=gs[:, 5:7]) + + # define subplots + # Edge/RP + ax1 = fig.add_subplot(gs00[1:4, 0:3]) + # distribution edge (ax1 top) + ax1t = fig.add_subplot(gs00[0, 0:3]) + # distribution RP (ax1 right) + ax1r = fig.add_subplot(gs00[1:4, 3]) + # Freq + ax2 = fig.add_subplot(gs01[1:4, :]) + # CSF + ax3 = fig.add_subplot(gs02[1:4, :]) + + # plot Freq + sns.boxplot(x="Motion", + y="Freq", + data=df, + ax=ax2, + palette=[colortrue, colorfalse]) + ax2.hlines(0.35, -1, 2, zorder=0, linestyles='dotted', linewidth=0.5) + ax2.set_ylim([0, 1]) + ax2.set_xlabel('Classification', fontsize=14, labelpad=10) + ax2.set_ylabel('High-Frequency Content', fontsize=14) + ax2.set_xticklabels(['Motion', 'Other']) + ax2.tick_params(axis='both', labelsize=12) + sns.despine(ax=ax2) + + # plot CSF + sns.boxplot(x="Motion", + y="CSF", + data=df, + ax=ax3, + palette=[colortrue, colorfalse]) + ax3.hlines(0.1, -1, 2, zorder=0, linestyles='dotted', linewidth=0.5) + ax3.set_ylim([0, 1]) + ax3.set_xlabel('Classification', fontsize=14, labelpad=10) + ax3.set_ylabel('CSF Fraction', fontsize=14) + ax3.set_xticklabels(['Motion', 'Other']) + ax3.tick_params(axis='both', labelsize=12) + sns.despine(ax=ax3) + + # plot Edge/RP relationship + # obtain projection line + hyp = np.array([-19.9751070082159, 9.95127547670627, 24.8333160239175]) + a = -hyp[1] / hyp[2] + xx = np.linspace(0, 1) + yy = a * xx - hyp[0] / hyp[2] + # plot scatter and line + if len(df) > 100: + sizemarker = 6 + else: + sizemarker = 10 + ax1.scatter(x="RP", + y="Edge", + data=df.loc[df['Motion'] == "False"], + color=colorfalse, + s=sizemarker) + # plot true ones on top to see how much the go over the border + # this gives an indication for how many were selected using the + # two other features + ax1.scatter(x="RP", + y="Edge", + data=df.loc[df['Motion'] == "True"], + color=colortrue, + s=sizemarker) + # add decision boundary + ax1.plot(xx, yy, '.', color="k", markersize=1) + # styling + ax1.set_ylim([0, 1]) + ax1.set_xlim([0, 1]) + ax1.set_xlabel('Maximum RP Correlation', fontsize=14, labelpad=10) + ax1.set_ylabel('Edge Fraction', fontsize=14) + ax1.set_xticks(np.arange(0, 1.2, 0.2)) + ax1.set_yticks(np.arange(0, 1.2, 0.2)) + ax1.tick_params(axis='both', labelsize=12) + + # plot distributions + # RP + sns.distplot(df.loc[df['Motion'] == "True", "RP"], + ax=ax1t, + color=colortrue, + hist_kws={'alpha': 0.2}) + sns.distplot(df.loc[df['Motion'] == "False", "RP"], + ax=ax1t, + color=colorfalse, + hist_kws={'alpha': 0.2}) + ax1t.set_xlim([0, 1]) + + # Edge + sns.distplot(df.loc[df['Motion'] == "True", "Edge"], + ax=ax1r, + vertical=True, + color=colortrue, + hist_kws={'alpha': 0.2}) + sns.distplot(df.loc[df['Motion'] == "False", "Edge"], + ax=ax1r, + vertical=True, + color=colorfalse, + hist_kws={'alpha': 0.2}) + ax1r.set_ylim([0, 1]) + + # cosmetics + for myax in [ax1t, ax1r]: + myax.set_xticks([]) + myax.set_yticks([]) + myax.set_xlabel('') + myax.set_ylabel('') + myax.spines['right'].set_visible(False) + myax.spines['top'].set_visible(False) + myax.spines['bottom'].set_visible(False) + myax.spines['left'].set_visible(False) + + # bring tickmarks back + for myax in fig.get_axes(): + myax.tick_params(which="major", direction='in', length=3) + + # add figure title + plt.suptitle('Component Assessment', fontsize=20) + + # outtakes + plt.savefig(os.path.join(outDir, 'ICA_AROMA_component_assessment.pdf'), + bbox_inches='tight') + + return + + +# allow use of module on its own +if __name__ == '__main__': + + import argparse + + parser = argparse.ArgumentParser(description="""Plot component classification overview + similar to plot in the main AROMA paper""") + # Required options + reqoptions = parser.add_argument_group('Required arguments') + reqoptions.add_argument('-i', '-in', + dest='myinput', + required=True, + help="""Input query or filename. + Use quotes when specifying a query""") + + optoptions = parser.add_argument_group('Optional arguments') + optoptions.add_argument('-outdir', + dest='outDir', + required=False, + default='.', + help="""Specification of directory + where figure will be saved""") + optoptions.add_argument('-type', + dest='plottype', + required=False, + default='assessment', + help="""Specification of the type of plot you want. + Currently this is a placeholder option for + potential other plots that might be added + in the future.""") + # parse arguments + args = parser.parse_args() + + if args.plottype == 'assessment': + classification_plot(args.myinput, args.outDir)