statistics-surface not compatible with Matlab R2019b
Feedback from a Clinica user: 

When running `statistics-surface` with Matlab R2019b, the following error is raised:
```

Standard error:
MATLAB code threw an exception:
Error: File: /Users/alexandre.routier/Projects/Git/clinica_dep/clinica/lib/clinicasurfstat/clinicasurfstat.m Line: 158 Column: 23
Indexing with parentheses '()' must appear as the last operation of a valid indexing expression.
File:/private/var/folders/m_/j76n37kn4vs6zsj8fq0qcgy0000dn9/T/tmpga40_rup/statistics-surface/surfstat/pyscript.m
Name:pyscript
Line:6
Return code: 0
```

The line involved in `pyscript.m` is:
```matlab
            factor1 = char(term(tsvdata{indexunique})){1};
```

It turns out that using older version of Matlab (when I tested with R2018b), the pipeline runs without issue.

diff --git a/clinica/lib/clinicasurfstat/clinicasurfstat.m b/clinica/lib/clinicasurfstat/clinicasurfstat.m
index 63b723c0..e490bb45 100755
--- a/clinica/lib/clinicasurfstat/clinicasurfstat.m
+++ b/clinica/lib/clinicasurfstat/clinicasurfstat.m
@@ -42,7 +42,7 @@ while length(pp) >= 2,
     prop = pp{1};
     val = pp{2};
     pp = pp(3 : end);
-    if ischar(prop) 
+    if ischar(prop)
         switch prop
             case 'sizeoffwhm'
                 sizeoffwhm = val;
@@ -108,7 +108,7 @@ for indexsubject = 1 : nrsubject
         error('Unexpected dimension of Y in SurfStatReadData')
     end
     disp(['The subject ID is: ', subjectname, '_', sessionname] )
-    if indexsubject == 1 % check the header format of the input tsv 
+    if indexsubject == 1 % check the header format of the input tsv
         thicksubject = zeros(nrsubject, size(Y,2));
         %if startsWith(contrast, '-') this works for matlab2016b
         if strfind(contrast, '-') % check for negative correlation analysis
@@ -155,8 +155,9 @@ switch glmtype
         if with_intercation ~= 1 % case without interaction, contrast should be negative and positive
             contrastpos    = eval([contrast '(1)']) - eval([ contrast '(2)']); % use char(eval(contrast))
             contrasteffectgroupneg    = eval([contrast '(2)']) - eval([ contrast '(1)']); % use char(eval(contrast))
-            factor1 = char(term(tsvdata{indexunique})){1};
-            factor2 = char(term(tsvdata{indexunique})){2};
+            cell_factor = char(term(tsvdata{indexunique}));
+            factor1 = cell_factor{1};
+            factor2 = cell_factor{2};
 
             thicksubject = thicksubject';
             slmmodel = SurfStatLinMod(thicksubject, eval(designmatrix), averagesurface);
@@ -287,7 +288,7 @@ switch glmtype
                 categorical_contrast = strs_contrast{2};
             else
                 continue_contrast = strs_contrast{2};
-                categorical_contrast = strs_contrast{1} ;               
+                categorical_contrast = strs_contrast{1} ;
             end
             categorical_term = eval(categorical_contrast);
             continue_term = eval(continue_contrast);
@@ -309,14 +310,14 @@ switch glmtype
             ylabel('Yseed')
             set(gcf,'PaperPositionMode','auto');
             print(['Highest T value vertex Yseed versus ' continue_contrast] ,'-djpeg','-r0'); close
-            
+
             % F statistics are obtained by comparing nested models
             % test if the contrast has the same format like in designmatrix
             if strfind(designmatrix, contrast) ~= 0 % contrast is inside designmatrix with the same order
                 designmatrix_reduced_model = strrep(strrep(designmatrix, ' ', ''), strcat('+',contrast), '');
             else % if the order is not the same in designmatrix
                 contrast_inverse = strsplit(contrast, '*');
-                contrast = char(strcat(contrast_inverse(2), '*', contrast_inverse(1))); 
+                contrast = char(strcat(contrast_inverse(2), '*', contrast_inverse(1)));
                 designmatrix_reduced_model = strrep(strrep(designmatrix, ' ', ''), strcat('+',contrast), '');
             end
             slm_reduce = SurfStatLinMod( thicksubject, eval(char(designmatrix_reduced_model)), averagesurface );
@@ -326,8 +327,8 @@ switch glmtype
             print(['interaction-' contrast '_measure-' measurelabel '_fwhm-' num2str(sizeoffwhm) '_FStatistics'] ,'-djpeg','-r0'); close
             fvaluewithmask = slm.t .* mask;
             save(['interaction-' contrast '_measure-' measurelabel '_fwhm-' num2str(sizeoffwhm) '_FStatistics.mat'],'fvaluewithmask');
-            
-            % pvalue 
+
+            % pvalue
             [ pval, ~, clus ] = SurfStatP( slm , mask, clusterthreshold);
             pval.thresh = thresholdcorrectedpvalue;
             SurfStatView( pval, averagesurface, ['Corrected Pvalue for interaction ', contrast, ' (clusterthreshold = ' num2str(clusterthreshold) ')']);
@@ -337,7 +338,7 @@ switch glmtype
             print(['interaction-' contrast '_measure-' measurelabel '_fwhm-' num2str(sizeoffwhm) '_correctedPValue'] ,'-djpeg','-r0'); close
             correctedpvaluesstruct = pval;
             save(['interaction-' contrast '_measure-' measurelabel '_fwhm-' num2str(sizeoffwhm) '_correctedPValue.mat'],'correctedpvaluesstruct');
-            
+
             % display the details of cluster-wise correction
             disp('###')
             disp('After correction(Clusterwise Correction for Multiple Comparisons): ')
@@ -348,13 +349,13 @@ switch glmtype
                 disp('No cluster found!')
             end
             disp('###')
-            
+
             % FDR correction with 0.05 threshold, q value
             qval = SurfStatQ( slm , mask );
             SurfStatView( qval, averagesurface, ['False discovery rate for interaction ' contrast]);
             set(gcf,'PaperPositionMode','auto');
             print(['interaction-' contrast '_measure-' measurelabel '_fwhm-' num2str(sizeoffwhm) '_FDR'] ,'-djpeg','-r0'); close
-            disp('Interaction: FDR'); 
+            disp('Interaction: FDR');
             qvaluesstruct = qval;
             save(['interaction-' contrast '_measure-' measurelabel '_fwhm-' num2str(sizeoffwhm) '_FDR.mat'],'qvaluesstruct');
         end
@@ -437,4 +438,3 @@ switch glmtype
     otherwise
         error('Check out if you define the glmtype flag correctly, or define your own general linear model, e,g MGLM')
 end
-
diff --git a/clinica/pipelines/statistics_surface/info.json b/clinica/pipelines/statistics_surface/info.json
index 590a99c4..4ce35066 100644
--- a/clinica/pipelines/statistics_surface/info.json
+++ b/clinica/pipelines/statistics_surface/info.json
@@ -5,5 +5,11 @@
     "space_wd": "20M",
     "credits": ["Arnaud Marcoux","Alexandre Routier"],
     "version": "0.1",
-    "dependencies": []
+    "dependencies": [
+        {
+            "type": "software",
+            "name": "freesurfer",
+            "version": ">=6.0.0"
+        }
+    ]
 }
diff --git a/clinica/pipelines/statistics_surface/statistics_surface_cli.py b/clinica/pipelines/statistics_surface/statistics_surface_cli.py
index 746628b1..dff6a80f 100644
--- a/clinica/pipelines/statistics_surface/statistics_surface_cli.py
+++ b/clinica/pipelines/statistics_surface/statistics_surface_cli.py
@@ -18,62 +18,63 @@ class StatisticsSurfaceCLI(ce.CmdParser):
         """Define the sub-command arguments."""
         from clinica.engine.cmdparser import PIPELINE_CATEGORIES
         from colorama import Fore
-        # Clinica compulsory arguments (e.g. BIDS, CAPS, group_id)
-        clinica_comp = self._args.add_argument_group('%sMandatory arguments%s' % (Fore.BLUE, Fore.RESET))
+
+        # Clinica compulsory arguments
+        clinica_comp = self._args.add_argument_group(PIPELINE_CATEGORIES['CLINICA_COMPULSORY'])
         clinica_comp.add_argument("caps_directory",
                                   help='Path to the CAPS directory.')
+        clinica_comp.add_argument("group_label",
+                                  help='User-defined identifier for the provided group of subjects.')
+        clinica_comp.add_argument("orig_input_data",
+                                  help='Type of surface-based feature: type '
+                                  '\'t1-freesurfer\' to use cortical thickness, '
+                                  '\'pet-surface\' to use projected FDG-PET data or '
+                                  '\'custom-pipeline\' to use you own data in CAPS directory '
+                                  '(see Wiki for details).',
+                                  choices=['t1-freesurfer', 'pet-surface', 'custom-pipeline'])
+        clinica_comp.add_argument("glm_type",
+                                  help='String based on the GLM type for the hypothesis. '
+                                  'You can choose between \'group_comparison\' and \'correlation\'.',
+                                  choices=['group_comparison', 'correlation'],
+                                  )
         clinica_comp.add_argument("subject_visits_with_covariates_tsv",
                                   help='TSV file containing a list of subjects with their sessions and all '
                                        'the covariates and factors needed for the GLM.')
-        clinica_comp.add_argument("design_matrix",
-                                  help='String to define the design matrix that fits into the GLM, '
-                                       'e.g. 1 + group + sex + age.')
         clinica_comp.add_argument("contrast",
                                   help='String to define the contrast matrix for the GLM, e.g. group. Please note '
                                        'that, when you want to perform negative correlation, the sign is ignored '
                                        'by the command line.')
-        clinica_comp.add_argument("string_format",
-                                  help='String to define the format of the columns in the TSV file, e.g., '
-                                       '%%s %%s %%s %%f if the columns contain a string, a string, a string '
-                                       'and a number, respectively.')
-        clinica_comp.add_argument("group_id",
-                                  help='User-defined identifier for the provided group of subjects.')
-        clinica_comp.add_argument("glm_type",
-                                  help='String based on the GLM type for the hypothesis. You can choose '
-                                       'between group_comparison and correlation.')
         # Optional arguments (e.g. FWHM)
         optional = self._args.add_argument_group(PIPELINE_CATEGORIES['OPTIONAL'])
+        optional.add_argument("-c", "--covariates",
+                              type=str, default=None, metavar='\'covariate_1 ... covariate_n\'',
+                              help='Covariates of the form --covariates \'covariate_1 ... covariate_n\'. '
+                                   'Each covariate must match the column name of the TSV file. '
+                                   'On default, no covariates are taken.')
         optional.add_argument("-fwhm", "--full_width_at_half_maximum",
                               type=int, default=20,
                               help='FWHM for the surface smoothing '
                                    '(default: --full_width_at_half_maximum %(default)s).')
