##compute the number of ADDITIONAL threads to be used by samtools (-@)
        addthreads=\${GALAXY_SLOTS:-1} && (( addthreads-- )) &&
    
        
        ##prepare input and indices
        ln -s '$input' infile &&
        #if $input.is_of_type('bam'):
            #if str( $input.metadata.bam_index ) != "None":
                ln -s '${input.metadata.bam_index}' infile.bai &&
            #else:
                samtools index infile infile.bai &&
            #end if
        #elif $input.is_of_type('cram'):
            #if str( $input.metadata.cram_index ) != "None":
                ln -s '${input.metadata.cram_index}' infile.crai &&
            #else:
                samtools index infile infile.crai &&
            #end if
        #end if
    
        
        ## Make the user-selected reference genome, if any, accessible through
        ## a shell variable $reffa, index the reference if necessary, and make
        ## the fai-index file available through a shell variable $reffai.

        ## For a cached genome simply sets the shell variables to point to the
        ## genome file and its precalculated index.
        ## For a genome from the user's history, if that genome is a plain
        ## fasta file, the code creates a symlink in the pwd, creates the fai
        ## index file next to it, then sets the shell variables to point to the
        ## symlink and its index.
        ## For a fasta.gz dataset from the user's history, it tries the same,
        ## but this will only succeed if the file got compressed with bgzip.
        ## For a regular gzipped file samtools faidx will fail, in which case
        ## the code falls back to decompressing to plain fasta before
        ## reattempting the indexing.
        ## Indexing of a bgzipped file produces a regular fai index file *and*
        ## a compressed gzi file. The former is identical to the fai index of
        ## the uncompressed fasta.

        ## If the user has not selected a reference (it's an optional parameter
        ## in some samtools wrappers), a cheetah boolean use_ref is set to
        ## False to encode that fact.

        #set use_ref=True
        #if $addref_cond.addref_select == "history":
            #if $addref_cond.ref.is_of_type('fasta'):
                reffa="reference.fa" &&
                ln -s '${addref_cond.ref}' \$reffa &&
                samtools faidx \$reffa &&
            #else:
                reffa="reference.fa.gz" &&
                ln -s '${addref_cond.ref}' \$reffa &&
                {
                    samtools faidx \$reffa ||
                    {
                        echo "Failed to index compressed reference. Trying decompressed ..." 1>&2 &&
                        gzip -dc \$reffa > reference.fa &&
                        reffa="reference.fa" &&
                        samtools faidx \$reffa;
                    }
                } &&
            #end if
            reffai=\$reffa.fai &&
        #elif $addref_cond.addref_select == "cached":
            ## in case of cached the absolute path is used which allows to read 
            ## a cram file  without specifying the reference
            reffa='${addref_cond.ref.fields.path}' &&
            reffai=\$reffa.fai &&
        #else
            #set use_ref=False
        #end if
    
        samtools stats
            #if $coverage_cond.coverage_select == 'yes':
                --coverage ${coverage_cond.coverage_min},${coverage_cond.coverage_max},${coverage_cond.coverage_step}
            #end if
            ${remove_dups}
            #if str( $filter_by_flags.filter_flags ) == "filter":
                #set $filter = $filter_by_flags.require_flags
                
        #set $flags = 0
        #if $filter
            #set $flags = sum(map(int, str($filter).split(',')))
        #end if
    
                --required-flag $flags
                #set $filter = $filter_by_flags.exclude_flags
                
        #set $flags = 0
        #if $filter
            #set $flags = sum(map(int, str($filter).split(',')))
        #end if
    
                --filtering-flag $flags
            #end if
            #if str($gc_depth):
                --GC-depth ${gc_depth}
            #end if
            #if str($insert_size):
                --insert-size ${insert_size}
            #end if
            ## #if $read_group
            ##    -I '$read_group'
            ## #end if
            #if str($read_length):
                --read-length ${read_length}
            #end if
            #if str($most_inserts):
                --most-inserts ${most_inserts}
            #end if
            #if str($trim_quality):
                --trim-quality ${trim_quality}
            #end if

            #if $use_ref:
                --ref-seq "\$reffa"
            #end if
            ## currently not implemented in Galaxy
            ## generates STR_VALUE.bamstat where STR is given by -P and VALUE is a value of the TAG given by -S
            ## needs some discover data sets action...
            ## -P, --split-prefix STR
            ## -S, --split TAG
            $sparse
            
        #if $cond_region.select_region == 'tab':
            -t '$cond_region.targetregions'
        #end if
    
            $remove_overlaps
            #if str($cov_threshold):
                -g $cov_threshold
            #end if
            -@ \$addthreads
            infile
            
        #if $cond_region.select_region == 'text':
            #for $i, $x in enumerate($cond_region.regions_repeat):
               '${x.region}'
            #end for
        #end if
    
            > '$output'

            #if $split_output_cond.split_output_selector == "yes":
                #set outputs_to_split = str($split_output_cond.generate_tables).split(',')
                && mkdir split
                #for s in str($split_output_cond.generate_tables).split(','):
                    && name=`cat '$output' | grep '\^$s' | cut -d'.' -f 1 | sed 's/^# //'`
                    && awk '/\^/{out=0} /\^$s/{out=1} {if(out==1){print $0}}' '$output' | sed 's/Use `grep .*` to extract this part.//' | sed 's/$s\t//' > "split/\$name.tab"
                #end for
            #end if