-        optional.add_argument("-ft", "--feature_type",
-                              type=str, default='cortical_thickness',
-                              help='Type of surface-based feature: cortical_thickness (from t1-freesurfer pipeline)'
-                                   'or pet_fdg_projection (from pet-surface pipeline)'
-                                   '(default: --feature_type %(default)s).')
+        # Optional arguments for custom pipeline
+        opt_custom_input = self._args.add_argument_group(
+            '%sPipeline options if you selected custom-pipeline%s' % (Fore.BLUE, Fore.RESET)
+        )
+        opt_custom_input.add_argument("-cf", "--custom_file",
+                                      type=str, default=None,
+                                      help='Pattern of file inside CAPS directory using @subject, @session, '
+                                           '@fwhm, @hemi. '
+                                           'This flag must be specified with the --measure_label flag). '
+                                           'See Wiki for an example.')
+        opt_custom_input.add_argument("-ml", "--measure_label",
+                                      type=str, default=None,
+                                      help='Name of the feature type, it will be saved on the CAPS '
+                                           '_measure-FEATURE_LABEL key-value association. '
+                                           'This flag must be specified with the --custom_file flag). '
+                                           'See Wiki for an example.')
         # Clinica standard arguments (e.g. --n_procs)
         self.add_clinica_standard_arguments(add_tsv_flag=False)
         # Advanced arguments (i.e. tricky parameters)
         advanced = self._args.add_argument_group(PIPELINE_CATEGORIES['ADVANCED'])
-        advanced.add_argument("-cf", "--custom_file",
-                              type=str, default=None,
-                              help='Pattern of file inside CAPS directory using @subject, @session, @fwhm, @hemi. '
-                                   'No --feature_type must be specified in order to use this flag. '
-                                   'If you use this flag, you must specify a label with the --feature_label flag). '
-                                   'See Wiki for an example.')
-        advanced.add_argument("-fl", "--feature_label",
-                              type=str, default=None,
-                              help='Name of the feature type, it will be saved on the CAPS _measure-FEATURE_LABEL '
-                                   'key-value association.')
-        advanced.add_argument("-tup", "--threshold_uncorrected_pvalue",
-                              type=float, default=0.001,
-                              help='Threshold to display the uncorrected p-value '
-                                   '(--threshold_uncorrected_pvalue %(default)s).')
-        advanced.add_argument("-tcp", "--threshold_corrected_pvalue",
-                              type=float, default=0.05,
-                              help='Threshold to display the corrected p-value '
-                                   '(default: --threshold_corrected_pvalue %(default)s)')
         advanced.add_argument("-ct", "--cluster_threshold",
                               type=float, default=0.001,
                               help='Threshold to define a cluster in the process of cluster-wise correction '
@@ -81,62 +82,34 @@ class StatisticsSurfaceCLI(ce.CmdParser):
 
     def run_command(self, args):
         """Run the pipeline with defined args."""
-        import os
         from networkx import Graph
         from .statistics_surface_pipeline import StatisticsSurface
-        from .statistics_surface_utils import (check_inputs,
-                                               get_t1_freesurfer_custom_file,
+        from .statistics_surface_utils import (get_t1_freesurfer_custom_file,
                                                get_fdg_pet_surface_custom_file)
-        from clinica.utils.stream import cprint
         from clinica.utils.ux import print_end_pipeline, print_crash_files_and_exit
         from clinica.utils.exceptions import ClinicaException
 
-        if args.feature_type is not None:
-            if args.custom_file is not None:
-                raise ClinicaException('--feature_type and --custom_file are mutually exclusive: you must choose '
-                                       'between one or the other. See documentation for more information.')
-            if args.feature_label is not None:
-                raise ClinicaException('--feature_label should not be used with --feature_type.')
-            # FreeSurfer cortical thickness
-            if args.feature_type == 'cortical_thickness':
-                args.custom_file = get_t1_freesurfer_custom_file()
-                args.feature_label = 'ct'
-            # PET cortical projection
-            elif args.feature_type == 'pet_fdg_projection':
-                args.custom_file = get_fdg_pet_surface_custom_file()
-                args.feature_label = 'fdg'
-            else:
-                raise ClinicaException('Feature type ' + args.feature_type + ' not recognized. Use --custom_file '
-                                       'to specify your own files (without --feature_type).')
-        elif args.feature_type is None:
-            if args.custom_file is None:
-                cprint('No feature type selected: using cortical thickness as default value')
-                args.custom_file = get_t1_freesurfer_custom_file()
-                args.feature_label = 'ct'
-            else:
-                cprint('Using custom features.')
-                if args.feature_label is None:
-                    raise ClinicaException('You must specify a --feature_label when using the --custom_files flag.')
+        # FreeSurfer cortical thickness
+        if args.orig_input_data == 't1-freesurfer':
+            args.custom_file = get_t1_freesurfer_custom_file()
+            args.measure_label = 'ct'
+        # PET cortical projection
+        elif args.orig_input_data == 'pet-surface':
+            args.custom_file = get_fdg_pet_surface_custom_file()
+            args.measure_label = 'fdg'
+        else:
+            if (args.custom_file is None) or (args.measure_label is None):
+                raise ClinicaException('You must set --measure_label and --custom_file flags.')
 
-        # Check if the group label has been existed, if yes, give an error to the users
-        # Note(AR): if the user wants to compare Cortical Thickness measure with PET measure
-        # using the group_id, Clinica won't allow it.
-        # TODO: Modify this behaviour
-        if os.path.exists(os.path.join(os.path.abspath(self.absolute_path(args.caps_directory)), 'groups', 'group-' + args.group_id)):
-            error_message = 'group_id: ' + args.group_id + ' already exists, please choose another one or delete ' \
-                            'the existing folder and also the working directory and rerun the pipeline'
-            raise ClinicaException(error_message)
         parameters = {
-            'group_label': args.group_id,
-            'design_matrix': args.design_matrix,
+            'orig_input_data': args.orig_input_data,
+            'group_label': args.group_label,
+            'covariates': args.covariates,
             'contrast': args.contrast,
-            'str_format': args.string_format,
             'glm_type': args.glm_type,
             'custom_file': args.custom_file,
-            'feature_label': args.feature_label,
+            'measure_label': args.measure_label,
             'full_width_at_half_maximum': args.full_width_at_half_maximum,
-            'threshold_uncorrected_pvalue': args.threshold_uncorrected_pvalue,
-            'threshold_corrected_pvalue': args.threshold_corrected_pvalue,
             'cluster_threshold': args.cluster_threshold,
         }
         pipeline = StatisticsSurface(
@@ -147,14 +120,6 @@ class StatisticsSurfaceCLI(ce.CmdParser):
             name=self.name
         )
 
-        check_inputs(pipeline.caps_directory,
-                     pipeline.parameters['custom_file'],
-                     pipeline.parameters['full_width_at_half_maximum'],
-                     pipeline.tsv_file)
-
-        cprint("Parameters used for this pipeline:")
-        cprint(pipeline.parameters)
-
         if args.n_procs:
             exec_pipeline = pipeline.run(plugin='MultiProc',
                                          plugin_args={'n_procs': args.n_procs})
diff --git a/clinica/pipelines/statistics_surface/statistics_surface_pipeline.py b/clinica/pipelines/statistics_surface/statistics_surface_pipeline.py
index 23f89282..7e0d6484 100644
--- a/clinica/pipelines/statistics_surface/statistics_surface_pipeline.py
+++ b/clinica/pipelines/statistics_surface/statistics_surface_pipeline.py
@@ -5,49 +5,15 @@ import clinica.pipelines.engine as cpe
 
 class StatisticsSurface(cpe.Pipeline):
     """
-    Based on the Matlab toolbox [SurfStat](http://www.math.mcgill.ca/keith/surfstat/), which performs statistical
-    analyses of univariate and multivariate surface and volumetric data using the generalized linear model (GLM),
-    this pipeline performs analyses including group comparison and correlation with surface-based features e.g.
-    cortical thickness from t1-freesurfer or map of activity from PET data from pet-surface pipeline.
-
-    TODO: Refactor StatisticsSurface
-        [ ] build_input_node
-            [ ] Remove current read_parameters_node
-            [ ] With the help of statistics_surface_utils.py::check_inputs, use new clinica_file_reader function
-                to check and extract surface-based features
-            [ ] Delete statistics_surface_utils.py::check_inputs function
-            [ ] Move statistics_surface_cli.py checks of input data in this method
-            [ ] Handle overwrite case
-            [ ] Display participants, covariates and info regarding the GLM.
-        [ ] build_core_nodes
-            [ ] Use surfstat.inputs.full_width_at_half_maximum = self.parameter['full_width_at_half_maximum']
-                instead of connecting read_parameters_node.inputs.full_width_at_half_maximum
-                to surfstat.inputs.full_width_at_half_maximum
-            [ ] Repeat for other keys
-            [ ] Use working directory
-        [ ] build_output_node
-            [ ] Remove path_to_matscript and freesurfer_home: it should be set in runmatlab function
-            [ ] Copy results from <WD> to <CAPS>
-        [ ] Clean/adapt statistics_surface_utils.py
-
-    Args:
-        caps_directory: str, the output folder of recon-all which will contain the result files: ?h.thickness.fwhm**.mgh.
-        tsv_file: str, Path to the tsv containing the information for GLM.
-        design_matrix: str, the linear model that fits into the GLM, for example '1+group'.
-        contrast: string, the contrast matrix for GLM, if the factor you choose is categorized variable, clinica_surfstat will create two contrasts,
-                  for example, contrast = 'Label', this will create contrastpos = Label.AD - Label.CN, contrastneg = Label.CN - Label.AD; if the fac-
-                  tory that you choose is a continuous factor, clinica_surfstat will just create one contrast, for example, contrast = 'Age', but note,
-                  the string name that you choose should be exactly the same with the columns names in your subjects_visits_tsv.
-        str_format: string, the str_format which uses to read your tsv file, the typy of the string should corresponds exactly with the columns in the tsv file.
-            Defaut parameters, we set these parameters to be some default values, but you can also set it by yourself:
-        group_label: current group name for this analysis
-        glm_type: based on the hypothesis, you should define one of the glm types, "group_comparison", "correlation"
-        full_width_at_half_maximum: fwhm for the surface smoothing, default is 20, integer.
-        threshold_uncorrected_pvalue: threshold to display the uncorrected Pvalue, float, default is 0.001.
-        threshold_corrected_pvalue: the threshold to display the corrected cluster, default is 0.05, float.
-        cluster_threshold: threshold to define a cluster in the process of cluster-wise correction, default is 0.001, float.
-        working_directory: define where to put the infomartion of the nipype workflow.
-        n_procs: define how many cores to run this workflow.
+    StatisticsSurface - Surface-based mass-univariate analysis with SurfStat
+
+    See documentation at http://www.clinica.run/doc/Pipelines/Stats_Surface/
+
+    Note:
+        The `tsv_file` attribute is overloaded for this pipeline. It must contain a list of subjects
+        with their sessions and all the covariates and factors needed for the GLM.
+
+        Pipeline parameters are explained in StatisticsSurfaceCLI.define_options()
 
     Returns:
         A clinica pipeline object containing the StatisticsSurface pipeline.
@@ -58,16 +24,16 @@ class StatisticsSurface(cpe.Pipeline):
         from clinica.utils.exceptions import ClinicaException
         from clinica.utils.group import check_group_label
 
+        if 'orig_input_data' not in self.parameters.keys():
+            raise KeyError('Missing compulsory orig_input_data key in pipeline parameter.')
+        if 'covariates' not in self.parameters.keys():
+            self.parameters['covariates'] = None
         if 'custom_file' not in self.parameters.keys():
             self.parameters['custom_file'] = get_t1_freesurfer_custom_file()
-        if 'feature_label' not in self.parameters.keys():
-            self.parameters['feature_label'] = 'ct',
+        if 'measure_label' not in self.parameters.keys():
+            self.parameters['measure_label'] = 'ct',
         if 'full_width_at_half_maximum' not in self.parameters.keys():
             self.parameters['full_width_at_half_maximum'] = 20
-        if 'threshold_uncorrected_pvalue' not in self.parameters.keys():
-            self.parameters['threshold_uncorrected_pvalue'] = 0.001
-        if 'threshold_corrected_pvalue' not in self.parameters.keys():
-            self.parameters['threshold_corrected_pvalue'] = 0.05,
         if 'cluster_threshold' not in self.parameters.keys():
             self.parameters['cluster_threshold'] = 0.001,
 
@@ -79,12 +45,6 @@ class StatisticsSurface(cpe.Pipeline):
             raise ClinicaException(
                 "FWHM for the surface smoothing you specified is wrong: it should be 0, 5, 10, 15 or 20 "
                 "(given value: %s)." % self.parameters['full_width_at_half_maximum'])
-        if self.parameters['threshold_uncorrected_pvalue'] < 0 or self.parameters['threshold_uncorrected_pvalue'] > 1:
-            raise ClinicaException("Uncorrected p-value threshold should be a lower than 1 "
-                                   "(given value: %s)." % self.parameters['threshold_uncorrected_pvalue'])
-        if self.parameters['threshold_corrected_pvalue'] < 0 or self.parameters['threshold_corrected_pvalue'] > 1:
-            raise ClinicaException("Corrected p-value threshold should be between 0 and 1 "
-                                   "(given value: %s)." % self.parameters['threshold_corrected_pvalue'])
         if self.parameters['cluster_threshold'] < 0 or self.parameters['cluster_threshold'] > 1:
             raise ClinicaException("Cluster threshold should be between 0 and 1 "
                                    "(given value: %s)." % self.parameters['cluster_threshold'])
@@ -100,18 +60,7 @@ class StatisticsSurface(cpe.Pipeline):
         Returns:
             A list of (string) input fields name.
         """
-
-        return ['design_matrix',
-                'contrast',
-                'str_format',
-                'group_label',
-                'glm_type',
-                'surface_file',
-                'full_width_at_half_maximum',
-                'threshold_uncorrected_pvalue',
-                'threshold_corrected_pvalue',
-                'cluster_threshold',
-                'feature_label']
+        return []
 
     def get_output_fields(self):
         """Specify the list of possible outputs of this pipelines.
@@ -119,134 +68,134 @@ class StatisticsSurface(cpe.Pipeline):
         Returns:
             A list of (string) output fields name.
         """
-
-        return []
+        return ['output_dir']
 
     def build_input_node(self):
         """Build and connect an input node to the pipelines.
         """
+        import os
+        from clinica.utils.inputs import clinica_file_reader
+        from clinica.utils.exceptions import ClinicaException
+        from clinica.utils.stream import cprint
+
+        # Check if already present in CAPS
+        # ================================
+        # Check if the group label has been existed, if yes, give an error to the users
+        # Note(AR): if the user wants to compare Cortical Thickness measure with PET measure
+        # using the group_id, Clinica won't allow it.
+        # TODO: Modify this behaviour
+        if os.path.exists(os.path.join(self.caps_directory, 'groups', 'group-' + self.parameters['group_label'])):
+            error_message = ('Group ID %s already exists, please choose another one or delete the existing folder and '
+                             'also the working directory and rerun the pipeline') % self.parameters['group_label']
+            raise ClinicaException(error_message)
+            # statistics_dir_tsv = os.path.join(input_directory, 'groups', group_id, 'statistics', 'participant.tsv')
+            # # Copy the subjects_visits_tsv to the result folder
+            # # First, check if the subjects_visits_tsv has the same info with the participant.tsv in the folder of statistics.
+            # # If the participant TSV does not exit, copy subjects_visits_tsv in the folder of statistics too,
+            # # if it is here, compare them.
+            # if not os.path.isfile(statistics_dir_tsv):
+            #     copy(subjects_visits_tsv, statistics_dir_tsv)
+            # else:
+            #     # Compare the two TSV files
+            #     if not have_same_subjects(statistics_dir_tsv, subjects_visits_tsv):
+            #         raise ValueError("It seems that this round of analysis does not contain the same subjects"
+            #                          "where you want to put the results, please check it!")
+
+        # Check input files before calling SurfStat with Matlab
+        # =====================================================
+        all_errors = []
+        # clinica_files_reader expects regexp to start at subjects/ so sub-*/ses-*/ is removed here
+        pattern_hemisphere = self.parameters['custom_file'].replace(
+            '@subject', 'sub-*').replace(
+            '@session', 'ses-*').replace(
+            '@fwhm', str(self.parameters['full_width_at_half_maximum'])).replace(
+            'sub-*/ses-*/', ''
+        )
+        # Files on left hemisphere
+        lh_surface_based_info = {
+            'pattern': pattern_hemisphere.replace('@hemi', 'lh'),
+            'description': 'surface-based features on left hemisphere at FWHM = %s' %
+                           self.parameters['full_width_at_half_maximum'],
+        }
+        try:
+            clinica_file_reader(self.subjects, self.sessions, self.caps_directory, lh_surface_based_info)
+        except ClinicaException as e:
+            all_errors.append(e)
+        rh_surface_based_info = {
+            'pattern': pattern_hemisphere.replace('@hemi', 'rh'),
+            'description': 'surface-based features on right hemisphere at FWHM = %s' %
+                           self.parameters['full_width_at_half_maximum'],
+        }
+        try:
+            clinica_file_reader(self.subjects, self.sessions, self.caps_directory, rh_surface_based_info)
+        except ClinicaException as e:
+            all_errors.append(e)
+        # Raise all errors if something happened
+        if len(all_errors) > 0:
+            error_message = 'Clinica faced errors while trying to read files in your CAPS directory.\n'
+            for msg in all_errors:
+                error_message += str(msg)
+            raise RuntimeError(error_message)
+
+        # Give pipeline info
+        # ==================
+        cprint('The pipeline will last a few minutes. Images generated by Matlab will popup during the pipeline.')
 
+    def build_output_node(self):
+        """Build and connect an output node to the pipeline."""
         import nipype.interfaces.utility as nutil
         import nipype.pipeline.engine as npe
-
-        read_parameters_node = npe.Node(name="LoadingCLIArguments",
-                                        interface=nutil.IdentityInterface(
-                                            fields=self.get_input_fields(),
-                                            mandatory_inputs=True))
-        read_parameters_node.inputs.design_matrix = self.parameters['design_matrix']
-        read_parameters_node.inputs.contrast = self.parameters['contrast']
-        read_parameters_node.inputs.str_format = self.parameters['str_format']
-        read_parameters_node.inputs.group_label = self.parameters['group_label']
-        read_parameters_node.inputs.glm_type = self.parameters['glm_type']
-        read_parameters_node.inputs.surface_file = self.parameters['custom_file']
-        read_parameters_node.inputs.full_width_at_half_maximum = self.parameters['full_width_at_half_maximum']
-        read_parameters_node.inputs.threshold_uncorrected_pvalue = self.parameters['threshold_uncorrected_pvalue']
-        read_parameters_node.inputs.threshold_corrected_pvalue = self.parameters['threshold_corrected_pvalue']
-        read_parameters_node.inputs.cluster_threshold = self.parameters['cluster_threshold']
-        read_parameters_node.inputs.feature_label = self.parameters['feature_label']
+        from .statistics_surface_utils import save_to_caps
+
+        # Writing results into CAPS
+        # =========================
+        save_to_caps = npe.Node(interface=nutil.Function(
+            input_names=['source_dir', 'caps_dir', 'overwrite_caps', 'pipeline_parameters'],
+            function=save_to_caps),
+            name='SaveToCaps')
+        save_to_caps.inputs.caps_dir = self.caps_directory
+        save_to_caps.inputs.overwrite_caps = self.overwrite_caps
+        save_to_caps.inputs.pipeline_parameters = self.parameters
 
         self.connect([
-            (read_parameters_node, self.input_node, [('design_matrix',                'design_matrix')]),  # noqa
-            (read_parameters_node, self.input_node, [('surface_file',                 'surface_file')]),  # noqa
-            (read_parameters_node, self.input_node, [('contrast',                     'contrast')]),  # noqa
-            (read_parameters_node, self.input_node, [('str_format',                   'str_format')]),  # noqa
-            (read_parameters_node, self.input_node, [('group_label',                  'group_label')]),  # noqa
-            (read_parameters_node, self.input_node, [('glm_type',                     'glm_type')]),  # noqa
-            (read_parameters_node, self.input_node, [('full_width_at_half_maximum',   'full_width_at_half_maximum')]),  # noqa
-            (read_parameters_node, self.input_node, [('threshold_uncorrected_pvalue', 'threshold_uncorrected_pvalue')]),  # noqa
-            (read_parameters_node, self.input_node, [('threshold_corrected_pvalue',   'threshold_corrected_pvalue')]),  # noqa
-            (read_parameters_node, self.input_node, [('cluster_threshold',            'cluster_threshold')]),  # noqa
-            (read_parameters_node, self.input_node, [('feature_label',                'feature_label')])
+            (self.output_node, save_to_caps, [('output_dir', 'source_dir')]),
         ])
 
-    def build_output_node(self):
-        """Build and connect an output node to the pipelines.
-        """
-
-        pass
-
     def build_core_nodes(self):
-        """Build and connect the core nodes of the pipelines.
-        """
-
-        import clinica.pipelines.statistics_surface.statistics_surface_utils as utils
+        """Build and connect the core nodes of the pipeline."""
+        import os
         import nipype.interfaces.utility as nutil
         import nipype.pipeline.engine as npe
-        from nipype.interfaces.io import JSONFileSink
-
-        # Node to fetch the input variables
-        data_prep = npe.Node(name='inputnode',
-                             interface=nutil.Function(
-                                 input_names=['input_directory', 'subjects_visits_tsv', 'group_label', 'glm_type'],
-                                 output_names=['path_to_matscript', 'surfstat_input_dir', 'output_directory',
-                                               'freesurfer_home', 'out_json'],
-                                 function=utils.prepare_data))
-        data_prep.inputs.input_directory = self.caps_directory
-        data_prep.inputs.subjects_visits_tsv = self.tsv_file
+        from .statistics_surface_utils import init_input_node
+        import clinica.pipelines.statistics_surface.statistics_surface_utils as utils
+
+        init_input = npe.Node(
+            interface=nutil.Function(
+                input_names=['parameters', 'base_dir', 'subjects_visits_tsv'],
+                output_names=['group_label', 'surfstat_results_dir'],
+                function=init_input_node),
+            name='0-InitPipeline')
+        init_input.inputs.parameters = self.parameters
+        init_input.inputs.base_dir = os.path.join(self.base_dir, self.name)
+        init_input.inputs.subjects_visits_tsv = self.tsv_file
 
         # Node to wrap the SurfStat matlab script
-        surfstat = npe.Node(name='surfstat',
+        surfstat = npe.Node(name='1-RunSurfStat',
                             interface=nutil.Function(
-                                input_names=['input_directory',
-                                             'output_directory',
+                                input_names=['caps_dir',
+                                             'output_dir',
                                              'subjects_visits_tsv',
-                                             'design_matrix',
-                                             'contrast',
-                                             'str_format',
-                                             'glm_type',
-                                             'group_label',
-                                             'freesurfer_home',
-                                             'surface_file',
-                                             'path_to_matscript',
-                                             'full_width_at_half_maximum',
-                                             'threshold_uncorrected_pvalue',
-                                             'threshold_corrected_pvalue',
-                                             'cluster_threshold',
-                                             'feature_label'],
-                                output_names=['out_images'],
-                                function=utils.runmatlab))
+                                             'pipeline_parameters',
+                                             ],
+                                output_names=['output_dir'],
+                                function=utils.run_matlab))
+        surfstat.inputs.caps_dir = self.caps_directory
         surfstat.inputs.subjects_visits_tsv = self.tsv_file
-
-        # Node to create the dictionary for JSONFileSink
-        jsondict = npe.Node(name='Jsondict',
-                            interface=nutil.Function(
-                                input_names=['glm_type', 'design_matrix', 'str_format', 'contrast', 'group_label', 'full_width_at_half_maximum',
-                                             'threshold_uncorrected_pvalue', 'threshold_corrected_pvalue', 'cluster_threshold'],
-                                output_names=['json_dict'],
-                                function=utils.json_dict_create))
-
-        # Node to write the GLM information into a JSON file
-        jsonsink = npe.Node(JSONFileSink(input_names=['out_file']), name='jsonsinker')
+        surfstat.inputs.pipeline_parameters = self.parameters
 
         # Connection
         # ==========
         self.connect([
-            (self.input_node, data_prep, [('group_label', 'group_label')]),  # noqa
-            (self.input_node, data_prep, [('glm_type', 'glm_type')]),  # noqa
-            (self.input_node, surfstat, [('design_matrix', 'design_matrix')]),  # noqa
-            (self.input_node, surfstat, [('contrast', 'contrast')]),  # noqa
-            (self.input_node, surfstat, [('str_format', 'str_format')]),  # noqa
-            (self.input_node, surfstat, [('glm_type', 'glm_type')]),  # noqa
-            (self.input_node, surfstat, [('group_label', 'group_label')]),  # noqa
-            (self.input_node, surfstat, [('full_width_at_half_maximum', 'full_width_at_half_maximum')]),  # noqa
-            (self.input_node, surfstat, [('threshold_uncorrected_pvalue', 'threshold_uncorrected_pvalue')]),  # noqa
-            (self.input_node, surfstat, [('threshold_corrected_pvalue', 'threshold_corrected_pvalue')]),  # noqa
-            (self.input_node, surfstat, [('cluster_threshold', 'cluster_threshold')]),  # noqa
-            (self.input_node, surfstat, [('surface_file', 'surface_file')]),  # noqa
-            (self.input_node, surfstat, [('feature_label', 'feature_label')]),
-            (data_prep, surfstat, [('surfstat_input_dir', 'input_directory')]),  # noqa
-            (data_prep, surfstat, [('path_to_matscript', 'path_to_matscript')]),  # noqa
-            (data_prep, surfstat, [('output_directory', 'output_directory')]),  # noqa
-            (data_prep, surfstat, [('freesurfer_home', 'freesurfer_home')]),  # noqa
-            (self.input_node, jsondict, [('glm_type', 'glm_type')]),  # noqa
-            (self.input_node, jsondict, [('design_matrix', 'design_matrix')]),  # noqa
-            (self.input_node, jsondict, [('str_format', 'str_format')]),  # noqa
-            (self.input_node, jsondict, [('contrast', 'contrast')]),  # noqa
-            (self.input_node, jsondict, [('group_label', 'group_label')]),  # noqa
-            (self.input_node, jsondict, [('full_width_at_half_maximum', 'full_width_at_half_maximum')]),  # noqa
-            (self.input_node, jsondict, [('threshold_uncorrected_pvalue', 'threshold_uncorrected_pvalue')]),  # noqa
-            (self.input_node, jsondict, [('threshold_corrected_pvalue', 'threshold_corrected_pvalue')]),  # noqa
-            (self.input_node, jsondict, [('cluster_threshold', 'cluster_threshold')]),  # noqa
-            (data_prep, jsonsink, [('out_json', 'out_file')]),  # noqa
-            (jsondict, jsonsink, [('json_dict', 'in_dict')]),  # noqa
+            (init_input, surfstat, [('surfstat_results_dir', 'output_dir')]),
+            (surfstat, self.output_node, [('output_dir', 'output_dir')]),
         ])
diff --git a/clinica/pipelines/statistics_surface/statistics_surface_utils.py b/clinica/pipelines/statistics_surface/statistics_surface_utils.py
index 85166da6..8fc396d0 100644
--- a/clinica/pipelines/statistics_surface/statistics_surface_utils.py
+++ b/clinica/pipelines/statistics_surface/statistics_surface_utils.py
@@ -27,246 +27,258 @@ def get_fdg_pet_surface_custom_file():
     return custom_file
 
 
-def prepare_data(input_directory, subjects_visits_tsv, group_label, glm_type):
-    """Fetch all the intermediate variables for this workflow.
+def init_input_node(parameters, base_dir, subjects_visits_tsv):
+    """Initialize the pipeline.
+
+    This function will:
+        - Create `surfstat_results_dir` in `base_dir`/<group_id> for SurfStat;
+        - Save pipeline parameters in JSON file;
+        - Copy TSV file with covariates;
+        - Print begin execution message.
+    """
+    import os
+    import errno
+    import json
+    import shutil
+    from clinica.utils.ux import print_begin_image
+    from clinica.pipelines.statistics_surface.statistics_surface_utils import create_glm_info_dictionary
+
+    group_id = 'group-' + parameters['group_label']
+
+    # Create surfstat_results_dir for SurfStat
+    surfstat_results_dir = os.path.join(base_dir, group_id)
+    try:
+        os.makedirs(surfstat_results_dir)
+    except OSError as e:
+        if e.errno != errno.EEXIST:  # EEXIST: folder already exists
+            raise e
+
+    # Save pipeline parameters in JSON file
+    glm_dict = create_glm_info_dictionary(subjects_visits_tsv, parameters)
+    json_filename = os.path.join(surfstat_results_dir, group_id + '_glm.json')
+    with open(json_filename, 'w') as json_file:
+        json.dump(glm_dict, json_file, indent=4)
+
+    # Copy TSV file with covariates
+    tsv_filename = os.path.join(surfstat_results_dir, group_id + '_covariates.tsv')
+    shutil.copyfile(subjects_visits_tsv, tsv_filename)
+
+    # Print begin message
+    list_keys = ['AnalysisType', 'Covariates', 'Contrast', 'FWHM', 'ClusterThreshold']
+    list_values = [
+        parameters['glm_type'],
+        parameters['covariates'],
+        parameters['contrast'],
+        str(parameters['full_width_at_half_maximum']),
+        str(parameters['cluster_threshold'])
+    ]
+    group_id = 'group-' + parameters['group_label']
+    print_begin_image(group_id, list_keys, list_values)
+
+    return parameters['group_label'], surfstat_results_dir
+
+
+def get_string_format_from_tsv(tsv_file):
+    """
+    Determine string format from TSV file.
+
+    If the TSV file is like:
+
+    participant_id  session_id  sex     group   age
+    sub-CLNC0001    ses-M00     Female  CN      71.1
+    sub-CLNC0002    ses-M00     Male    CN      81.3
+    sub-CLNC0003    ses-M00     Male    CN      75.4
+
+    The columns of the TSV file contains consecutively strings, strings,
+    strings, strings and float. The string_format is therefore "%s %s %s %s %f".
 
     Args:
-        input_directory: CAPS directory
-        subjects_visits_tsv: TSV file defining the GLM model
-        group_label: Current group name for this analysis
-        glm_type: Based on the hypothesis, you should define one of the glm types, "group_comparison", "correlation"
+        tsv_file: TSV file.
 
     Returns:
-
+        String formatting of the TSV file (e.g. "%s %s %s %s %f")
     """
-    import os
-    from shutil import copy
-    import clinica.pipelines as clp
-    import sys
     import pandas as pd
-    from clinica.utils.stream import cprint
-
-    path_to_matscript = os.path.join(os.path.dirname(clp.__path__[0]), 'lib', 'clinicasurfstat')
-
-    # CAPS input and output vars
-    input_directory = os.path.expanduser(input_directory)
-    surfstat_input_dir = os.path.join(input_directory, 'subjects')
-
-    group_id = 'group-' + group_label
-    statistics_dir_tsv = os.path.join(input_directory, 'groups', group_id, 'statistics', 'participant.tsv')
-
-    if glm_type == "group_comparison":
-        output_directory = os.path.join(input_directory, 'groups', group_id, 'statistics', 'surfstat_group_comparison')
-        if not os.path.exists(output_directory):
-            try:
-                os.makedirs(output_directory)
-            except BaseException:
-                raise OSError("SurfStat: can't create destination directory (%s)!" % (output_directory))
-    elif glm_type == "correlation":
-        output_directory = os.path.join(input_directory, 'groups', group_id, 'statistics', 'surfstat_correlation_analysis')
-        if not os.path.exists(output_directory):
-            try:
-                os.makedirs(output_directory)
-            except BaseException:
-                raise OSError("SurfStat: can't create destination directory (%s)!" % (output_directory))
-    else:
-        cprint("The other GLM situations have not been implemented in this pipeline.")
-        sys.exit()
-
-    # Copy the subjects_visits_tsv to the result folder
-    # First, check if the subjects_visits_tsv has the same info with the participant.tsv in the folder of statistics.
-    # If the participant tsv does not exit, cp subjects_visits_tsv in the folder of statistics too, if it is here, compare them.
-    if not os.path.isfile(statistics_dir_tsv):
-        copy(subjects_visits_tsv, statistics_dir_tsv)
-    else:
-        # Compare the two TSV files
-        participant_df = pd.io.parsers.read_csv(statistics_dir_tsv, sep='\t')
-        subjects_visits_tsv_df = pd.io.parsers.read_csv(subjects_visits_tsv, sep='\t')
-        participant_list = list(participant_df.participant_id)
-        subjects_visits_tsv_list = list(subjects_visits_tsv_df.participant_id)
-        dif_list = list(set(participant_list) - set(subjects_visits_tsv_list))
-        try:
-            len(dif_list) == 0
-        except BaseException:
-            raise ValueError("It seems that this round of analysis does not contain the same subjects where you want to put the results, please check it!")
-
-        group_tsv = 'group-' + group_label + '_participants.tsv'
-        copied_tsv = os.path.join(output_directory, group_tsv)
-        copy(subjects_visits_tsv, copied_tsv)
-    # Point to the path to the json file
-    out_json = os.path.join(output_directory, 'group-' + group_label + '_glm.json')
-
-    freesurfer_home = os.environ["FREESURFER_HOME"]
-
-    return path_to_matscript, surfstat_input_dir, output_directory, freesurfer_home, out_json
-
-
-def runmatlab(input_directory,
-              output_directory,
-              subjects_visits_tsv,
-              design_matrix, contrast,
-              str_format,
-              glm_type,
-              group_label,
-              freesurfer_home,
-              surface_file,
-              path_to_matscript,
-              full_width_at_half_maximum,
-              threshold_uncorrected_pvalue,
-              threshold_corrected_pvalue,
-              cluster_threshold,
-              feature_label):
-    """
-        a wrapper the matlab script of surfstat with nipype.
 
-    Args:
-        input_directory: surfstat_input_dir where containing all the subjects' output in CAPS directory
-        output_directory: output folder to contain the result in CAPS folder
-        subjects_visits_tsv: tsv file containing the glm information
-        design_matrix: str, the linear model that fits into the GLM, for example '1+group'.
-        contrast: string, the contrast matrix for GLM, if the factor you choose is categorized variable, clinica_surfstat will create two contrasts,
-                  for example, contrast = 'Label', this will create contrastpos = Label.AD - Label.CN, contrastneg = Label.CN - Label.AD; if the fac-
-                  tory that you choose is a continuous factor, clinica_surfstat will just create one contrast, for example, contrast = 'Age', but note,
-                  the string name that you choose should be exactly the same with the columns names in your subjects_visits_tsv.
-        str_format:string, the str_format which uses to read your tsv file, the type of the string should corresponds exactly with the columns in the tsv file.
-                  Defaut parameters, we set these parameters to be some default values, but you can also set it by yourself:
-        glm_type: based on the hypothesis, you should define one of the glm types, "group_comparison", "correlation"
-        group_label: current group name for this analysis
-        freesurfer_home: the environmental variable $FREESURFER_HOME
-        surface_file: Specify where to find the data surfaces file in the "CAPS/subject" directory, using specific keywords.
-                     For instance, to catch for each subject the cortical thickness, the string used will be:
-                     '@subject/@session/t1/freesurfer_cross_sectional/@subject_@session/surf/@hemi.thickness.fwhm@fwhm.fsaverage.mgh'
-                     More information is available on the documentation page of the surfstat pipelines. The keywords @subject @ session @hemi @fwhm
-                     represents the variable parts.
-        path_to_matscript: path to find the matlab script
-        full_width_at_half_maximum: fwhm for the surface smoothing, default is 20, integer.
-        threshold_uncorrected_pvalue: threshold to display the uncorrected Pvalue, float, default is 0.001.
-        threshold_corrected_pvalue: the threshold to display the corrected cluster, default is 0.05, float.
-        cluster_threshold: threshold to define a cluster in the process of cluster-wise correction, default is 0.001, float.
+    demographics_df = pd.read_csv(tsv_file, sep='\t')
+
+    def dtype_to_str_format(dtype):
+        """Convert pandas dtypes (e.g. int64) to string format (e.g. %d)"""
+        import numpy as np
+        if dtype == np.int64:
+            str_format = "%d"
+        elif dtype == np.float64:
+            str_format = "%f"
+        elif dtype == np.object:
+            str_format = "%s"
+        else:
+            raise ValueError("Unknown dtype (given: %s)" % dtype)
+        return str_format
+
+    list_str_format = []
+    for column in demographics_df.columns:
+        list_str_format.append(dtype_to_str_format(demographics_df[column].dtype))
+
+    return " ".join(list_str_format)
 
-    Returns:
 
+def covariates_to_design_matrix(contrast, covariates=None):
     """
-    from nipype.interfaces.matlab import MatlabCommand, get_matlab_command
-    from os.path import join
-    import sys
-    from clinica.utils.stream import cprint
+    Generate design matrix for SurfStat based on the contrast and the optional list of covariates.
+
+    Design matrix "1 + <contrast> + <covariate_1> + ... + <covariate_n>"
+
+    Example:
+        >>> from clinica.pipelines.statistics_surface.statistics_surface_utils import covariates_to_design_matrix
+        >>> covariates_to_design_matrix('group', 'age sex group')
+        1 + group + age + sex
+        >>> covariates_to_design_matrix('group', 'age')
+        1 + group + age
+        >>> covariates_to_design_matrix('group', None)
+        1 + group
+    """
+    if covariates:
+        # Convert string to list while handling case where several spaces are present
+        list_covariates = list(set(covariates.split(' ')))
+        try:
+            list_covariates.remove('')
+        except ValueError:
+            pass
+
+        if contrast in list_covariates:
+            design_matrix = "1 + " + " + ".join(covariate for covariate in list_covariates)
+        else:
+            design_matrix = "1 + " + contrast + " + " + " + ".join(covariate for covariate in list_covariates)
+    else:
+        design_matrix = "1 + " + contrast
 
-    MatlabCommand.set_default_matlab_cmd(
-        get_matlab_command())  # this is to set the matlab_path(os.environ) in your bashrc file, to choose which version of matlab do you wanna use
-    # here, set_default_matlab_cmd is a @classmethod
-    matlab = MatlabCommand()
+    return design_matrix
 
-    # add the dynamic traits
-    # openGL_trait = traits.Bool(True, argstr='-nosoftwareopengl', usedefault=True, desc='Switch on hardware openGL', nohash=True)
-    # matlab.input_spec.add_trait(matlab.input_spec(), 'nosoftwareopengl', openGL_trait() )
-    if sys.platform.startswith('linux'):
-        matlab.inputs.args = '-nosoftwareopengl'  # Bug, for my laptop, it does not work, but the command line does have the flag -nosoftwareopengl, we should try on other computer's matlab to check if this flag works!
-    matlab.inputs.paths = path_to_matscript  # CLINICA_HOME, this is the path to add into matlab, addpath
 
-    matlab.inputs.script = """
-    clinicasurfstat('%s', '%s', '%s', '%s', '%s', '%s', '%s', '%s', '%s', '%s', '%s', '%s', %d, '%s', %.3f, '%s', %.3f, '%s', %.3f);
-    """ % (input_directory, output_directory, subjects_visits_tsv, design_matrix, contrast, str_format, glm_type, group_label, freesurfer_home, surface_file, feature_label, 'sizeoffwhm',
-           full_width_at_half_maximum,
-           'thresholduncorrectedpvalue', threshold_uncorrected_pvalue, 'thresholdcorrectedpvalue',
-           threshold_corrected_pvalue, 'clusterthreshold',
-           cluster_threshold)  # here, we should define the inputs for the matlab function that you want to use
-    matlab.inputs.mfile = True  # this will create a file: pyscript.m , the pyscript.m is the default name
-    matlab.inputs.single_comp_thread = False  # this will stop runing with single thread
-    matlab.inputs.logfile = join(output_directory, "matlab_output.log")
-    cprint("Matlab logfile is located in the folder: %s" % matlab.inputs.logfile)
-    cprint("Matlab script command = %s" % matlab.inputs.script)
-    cprint("MatlabCommand inputs flag: single_comp_thread = %s" % matlab.inputs.single_comp_thread)
-    cprint("MatlabCommand choose which matlab to use(matlab_cmd): %s" % get_matlab_command())
-    if sys.platform.startswith('linux'):
-        cprint("MatlabCommand inputs flag: nosoftwareopengl = %s" % matlab.inputs.args)
-    out = matlab.run()
-    return out
-
-
-def json_dict_create(glm_type,
-                     design_matrix,
-                     str_format,
-                     contrast,
-                     group_label,
-                     full_width_at_half_maximum,
-                     threshold_uncorrected_pvalue,
-                     threshold_corrected_pvalue,
-                     cluster_threshold):
-    """Create a JSON file containing the GLM information.
+def run_matlab(caps_dir,
+               output_dir,
+               subjects_visits_tsv,
+               pipeline_parameters):
+    """
+    Wrap the call of SurfStat using clinicasurfstat.m Matlab script.
 
     Args:
-        design_matrix: str, the linear model that fits into the GLM, for example '1+group'.
-        contrast: string, the contrast matrix for GLM, if the factor you choose is categorized variable, clinica_surfstat will create two contrasts,
-                  for example, contrast = 'Label', this will create contrastpos = Label.AD - Label.CN, contrastneg = Label.CN - Label.AD; if the fac-
-                  tory that you choose is a continuous factor, clinica_surfstat will just create one contrast, for example, contrast = 'Age', but note,
-                  the string name that you choose should be exactly the same with the columns names in your subjects_visits_tsv.
-        str_format:string, the str_format which uses to read your tsv file, the typy of the string should corresponds exactly with the columns in the tsv file.
-                  Defaut parameters, we set these parameters to be some default values, but you can also set it by yourself:
-        glm_type: based on the hypothesis, you should define one of the glm types, "group_comparison", "correlation"
-        group_label: current group name for this analysis
-        full_width_at_half_maximum: fwhm for the surface smoothing, default is 20, integer.
-        threshold_uncorrected_pvalue: threshold to display the uncorrected Pvalue, float, default is 0.001.
-        threshold_corrected_pvalue: the threshold to display the corrected cluster, default is 0.05, float.
-        cluster_threshold: threshold to define a cluster in the process of cluster-wise correction, default is 0.001, float.
+        caps_dir (str): CAPS directory containing surface-based features
+        output_dir (str): Output directory that will contain outputs of clinicasurfstat.m
+        subjects_visits_tsv (str): TSV file containing the GLM information
+        pipeline_parameters (dict): parameters of StatisticsSurface pipeline
+    """
+    import os
+    import sys
+    from nipype.interfaces.matlab import MatlabCommand, get_matlab_command
+    import clinica.pipelines as clinica_pipelines
+    from clinica.utils.check_dependency import check_environment_variable
+    from clinica.pipelines.statistics_surface.statistics_surface_utils import covariates_to_design_matrix, get_string_format_from_tsv
 
-    Returns:
+    path_to_matlab_script = os.path.join(os.path.dirname(clinica_pipelines.__path__[0]), 'lib', 'clinicasurfstat')
+    freesurfer_home = check_environment_variable('FREESURFER_HOME', 'FreeSurfer')
 
-    """
-    json_dict = {
-        'AnalysisType': glm_type,
-        'DesignMatrix': design_matrix,
-        'StringFormatTSV': str_format,
-        'Contrast': contrast,
-        'GroupLabel': group_label,
-        'FWHM': full_width_at_half_maximum,
-        'ThresholdUncorrectedPvalue': threshold_uncorrected_pvalue,
-        'ThresholdCorrectedPvalue': threshold_corrected_pvalue,
-        'ClusterThreshold': cluster_threshold
+    MatlabCommand.set_default_matlab_cmd(
+        get_matlab_command()
+    )
+    matlab = MatlabCommand()
+    matlab.inputs.paths = path_to_matlab_script
+    matlab.inputs.script = """
+    clinicasurfstat('%s', '%s', '%s', '%s', '%s', '%s', '%s', '%s', '%s', '%s', '%s', '%s', %d, '%s', %.3f, '%s', %.3f, '%s', %.3f);
+    """ % (
+        os.path.join(caps_dir, 'subjects'),
+        output_dir,
+        subjects_visits_tsv,
+        covariates_to_design_matrix(
+            pipeline_parameters['contrast'],
+            pipeline_parameters['covariates']
+        ),
+        pipeline_parameters['contrast'],
+        get_string_format_from_tsv(subjects_visits_tsv),
+        pipeline_parameters['glm_type'],
+        pipeline_parameters['group_label'],
+        freesurfer_home,
+        pipeline_parameters['custom_file'],
+        pipeline_parameters['measure_label'],
+        'sizeoffwhm', pipeline_parameters['full_width_at_half_maximum'],
+        'thresholduncorrectedpvalue', 0.001,
+        'thresholdcorrectedpvalue', 0.05,
+        'clusterthreshold', pipeline_parameters['cluster_threshold']
+        )
+    # This will create a file: pyscript.m , the pyscript.m is the default name
+    matlab.inputs.mfile = True
+    # This will stop running with single thread
+    matlab.inputs.single_comp_thread = False
+    matlab.inputs.logfile = 'group-' + pipeline_parameters['group_label'] + '_matlab.log'
+
+    # cprint("Matlab logfile is located at the following path: %s" % matlab.inputs.logfile)
+    # cprint("Matlab script command = %s" % matlab.inputs.script)
+    # cprint("MatlabCommand inputs flag: single_comp_thread = %s" % matlab.inputs.single_comp_thread)
+    # cprint("MatlabCommand choose which matlab to use(matlab_cmd): %s" % get_matlab_command())
+    matlab.run()
+
+    return output_dir
+
+
+def create_glm_info_dictionary(tsv_file, pipeline_parameters):
+    """Create dictionary containing the GLM information that will be stored in a JSON file."""
+    out_dict = {
+        'AnalysisType': pipeline_parameters['glm_type'],
+        'DesignMatrix':  covariates_to_design_matrix(
+            pipeline_parameters['contrast'],
+            pipeline_parameters['covariates']
+        ),
+        'StringFormatTSV': get_string_format_from_tsv(tsv_file),
+        'Contrast': pipeline_parameters['contrast'],
+        'Covariates': pipeline_parameters['covariates'],
+        'GroupLabel': pipeline_parameters['group_label'],
+        'FWHM': pipeline_parameters['full_width_at_half_maximum'],
+        'ThresholdUncorrectedPvalue': 0.001,
+        'ThresholdCorrectedPvalue': 0.05,
+        'ClusterThreshold': pipeline_parameters['cluster_threshold']
     }
+    return out_dict
 
-    return json_dict
 
+def save_to_caps(source_dir, caps_dir, overwrite_caps, pipeline_parameters):
+    """Save `source_dir`/ to CAPS folder.
 
-def check_inputs(caps,
-                 custom_filename,
-                 fwhm,
-                 tsv_file):
-    """
-        Simply checks if the custom strings provided find the correct files. If files are missing, it is easier to
-        raise an exception here rather than doing it in the MATLAB script.
-    Args:
-        caps: The CAPS folder
-        custom_filename: The string defining where are the files needed for the surfstat analysis
-        fwhm: Full Width at Half Maximum used for the smoothing of the data, needed to catch the file
-        tsv_file: tsv file that contains the subject and session lists, along with the other covariates
+    This function copies outputs of `source_dir`/ to
+    `caps_dir`/groups/<group_id>/<statistics>/surfstat_<glm_type>/
 
-    Returns:
+    The `source_dir`/ folder should contain the following elements:
+        - group-<group_label>_<group_1_or_2>-lt-<group_1_or_2>_measure-<measure>_fwhm-<label>_suffix.ext
+    or
+        - group-<group_label>_correlation-<label>_contrast-{-|+}_measure-<measure>_fwhm-<label>_suffix.ext
+    and
+        - group-<group_label>_covariates.tsv
+        - group-<group_label>_glm.json
 
+    Raise:
+        NotImplementedError: If overwrite_caps=True.
     """
     import os
-    from clinica.utils.stream import cprint
-    import pandas as pd
+    import shutil
+    from clinica.utils.ux import print_end_image
+
+    group_id = 'group-' + pipeline_parameters['group_label']
+
+    if pipeline_parameters['glm_type'] == "group_comparison":
+        surfstat_folder = 'surfstat_' + pipeline_parameters['glm_type']
+    elif pipeline_parameters['glm_type'] == "correlation":
+        surfstat_folder = 'surfstat_' + pipeline_parameters['glm_type'] + '_analysis'
+    else:
+        raise NotImplementedError("The other GLM situations have not been implemented in this pipeline.")
+
+    destination_dir = os.path.join(
+        os.path.expanduser(caps_dir),
+        'groups',
+        group_id,
+        'statistics',
+        surfstat_folder
+    )
 
-    subjects_visits = pd.io.parsers.read_csv(tsv_file, sep='\t')
-    subjects = list(subjects_visits.participant_id)
-    sessions = list(subjects_visits.session_id)
-
-    missing_files = []
-    for idx in range(len(subjects)):
-        fullpath = os.path.join(caps,
-                                'subjects',
-                                custom_filename.replace('@subject', subjects[idx]).replace('@session', sessions[idx]).replace('@fwhm', str(fwhm)))
-        left_hemi = fullpath.replace('@hemi', 'lh')
-        right_hemi = fullpath.replace('@hemi', 'rh')
-
-        if not os.path.exists(left_hemi):
-            missing_files.append(left_hemi)
-        if not os.path.exists(right_hemi):
-            missing_files.append(right_hemi)
-
-    if len(missing_files) > 0:
-        cprint(' ** Missing files **')
-        for l in missing_files:
-            cprint('Not found: ' + l)
-        raise Exception(str(len(missing_files)) + ' files not found !')
+    if overwrite_caps:
+        raise NotImplementedError('save_to_caps(overwrite_caps=True) not implemented')
+    shutil.copytree(source_dir, destination_dir, symlinks=True)
+    print_end_image(group_id)
diff --git a/clinica/pipelines/statistics_volume/statistics_volume_cli.py b/clinica/pipelines/statistics_volume/statistics_volume_cli.py
index 02d5f32c..e97db186 100644
--- a/clinica/pipelines/statistics_volume/statistics_volume_cli.py
+++ b/clinica/pipelines/statistics_volume/statistics_volume_cli.py
@@ -17,35 +17,53 @@ class StatisticsVolumeCLI(ce.CmdParser):
     def define_options(self):
         """Define the sub-command arguments."""
         from clinica.engine.cmdparser import PIPELINE_CATEGORIES
+        from colorama import Fore
 
         # Clinica compulsory arguments
         clinica_comp = self._args.add_argument_group(PIPELINE_CATEGORIES['CLINICA_COMPULSORY'])
         clinica_comp.add_argument("caps_directory",
                                   help='Path to the CAPS directory.')
+        clinica_comp.add_argument("group_label",
+                                  help='User-defined identifier for the provided group of subjects.')
+        clinica_comp.add_argument("orig_input_data",
+                                  help='Type of volume-based feature: type '
+                                  '\'t1-volume\' to use gray matter maps, '
+                                  '\'pet-volume\' to use FDG-PET data or '
+                                  '\'custom-pipeline\' to use you own data in CAPS directory '
+                                  '(see Wiki for details).',
+                                  choices=['t1-volume', 'pet-volume', 'custom-pipeline'])
         clinica_comp.add_argument("subject_visits_with_covariates_tsv",
                                   help='TSV file containing a list of subjects with their sessions and all '
                                        'the covariates and factors needed for the GLM.')
         clinica_comp.add_argument("contrast",
                                   help='Defines the contrast. Must be one of the column names form the TSV file.')
-        clinica_comp.add_argument("feature_type",
-                                  help='Type of volume-based feature: graymatter (from t1-volume pipeline) or '
-                                       'fdg (from pet-volume pipeline). Use your own if you want to use '
-                                       'the --custom_file flag. ')
-        clinica_comp.add_argument("group_id",
-                                  help='User-defined identifier for the provided group of subjects.')
 
         # Optional arguments (e.g. FWHM)
         optional = self._args.add_argument_group(PIPELINE_CATEGORIES['OPTIONAL'])
-        optional.add_argument("--custom_files", "-cf", type=str, default=None,
-                              help='Custom file string. Specify filename using * when the subject or session name '
-                                   'appears e.g. \'*_task-rest_acq-fdg_pet_space-Ixi549Space_pet.nii.gz\' will grab '
-                                   'the corresponding file in all the subjects/sessions.')
         optional.add_argument("-gic", "--group_id_caps", type=str, default=None,
                               help='Name of the group that Clinica needs to use to grab input file.')
         optional.add_argument("-fwhm", "--full_width_at_half_maximum", type=int, default=8,
                               help='Full Width at Half Maximum (FWHM) of the smoothing used in your input file '
                                    '(default: --full_width_at_half_maximum %(default)s).')
 
+        # Optional arguments for custom pipeline
+        opt_custom_input = self._args.add_argument_group(
+            '%sPipeline options if you selected custom-pipeline%s' % (Fore.BLUE, Fore.RESET)
+        )
+        opt_custom_input.add_argument("-cf", "--custom_file",
+                                      type=str, default=None,
+                                      help='Custom file string. Specify filename using * when the subject or session name '
+                                      'appears e.g. \'*_task-rest_acq-fdg_pet_space-Ixi549Space_pet.nii.gz\' will grab '
+                                      'the corresponding file in all the subjects/sessions. '
+                                      'This flag must be specified with the --measure_label flag). '
+                                      'See Wiki for an example.')
+        opt_custom_input.add_argument("-ml", "--measure_label",
+                                      type=str, default=None,
+                                      help='Name of the feature type, it will be saved on the CAPS '
+                                           '_measure-FEATURE_LABEL key-value association. '
+                                           'This flag must be specified with the --custom_file flag). '
+                                           'See Wiki for an example.')
+
         # Clinica standard arguments (e.g. --n_procs)
         self.add_clinica_standard_arguments(add_tsv_flag=False)
 
@@ -60,11 +78,17 @@ class StatisticsVolumeCLI(ce.CmdParser):
         from networkx import Graph
         from .statistics_volume_pipeline import StatisticsVolume
         from clinica.utils.ux import print_end_pipeline, print_crash_files_and_exit
+        from clinica.utils.exceptions import ClinicaException
+
+        # Custom pipeline
+        if args.orig_input_data == 'custom-pipeline':
+            if (args.custom_file is None) or (args.measure_label is None):
+                raise ClinicaException('You must set --measure_label and --custom_file flags.')
 
         parameters = {
             'contrast': args.contrast,
-            'feature_type': args.feature_type,
-            'group_id': args.group_id,
+            'orig_input_data': args.orig_input_data,
+            'group_label': args.group_id,
             'custom_files': args.custom_files,
             'cluster_threshold': args.cluster_threshold,
             'group_id_caps': args.group_id_caps,
diff --git a/clinica/pipelines/statistics_volume/statistics_volume_pipeline.py b/clinica/pipelines/statistics_volume/statistics_volume_pipeline.py
index feef9f05..2090df50 100644
--- a/clinica/pipelines/statistics_volume/statistics_volume_pipeline.py
+++ b/clinica/pipelines/statistics_volume/statistics_volume_pipeline.py
@@ -13,6 +13,9 @@ class StatisticsVolume(cpe.Pipeline):
         """Check pipeline parameters."""
         from clinica.utils.exceptions import ClinicaException
 
+        if 'orig_input_data' not in self.parameters.keys():
+            raise KeyError('Missing compulsory orig_input_data key in pipeline parameter.')
+
         if self.parameters['cluster_threshold'] < 0 or self.parameters['cluster_threshold'] > 1:
             raise ClinicaException("Cluster threshold should be between 0 and 1 "
                                    "(given value: %s)." % self.parameters['cluster_threshold'])
@@ -48,15 +51,12 @@ class StatisticsVolume(cpe.Pipeline):
 
     def build_input_node(self):
         """Build and connect an input node to the pipeline."""
-        import os
         from colorama import Fore
         import nipype.interfaces.utility as nutil
         import nipype.pipeline.engine as npe
-        from clinica.iotools.utils.data_handling import check_volume_location_in_world_coordinate_system
-        from clinica.utils.exceptions import ClinicaBIDSError, ClinicaException
-        from clinica.utils.filemanip import extract_subjects_sessions_from_filename, save_participants_sessions
+        from clinica.utils.exceptions import ClinicaException
         from clinica.utils.inputs import clinica_file_reader
-        from clinica.utils.input_files import T1W_NII
+        from clinica.utils.input_files import t1_volume_template_tpm_in_mni
         from clinica.utils.stream import cprint
         from clinica.utils.ux import print_images_to_process, print_begin_image
 
@@ -65,46 +65,40 @@ class StatisticsVolume(cpe.Pipeline):
             gic = self.parameters['group_id_caps']
 
         all_errors = []
-        if self.parameters['feature_type'] == 'fdg':
-            try:
-                input_files = clinica_file_reader(self.subjects,
-                                                  self.sessions,
-                                                  self.caps_directory,
-                                                  {'pattern': '*_pet_space-Ixi549Space_suvr-pons_mask-brain_fwhm-' + str(self.parameters['full_width_at_half_maximum']) + 'mm_pet.nii*',
-                                                   'description': 'pons normalized FDG PET image in MNI space (brain masked)',
-                                                   'needed_pipeline': 'pet-volume'})
-            except ClinicaException as e:
-                all_errors.append(e)
-        elif self.parameters['feature_type'] == 'graymatter':
-            try:
-                input_files = clinica_file_reader(self.subjects,
-                                                  self.sessions,
-                                                  self.caps_directory,
-                                                  {'pattern': 't1/spm/dartel/group-' + gic + '/*_T1w_segm-graymatter_space-Ixi549Space_modulated-on_fwhm-' + str(self.parameters['full_width_at_half_maximum']) + 'mm_probability.nii.*',
-                                                   'description': 'probability map of gray matter segmentation based on T1w image in MNI space',
-                                                   'needed_pipeline': 't1-volume or t1-volume-existing-template'})
-            except ClinicaException as e:
-                all_errors.append(e)
+        if self.parameters['orig_input_data'] == 'pet-volume':
+            self.parameters['measure_label'] = 'fdg'
+            information_dict = {
+                'pattern': '*_pet_space-Ixi549Space_suvr-pons_mask-brain_fwhm-' + str(self.parameters['full_width_at_half_maximum']) + 'mm_pet.nii*',
+                'description': 'pons normalized FDG PET image in MNI space (brain masked)',
+                'needed_pipeline': 'pet-volume'
+            }
+        elif self.parameters['orig_input_data'] == 't1-volume':
+            self.parameters['measure_label'] = 'graymatter'
+            information_dict = t1_volume_template_tpm_in_mni(gic, 0, True)
 
         else:
-            if not self.parameters['custom_files']:
-                raise ClinicaException(Fore.RED + '[Error] You did not specify the --custom_files flag in the command line for the feature type '
-                                       + Fore.Blue + self.parameters['feature_type'] + Fore.RED + '! Clinica can\'t '
+            if not self.parameters['custom_file']:
+                raise ClinicaException(Fore.RED + '[Error] You did not specify the --custom_file flag in the command line for the feature type '
+                                       + Fore.Blue + self.parameters['measure_label'] + Fore.RED + '! Clinica can\'t '
                                        + 'know what file to use in your analysis ! Type: \n\t' + Fore.BLUE + 'clinica run statistics-volume\n'
                                        + Fore.RED + ' to have help on how to use the command line.' + Fore.RESET)
-            try:
-                # If custom file are grabbed, information of fwhm is irrelevant and should not appear on final filenames
-                self.parameters['full_width_at_half_maximum'] = None
-                input_files = clinica_file_reader(self.subjects,
-                                                  self.sessions,
-                                                  self.caps_directory,
-                                                  {'pattern': self.parameters['custom_files'],
-                                                   'description': 'custom file provided by user'})
-            except ClinicaException as e:
-                all_errors.append(e)
+            # If custom file are grabbed, information of fwhm is irrelevant and should not appear on final filenames
+            self.parameters['full_width_at_half_maximum'] = None
+            information_dict = {
+                'pattern': self.parameters['custom_file'],
+                'description': 'custom file provided by user'
+            }
+
+        try:
+            input_files = clinica_file_reader(self.subjects,
+                                              self.sessions,
+                                              self.caps_directory,
+                                              information_dict)
+        except ClinicaException as e:
+            all_errors.append(e)
 
         if len(all_errors) > 0:
-            error_message = 'Clinica faced errors while trying to read files in your BIDS or CAPS directories.\n'
+            error_message = 'Clinica faced errors while trying to read files in your CAPS directories.\n'
             for msg in all_errors:
                 error_message += str(msg)
             raise ClinicaException(error_message)
@@ -119,7 +113,7 @@ class StatisticsVolume(cpe.Pipeline):
         if len(self.subjects):
             print_images_to_process(self.subjects, self.sessions)
             cprint('The pipeline will last a few minutes. Images generated by SPM will popup during the pipeline.')
-            print_begin_image('group-' + self.parameters['group_id'])
+            print_begin_image('group-' + self.parameters['group_label'])
 
         self.connect([
             (read_parameters_node,      self.input_node,    [('input_files',    'input_files')])
@@ -132,9 +126,9 @@ class StatisticsVolume(cpe.Pipeline):
         from os.path import join, pardir
 
         relative_path = join('groups',
-                             'group-' + self.parameters['group_id'],
+                             'group-' + self.parameters['group_label'],
                              'statistics_volume',
-                             'group_comparison_measure-' + self.parameters['feature_type'])
+                             'group_comparison_measure-' + self.parameters['measure_label'])
 
         datasink = npe.Node(nio.DataSink(), name='sinker')
         datasink.inputs.base_directory = join(self.caps_directory, relative_path)
@@ -152,11 +146,11 @@ class StatisticsVolume(cpe.Pipeline):
 
                 # resels per voxels
                 (join(self.caps_directory, relative_path) + '/resels_per_voxels/resels_per_voxel.nii',
-                 join(self.caps_directory, relative_path) + '/group-' + self.parameters['group_id'] + '_RPV.nii'),
+                 join(self.caps_directory, relative_path) + '/group-' + self.parameters['group_label'] + '_RPV.nii'),
 
                 # mask
                 (join(self.caps_directory, relative_path) + '/mask/included_voxel_mask.nii',
-                 join(self.caps_directory, relative_path) + '/group-' + self.parameters['group_id'] + '_mask.nii'),
+                 join(self.caps_directory, relative_path) + '/group-' + self.parameters['group_label'] + '_mask.nii'),
 
                 # variance of error
                 (join(self.caps_directory, relative_path) + '/variance_of_error/(.*)',
@@ -164,7 +158,7 @@ class StatisticsVolume(cpe.Pipeline):
 
                 # tsv file
                 (join(self.caps_directory, relative_path) + r'/tsv_file/.*',
-                 join(self.caps_directory, relative_path) + '/' + pardir + '/group-' + self.parameters['group_id'] + '_participants.tsv'),
+                 join(self.caps_directory, relative_path) + '/' + pardir + '/group-' + self.parameters['group_label'] + '_participants.tsv'),
 
                 # report (figures)
                 (join(self.caps_directory, relative_path) + r'/figures/(.*)',
@@ -172,8 +166,8 @@ class StatisticsVolume(cpe.Pipeline):
 
                 # regression coefficient
                 (join(self.caps_directory, relative_path) + r'/regression_coeff/(.*).nii',
-                 join(self.caps_directory, relative_path) + '/group-' + self.parameters['group_id'] + r'_covariate-\1'
-                 + '_measure-' + self.parameters['feature_type'] + '_fwhm-'
+                 join(self.caps_directory, relative_path) + '/group-' + self.parameters['group_label'] + r'_covariate-\1'
+                 + '_measure-' + self.parameters['measure_label'] + '_fwhm-'
                  + str(self.parameters['full_width_at_half_maximum']) + '_regressionCoefficient.nii'),
 
             ]
@@ -278,20 +272,15 @@ class StatisticsVolume(cpe.Pipeline):
         # Print result to txt file if spm
 
         # Export results to output node
-        read_output_node = npe.Node(nutil.Function(input_names=['spm_mat', 'class_names', 'covariates', 'group_id', 'fwhm', 'measure'],
-                                                   output_names=['spmT_0001',
-                                                                 'spmT_0002',
-                                                                 'spm_figures',
-                                                                 'variance_of_error',
-                                                                 'resels_per_voxels',
-                                                                 'mask',
-                                                                 'regression_coeff',
-                                                                 'contrasts'],
-                                                   function=utils.read_output),
-                                    name='read_output_node')
-        read_output_node.inputs.group_id = self.parameters['group_id']
+        read_output_node = npe.Node(nutil.Function(
+            input_names=['spm_mat', 'class_names', 'covariates', 'group_label', 'fwhm', 'measure'],
+            output_names=['spmT_0001', 'spmT_0002', 'spm_figures', 'variance_of_error',
+                          'resels_per_voxels', 'mask', 'regression_coeff', 'contrasts'],
+            function=utils.read_output),
+            name='read_output_node')
+        read_output_node.inputs.group_label = self.parameters['group_label']
         read_output_node.inputs.fwhm = self.parameters['full_width_at_half_maximum']
-        read_output_node.inputs.measure = self.parameters['feature_type']
+        read_output_node.inputs.measure = self.parameters['measure_label']
 
         # Connection
         # ==========
diff --git a/clinica/pipelines/statistics_volume/statistics_volume_utils.py b/clinica/pipelines/statistics_volume/statistics_volume_utils.py
index 632245d2..366c9fc8 100644
--- a/clinica/pipelines/statistics_volume/statistics_volume_utils.py
+++ b/clinica/pipelines/statistics_volume/statistics_volume_utils.py
@@ -360,7 +360,7 @@ def contrast(mat_file, template_file, covariates, class_names):
     return current_model_estimation
 
 
-def read_output(spm_mat, class_names, covariates, group_id, fwhm, measure):
+def read_output(spm_mat, class_names, covariates, group_label, fwhm, measure):
     """
         Once analysis is done, grab all the different filenames and rename them in current directory according to class
         names
@@ -368,7 +368,7 @@ def read_output(spm_mat, class_names, covariates, group_id, fwhm, measure):
         spm_mat: (str) path to the SPM.mat file of the SPM analysis
         class_names: (list) of str of length 2 that correspond to the 2 classes for the group comparison
         covariates: (list) of str: list of covariates
-        group_id: name of the group id
+        group_label: name of the group label
         fwhm: fwhm in mm used
         measure: measure used
 
@@ -400,7 +400,7 @@ def read_output(spm_mat, class_names, covariates, group_id, fwhm, measure):
     if len(figures) < 2:
         raise RuntimeError('[Error] Figures were not generated')
     fig_number = [int(f[-7:-4]) for f in figures]
-    spm_figures = [abspath('./' + 'group-' + group_id + '_report-' + str(i) + '.png') for i in fig_number]
+    spm_figures = [abspath('./' + 'group-' + group_label + '_report-' + str(i) + '.png') for i in fig_number]
     for old_name, new_name in zip(figures, spm_figures):
         copyfile(old_name, new_name)
 
@@ -410,15 +410,15 @@ def read_output(spm_mat, class_names, covariates, group_id, fwhm, measure):
     if len(spm_T) != 2:
         raise RuntimeError('[Error] ' + str(len(spm_T)) + ' SPM t-map(s) were found')
     if fwhm:
-        spmT_0001 = abspath('group-' + group_id + '_' + class_names[0] + '-lt-' + class_names[1] + '_measure-' + measure + '_fwhm-' + str(int(fwhm)) + '_TStatistics.nii')
-        spmT_0002 = abspath('group-' + group_id + '_' + class_names[1] + '-lt-' + class_names[0] + '_measure-' + measure + '_fwhm-' + str(int(fwhm)) + '_TStatistics.nii')
+        spmT_0001 = abspath('group-' + group_label + '_' + class_names[0] + '-lt-' + class_names[1] + '_measure-' + measure + '_fwhm-' + str(int(fwhm)) + '_TStatistics.nii')
+        spmT_0002 = abspath('group-' + group_label + '_' + class_names[1] + '-lt-' + class_names[0] + '_measure-' + measure + '_fwhm-' + str(int(fwhm)) + '_TStatistics.nii')
     else:
-        spmT_0001 = abspath('group-' + group_id + '_' + class_names[0] + '-lt-' + class_names[1] + '_measure-' + measure + '_TStatistics.nii')
-        spmT_0002 = abspath('group-' + group_id + '_' + class_names[1] + '-lt-' + class_names[0] + '_measure-' + measure + '_TStatistics.nii')
+        spmT_0001 = abspath('group-' + group_label + '_' + class_names[0] + '-lt-' + class_names[1] + '_measure-' + measure + '_TStatistics.nii')
+        spmT_0002 = abspath('group-' + group_label + '_' + class_names[1] + '-lt-' + class_names[0] + '_measure-' + measure + '_TStatistics.nii')
     copyfile(join(dirname(spm_mat), 'spmT_0001.nii'), spmT_0001)
     copyfile(join(dirname(spm_mat), 'spmT_0002.nii'), spmT_0002)
 
-    variance_of_error = abspath('./group-' + group_id + '_VarianceError.nii')
+    variance_of_error = abspath('./group-' + group_label + '_VarianceError.nii')
     copyfile(abspath(join(dirname(spm_mat), 'ResMS.nii')), variance_of_error)
 
     resels_per_voxels = abspath('./resels_per_voxel.nii')
@@ -444,8 +444,8 @@ def read_output(spm_mat, class_names, covariates, group_id, fwhm, measure):
     con_files = [abspath(join(dirname(spm_mat), f)) for f in list_files if f.startswith('con_')]
     if len(con_files) != 2:
         raise RuntimeError('There must exists only 2 contrast files !')
-    contrasts = [abspath('group-' + group_id + '_' + class_names[0] + '-lt-' + class_names[1] + '_measure-' + measure + '_contrast.nii'),
-                 abspath('group-' + group_id + '_' + class_names[1] + '-lt-' + class_names[0] + '_measure-' + measure + '_contrast.nii')]
+    contrasts = [abspath('group-' + group_label + '_' + class_names[0] + '-lt-' + class_names[1] + '_measure-' + measure + '_contrast.nii'),
+                 abspath('group-' + group_label + '_' + class_names[1] + '-lt-' + class_names[0] + '_measure-' + measure + '_contrast.nii')]
     for con, contrast in zip(con_files, contrasts):
         copyfile(con, contrast)
 
diff --git a/clinica/utils/inputs.py b/clinica/utils/inputs.py
index 8bfb3bb7..cd84ec0c 100644
--- a/clinica/utils/inputs.py
+++ b/clinica/utils/inputs.py
@@ -131,7 +131,7 @@ def clinica_file_reader(subjects,
         subjects: list of subjects
         sessions: list of sessions (must be same size as subjects, and must correspond )
         input_directory: location of the bids or caps directory
-        information: dictionnary containg all the relevant information to look for the files. Dict must contains the
+        information: dictionary containing all the relevant information to look for the files. Dict must contains the
                      following keys : pattern, description. The optional key is: needed_pipeline
                              pattern: define the pattern of the final file
                              description: string to describe what the file is
diff --git a/clinica/utils/participant.py b/clinica/utils/participant.py
index 4656521e..b16f75d6 100644
--- a/clinica/utils/participant.py
+++ b/clinica/utils/participant.py
@@ -108,3 +108,25 @@ def get_subject_session_list(input_dir, ss_file=None, is_bids_dir=True, use_sess
 
     participant_ids, session_ids = read_participant_tsv(ss_file)
     return session_ids, participant_ids
+
+
+def have_same_subjects(tsv_file_1, tsv_file_2):
+    """Return True if `tsv_file_1` and `tsv_file_2` have the same subjects, False otherwise."""
+    import pandas as pd
+
+    tsv_df_1 = pd.io.parsers.read_csv(tsv_file_1, sep='\t')
+    tsv_df_2 = pd.io.parsers.read_csv(tsv_file_2, sep='\t')
+    image_ids_1 = [
+        p_id + '_' + s_id
+        for (p_id, s_id) in zip(list(tsv_df_1.participant_id), list(tsv_df_1.session_id))
+    ]
+    image_ids_2 = [
+        p_id + '_' + s_id
+        for (p_id, s_id) in zip(list(tsv_df_2.participant_id), list(tsv_df_2.session_id))
+    ]
+    diff_image_ids = list(set(image_ids_1) - set(image_ids_2))
+
+    if len(diff_image_ids) == 0:
+        return True
+    else:
+        return False
diff --git a/test/instantiation/test_instantiate_all_pipelines.py b/test/instantiation/test_instantiate_all_pipelines.py
index c2a473d2..677a3309 100644
--- a/test/instantiation/test_instantiate_all_pipelines.py
+++ b/test/instantiation/test_instantiate_all_pipelines.py
@@ -220,16 +220,14 @@ def test_instantiate_StatisticsSurface():
     root = dirname(abspath(join(abspath(__file__), pardir)))
     root = join(root, 'data', 'StatisticsSurface')
     parameters = {
-        'design_matrix': '1 + group + age + sex',
+        'orig_input_data': 't1-freesurfer',
+        'covariates': 'age + sex',
         'contrast': 'group',
-        'str_format': '%s %s %s %f %s',
         'group_label': 'UnitTest',
         'glm_type': 'group_comparison',
         'custom_file': '@subject/@session/t1/freesurfer_cross_sectional/@subject_@session/surf/@hemi.thickness.fwhm@fwhm.fsaverage.mgh',
-        'feature_label': 'cortical_thickness',
+        'measure_label': 'ct',
         'full_width_at_half_maximum': 20,
-        'threshold_uncorrected_pvalue': 0.001,
-        'threshold_corrected_pvalue': 0.05,
         'cluster_threshold': 0.001
     }
     pipeline = StatisticsSurface(
@@ -360,9 +358,10 @@ def test_instantiate_StatisticsVolume():
 
     # Instantiate pipeline and run()
     parameters = {
+        'orig_input_data': 'pet-volume',
         'contrast': 'group',
-        'feature_type': 'fdg',
-        'group_id': 'UnitTest',
+        'measure_label': 'fdg',
+        'group_label': 'UnitTest',
         'cluster_threshold': 0.001,
         'group_id_caps': None,
         'full_width_at_half_maximum': 8
diff --git a/test/nonregression/test_run_pipelines.py b/test/nonregression/test_run_pipelines.py
index 379436d7..d17317f4 100644
--- a/test/nonregression/test_run_pipelines.py
+++ b/test/nonregression/test_run_pipelines.py
@@ -527,16 +527,14 @@ def test_run_StatisticsSurface(cmdopt):
     shutil.copytree(join(root, 'in', 'caps'), join(root, 'out', 'caps'))
 
     parameters = {
-        'design_matrix': '1 + group + age + sex',
+        'orig_input_data': 't1-freesurfer',
+        'covariates': 'age sex',
         'contrast': 'group',
-        'str_format': '%s %s %s %f %s',
         'group_label': 'UnitTest',
         'glm_type': 'group_comparison',
         'custom_file': '@subject/@session/t1/freesurfer_cross_sectional/@subject_@session/surf/@hemi.thickness.fwhm@fwhm.fsaverage.mgh',
-        'feature_label': 'cortical_thickness',
+        'measure_label': 'ct',
         'full_width_at_half_maximum': 20,
-        'threshold_uncorrected_pvalue': 0.001,
-        'threshold_corrected_pvalue': 0.05,
         'cluster_threshold': 0.001
     }
     pipeline = StatisticsSurface(
@@ -549,8 +547,9 @@ def test_run_StatisticsSurface(cmdopt):
     pipeline.run(plugin='MultiProc', plugin_args={'n_procs': 8}, bypass_check=True)
 
     # Check files
-    out_file = join(root, 'out/caps/groups/group-UnitTest/statistics/surfstat_group_comparison/group-UnitTest_AD-lt-CN_measure-cortical_thickness_fwhm-20_correctedPValue.mat')
-    ref_file = join(root, 'ref/group-UnitTest_AD-lt-CN_measure-cortical_thickness_fwhm-20_correctedPValue.mat')
+    filename = 'group-UnitTest_AD-lt-CN_measure-ct_fwhm-20_correctedPValue.mat'
+    out_file = join(root, 'out', 'caps', 'groups', 'group-UnitTest', 'statistics', 'surfstat_group_comparison', filename)
+    ref_file = join(root, 'ref', filename)
 
     out_file_mat = loadmat(out_file)['correctedpvaluesstruct']
     ref_file_mat = loadmat(ref_file)['correctedpvaluesstruct']
@@ -720,9 +719,10 @@ def test_run_StatisticsVolume(cmdopt):
 
     # Instantiate pipeline and run()
     parameters = {
+        'orig_input_data': 'pet-volume',
         'contrast': 'group',
-        'feature_type': 'fdg',
-        'group_id': 'UnitTest',
+        'measure_label': 'fdg',
+        'group_label': 'UnitTest',
         'cluster_threshold': 0.001,
         'group_id_caps': None,
         'full_width_at_half_maximum': 8
['clinica/pipelines/statistics_volume/statistics_volume_pipeline.py', 'clinica/utils/participant.py', 'clinica/pipelines/statistics_surface/statistics_surface_pipeline.py', 'test/instantiation/test_instantiate_all_pipelines.py', 'test/nonregression/test_run_pipelines.py', 'clinica/pipelines/statistics_surface/info.json', 'clinica/pipelines/statistics_surface/statistics_surface_cli.py', 'clinica/pipelines/statistics_surface/statistics_surface_utils.py', 'clinica/pipelines/statistics_volume/statistics_volume_cli.py', 'clinica/utils/inputs.py', 'clinica/pipelines/statistics_volume/statistics_volume_utils.py', 'clinica/lib/clinicasurfstat/clinicasurfstat.m']