hicup

#!/usr/bin/perl

use strict;
use warnings;
use Time::HiRes;
use Getopt::Long;
use File::Temp qw/ tempdir /;
use File::Path qw(remove_tree);
use File::Basename;
use FindBin '$Bin';
use lib $Bin;
use Cwd qw(cwd);
use hicup_module;
use hicup_module qw(hashVal check_no_duplicate_filename checkAligner checkAlignerIndices
  quality_checker fileNamer datestampGenerator determineAlignerFormat);

use Data::Dumper;

###################################################################################
###################################################################################
##This file is Copyright (C) 2023, Steven Wingett                                ##
##                                                                               ##
##                                                                               ##
##This file is part of HiCUP.                                                    ##
##                                                                               ##
##HiCUP is free software: you can redistribute it and/or modify                  ##
##it under the terms of the GNU General Public License as published by           ##
##the Free Software Foundation, either version 3 of the License, or              ##
##(at your option) any later version.                                            ##
##                                                                               ##
##HiCUP is distributed in the hope that it will be useful,                       ##
##but WITHOUT ANY WARRANTY; without even the implied warranty of                 ##
##MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the                  ##
##GNU General Public License for more details.                                   ##
##                                                                               ##
##You should have received a copy of the GNU General Public License              ##
##along with HiCUP.  If not, see <http://www.gnu.org/licenses/>.                 ##
###################################################################################
###################################################################################

###################################################################################
###################################################################################
##                                                                               ##
##       This file has been modified to support the following aligners           ##
##                           Dragen, HiSAT2, STAR                                ##
##                                                                               ##
##                           HiCUP+ (HiCUP-Plus)                                 ##
##   Maintained by S. Thomas Kelly (simonthomas.kelly [at] hugp [dot] com)       ##
##                                                                               ##
##         Changes: supports addition configuration options                      ##
##                                                                               ##
###################################################################################
###################################################################################

##########################################################
#Get user-supplied parameters
#Option variables
my %config = (
    ambiguous => '',
    config    => '',
    datestamp => '',
    example   => '',
    help      => '',
    keep      => '',
    nofill    => '',
    outdir    => '',
    quiet     => '',
    version   => '',
    zip       => '',
    bowtie    => '',
    bowtie2   => '',
    digest    => '',
    dragen    => '',
    format    => '',
    hisat2    => '',
    index     => '',
    longest   => '',
    re1       => '',
    shortest  => '',
    samtools  => '',
    star      => '',
    temp      => '',
    threads   => '',
    filenames => '',
    sequence  => '',
    r         => '',
    aligner   => ''    #Not-specified by user - name of aligner to use
);

my $config_result = GetOptions(    #Stores parameters
    "ambiguous"  => \$config{ambiguous},
    "config=s"   => \$config{config},      #Keep, despite being the default
    "datestamp"  => \$config{datestamp},
    "example"    => \$config{example},
    "help"       => \$config{help},
    "keep"       => \$config{keep},
    "nofill"     => \$config{nofill},
    "outdir=s"   => \$config{outdir},
    "version"    => \$config{version},
    "quiet"      => \$config{quiet},
    "zip"        => \$config{zip},
    "bowtie=s"   => \$config{bowtie},
    "bowtie2=s"  => \$config{bowtie2},
    "digest=s"   => \$config{digest},
    "dragen=s"   => \$config{dragen},
    "format=s"   => \$config{format},
    "hisat2=s"   => \$config{hisat2},
    "index=s"    => \$config{index},
    "longest=s"  => \$config{longest},
    "re1=s"      => \$config{re1},
    "shortest=s" => \$config{shortest},
    "samtools=s" => \$config{samtools},
    "star=s"     => \$config{star},
    "temp=s"     => \$config{temp},
    "threads=s"  => \$config{threads},
    "filenames"  => \$config{filenames},
    "r=s"        => \$config{r},
    "sequence"   => \$config{sequence}
);

die "Command line options need to be in the correct format (hicup -help for more details).\n" unless ($config_result);

$config{help} = 1 unless ( hashVal(%config) );    #Print help and exit if no command line parameters

if ( $config{help} ) {
    print while (<DATA>);
    exit(0);
}

#Print version and exit
if ( $config{version} ) {
    print "HiCUP+ v$hicup_module::VERSION\n";
    exit(0);
}

if ( $config{example} ) {
    print_example_config_file('hicup_example.conf');
    exit(0);
}

if( hasval($config{config}) ){
   die "Configuration file '$config{config}' does not exist\n" unless(-e $config{config});
}

print "Starting HiCUP+ pipeline (v$hicup_module::VERSION)\n" unless $config{quiet};
print "PLEASE NOTE: HICUP+ REQUIRES THE R PACKAGES TIDYVERSE AND PLOTLY INSTALLED\n";
print "SEE DOCUMENTATION FOR MORE DETAILS\n";

my @filenames;
if( hasval($config{config}) ){
    @filenames = process_config( $config{config}, \%config );    #Modifies %config and returns an array of the filenames
    if ( scalar @filenames % 2 ) {
        die "There needs to be an even number of files in the configuration file, see hicup --help for more details.\n";
    }
}

if (@ARGV) {
    if ( scalar @ARGV % 2 ) {
        die "There needs to be an even number of files specified in the command line, see hicup --help for more details.\n";
    }
    push( @filenames, @ARGV );                                                               #Add filenames specified in the command line to those in the configuration file
}


#Check files specified and exist
if(@filenames){
   unless ( check_files_exist( \@filenames, 'EXISTS' ) ) {
       die "Please adjust configuration.\n";
   }
}else{
   die "Please specify files to process.\n";
}


my %files = @filenames;                                                                      #%files : hash of paired forward and reverse files
die "Only unique filenames are allowed (irrespective of the path to the file).\n" unless check_no_duplicate_filename(%files);    #Check hash contains no duplicate filenames (irrespective of location)

##########################################################
#Check user-supplied parameters are ok
unless ( check_parameters() ) {
    die "Please change configuration file and/or command-line parameters and/or installation accordingly\n";
}


##########################################################
#Set datestamp
$config{datestamp} = datestampGenerator();   
if( hasval ($config{outdir}) ){
    $config{datestamp} = versioner( $config{datestamp}, $config{outdir}, $config{temp} );
}else{
    $config{datestamp} = versioner( $config{datestamp}, './', $config{temp} );    #Search cwd
}


#########################################################
#Check input files exist and outputfiles do not
my @hicup_final_outfiles = fileNamer( \@filenames, \%config, 'all', 1, 1, 1, 1, 1 );
foreach my $file (@hicup_final_outfiles) {
    $file = $config{outdir} . $file;
}

unless ( check_files_exist( \@hicup_final_outfiles, 'NOT_EXISTS' ) ) {
     die "HiCUP+ will not run until these files have been removed.\n";
}

#########################################################################
#Check files exist
#Check the input files exist
#Sequences, genome, digest
my $files_exist = 1;
unless ( -e $config{digest} ) {
    warn "File '$config{digest}' does not exist\n";
    $files_exist = 0;
}

foreach my $filenameF ( keys %files ) {
    my $filenameR = $files{$filenameF};
    unless ( -e $filenameF ) {
        warn "Sequence file '$filenameF' does not exist\n";
        $files_exist = 0;
    }
    unless ( -e $filenameR ) {
        warn "Sequence file '$filenameR' does not exist\n";
        $files_exist = 0;
    }
}

unless ($files_exist) {
    die "Please adjust the configuration file or add specified files.\n";
}

#Modify the outdir/temp directory setting as required
if ( hasval( $config{temp} ) ) {
    if ( hasval( $config{outdir} ) ) {
        $config{originalOutdir} = $config{outdir};
    } else {
        $config{originalOutdir} = './';
    }

    ( $config{temp} ) = tempdir( 'HiCUP_temp_directory_XXXXXX', DIR => $config{temp} );
    $config{temp}   = $config{temp} . '/';
    $config{outdir} = $config{temp};         #Write future output to the temp directory
}

########################################################################
#Start the pipeline
my $sonication_protocol = 0;                 #Is this a sonication or a double-digest protocol

if ( $config{digest} =~ /\.gz$/ ) {
    open( DIGEST, "gunzip -c $config{digest} |" ) or die "Couldn't read file '$config{digest}' : $!";
} else {
    open( DIGEST, $config{digest} ) or die "Cannot read file read file '$config{digest}' : $!";
}

my $header_line = scalar <DIGEST>;
unless ( $header_line =~ /^.+\t.+\t(.+)\t.+$/ ) {
    die "Genome digest file header is not in the correct format\n";
}

#Obtain the restriction enzyme 1 sequence from the digest file
if ( $1 eq "Restriction_Enzyme2:None" ) {    #Sonication protocol
    $sonication_protocol = 1;
}

my ($re1_string, $re2_string) = ('', '');
if ( $config{sequence} ne '' ) {
    $config{re1} = '';
} else {
    print "Reading genome digest file '$config{digest}' to determine Hi-C restriction enzyme\n" unless $config{quiet};

    unless ( $header_line =~ /^.+\tRestriction_Enzyme1:(.+\[.+\])+\tRestriction_Enzyme2:(.+)\t.+/ ) {
        die "Genome digest file $config{digest} is not in the correct format.\n";
    }
    my ($re1, $re2) = ($1, $2);

    while ($re1 =~ m/((\S+)\s\[(\S+)\])/g) {
        my ($enz_name, $rest_site) = ($2, $3);
        unless ( $rest_site =~ /^[ATCGN\^]+$/ ) {
            die "The restriction site ($rest_site) of $enz_name needs to be a valid DNA sequence\n";
        }
        unless ( ( $rest_site =~ tr/\^// ) == 1 ) {
            die "The restriction site ($rest_site) of $enz_name should contain one cut position, denoted by '^'.\n";
        }
        $re1_string .= $rest_site . "," . $enz_name . ":";
    }
    $re1_string =~ s/\:$//;
    $config{re1} = ($re1_string =~ /:/) ? "Restriction enzymes: " . $re1 : "Restriction enzyme: " . $re1;

    if ($re2 ne "None" ) {
        unless ($re2 =~ /((.+)\s\[(.+?)\])+/) {
            die "The restriction enzyme 2 $re2 is not in the correct format.\n";
        }
        while ($re2 =~ /((\S+)\s\[(\S+)\])/g) {
            my ($enz_name, $rest_site) = ($2, $3);
            unless ( $rest_site =~ /^[ATCGN\^]+$/ ) {
                die "The restriction site ($rest_site) of $enz_name needs to be a valid DNA sequence\n";
            }
            unless ( ( $rest_site =~ tr/\^// ) == 1 ) {
                die "The restriction site ($rest_site) of $enz_name should contain one cut position, denoted by '^'.\n";
            }
            $re2_string .= $rest_site . "," . $enz_name . ":";
        }
        $re2_string =~ s/\:$//;
        $config{re1} .= ($re2_string =~ /:/) ? " - Sonication enzymes: " . $re2 : " - Sonication enzyme: " . $re2;
    }
}
close DIGEST;

#Add command line flags to relevant arguments
$config{ambiguous} = "--ambiguous"                if ( $config{ambiguous} ne '' );
$config{bowtie}    = "--bowtie $config{bowtie}"   if ( $config{bowtie} ne '' );
$config{bowtie2}   = "--bowtie2 $config{bowtie2}" if ( $config{bowtie2} ne '' );
$config{dragen}    = "--dragen $config{dragen}"   if ( $config{dragen} ne '' );
$config{hisat2}    = "--hisat2 $config{hisat2}"   if ( $config{hisat2} ne '' );
$config{star}      = "--star $config{star}"       if ( $config{star} ne '' );
my $datestamp_flag = "--datestamp $config{datestamp}";    #Enables all summary files to have same datestamp
$config{digest} = "--digest $config{digest}" if ( $config{digest} ne '' );
$config{format} = "--format $config{format}" if ( $config{format} ne '' );
if ( $config{longest} ne '' ) {
    $config{longest} = "--longest $config{longest}";
}
my $outdir_flag = '';
$outdir_flag       = "--outdir $config{outdir}"     if ( $config{outdir} );           #This is the flag, not outdir
$config{nofill}    = "--nofill"                     if ( $config{nofill} ne '' );
$re1_string        = "--re1 $re1_string"            if ( $re1_string ne '' );
## $config{samtools}    = "--samtools $config{samtools}"   if ( $config{samtools} ne '' );
$config{seq_trunc} = "--sequence $config{sequence}" if ( $config{sequence} ne '' );
if ( $config{shortest} ne '' ) {
    $config{shortest} = "--shortest $config{shortest}";
}

$config{quiet} = "--quiet" if ( $config{quiet} ne '' );
$config{zip}   = "--zip"   if ( $config{zip} );
$config{index} = "--index $config{index}";
$config{r}     = "-r $config{r}";

#Run each hicup script in turn
#Pipeline will delete intermediate files, unless specified otherwise i.e. flag 'keep:'
#Pipeline will keep all summary files
#Pipeline will zip final and intermediate data output files if 'zip:' selected
#All files will be kept in current folder, but a new folder shall be created for hicup_filter reject sequences

#Truncate sequences
my $filenames_string = filename_string_generator(%files);

#$config{outdir} = $outdir_temp;
my $truncater_version;    #Stores the truncater version number

$truncater_version = `$Bin/hicup_truncater -v` or die "Can't run hicup_truncater to determine version number\n.";

# start timing step
print STDOUT "\nstarting HiCUP truncater\n";
my $start_time1 = [Time::HiRes::gettimeofday()];

if ( $config{zip} ) {
    !system("$Bin/hicup_truncater $config{r} $datestamp_flag $config{nofill} $outdir_flag $config{quiet} $re1_string $config{sequence} $config{threads} $config{zip} $filenames_string")
      or die "Can't run hicup_truncater\n.";
} else {
    !system("$Bin/hicup_truncater $config{r} $datestamp_flag $config{nofill} $outdir_flag $config{quiet} $re1_string $config{sequence} $config{threads} $filenames_string")
      or die "Can't run hicup_truncater\n.";
}

#report timing results
my $diff1 = Time::HiRes::tv_interval($start_time1);
print STDOUT "\ncompleted HiCUP truncater\n";
print STDOUT "\n\ntruncater time: $diff1\n";

#Map and pair sequences
my $ligation = "\@PG\tID:HiCUP+ Truncater\tVN:" . $hicup_module::VERSION;       #Create the hidden ligation flag to send to HiCUP mapper to print the hicup_truncation settings to the output file
if ( $config{re1} ) {
   #We cannot have a ":" in the SAM header at this point, so remove it
   my $colonless_re1 = $config{re1};
   $colonless_re1 =~ s/://g;
    $ligation = $ligation . "\tDS:\"" . "$colonless_re1\"";
} else {
    $ligation = $ligation . "\tDS\"" . "Ligation junction sequences $config{sequence}\"";
}
$ligation = "--ligation \'$ligation\'";

my @truncater_outputfiles = fileNamer( \@filenames, \%config, 'truncater' );
@truncater_outputfiles = arrayAppend( \@truncater_outputfiles, $config{outdir} );
my %truncater_outputfiles_paired = @truncater_outputfiles;
$filenames_string = filename_string_generator(%truncater_outputfiles_paired);

# start timing step
print STDOUT "\nstarting HiCUP mapper\n";
my $start_time2 = [Time::HiRes::gettimeofday()];

print "\nindex file: $config{index} \n" unless $config{quiet};
print "\n$Bin/hicup_mapper $config{r} $config{bowtie} $config{bowtie2} $config{dragen} $config{hisat2} $config{star} $config{ambiguous} $datestamp_flag $config{format} $outdir_flag $config{index} $ligation $config{quiet} $config{threads} $config{zip} $filenames_string \n"  unless $config{quiet};

if ( $config{zip} ) {
    !system(
"$Bin/hicup_mapper $config{r} $config{bowtie} $config{bowtie2} $config{dragen} $config{hisat2} $config{star} $config{ambiguous} $datestamp_flag $config{format} $outdir_flag $config{index} $ligation $config{quiet} $config{threads} $config{zip} $filenames_string"
    ) or die "Can't run hicup_mapper\n.";
} else {
    !system(
"$Bin/hicup_mapper $config{r} $config{bowtie} $config{bowtie2} $config{dragen} $config{hisat2} $config{star} $config{ambiguous} $datestamp_flag $config{format} $outdir_flag $config{index} $ligation $config{quiet} $config{samtools} $config{samtools} $config{threads} $filenames_string"
    ) or die "Can't run hicup_mapper\n.";
}


#report timing results
my $diff2 = Time::HiRes::tv_interval($start_time2);
print STDOUT "\ncompleted HiCUP mapper\n";
print STDOUT "\n\nmapper time: $diff2\n";

#Delete truncated intermediates
unless ( $config{keep} ) {
    foreach my $key ( keys %truncater_outputfiles_paired ) {
        unlink $key or warn "Cannot delete '$key'.\n";                                                                  #First in the pair
        unlink $truncater_outputfiles_paired{$key} or warn "Cannot delete '$truncater_outputfiles_paired{$key}'.\n";    #Second in the pair
    }
}

#Filter files
my @mapper_outputfiles = fileNamer( \@truncater_outputfiles, \%config, 'mapper' );
@mapper_outputfiles = arrayAppend( \@mapper_outputfiles, $config{outdir} );

$filenames_string = filename_string_generator_unpaired(@mapper_outputfiles);                                            #Use this subroutine for unpaired filenames

# start timing step
print STDOUT "\nstarting HiCUP filter\n";
my $start_time3 = [Time::HiRes::gettimeofday()];

if ( $config{zip} ) {
    !system("$Bin/hicup_filter $config{r} $datestamp_flag $config{digest} $config{longest} $outdir_flag $config{quiet} $config{shortest} $config{samtools} $config{threads} $config{zip} $filenames_string")
      or die "Can't run hicup_filter.\n";
} else {
    !system("$Bin/hicup_filter $config{r} $datestamp_flag $config{digest} $config{longest} $outdir_flag $config{quiet} $config{shortest} $config{samtools} $config{threads} $filenames_string")
      or die "Can't run hicup_filter.\n";
}


#report timing results
my $diff3 = Time::HiRes::tv_interval($start_time3);
print STDOUT "\ncompleted HiCUP filter\n";
print STDOUT "\n\nfilter time: $diff3\n";

#Delete paired intermediates
unless ( $config{keep} ) {
    foreach (@mapper_outputfiles) {
        my $file_to_delete = $_;
        if(-e $file_to_delete){   #May not be present if no reads mapped
            unlink $file_to_delete or warn "Cannot delete '$file_to_delete" . "$_'.\n";
        } 
    }
}

#De-duplicate files
my @filter_outputfiles = fileNamer( \@mapper_outputfiles, \%config, 'filter' );    #Subroutine processing hash as an array
@filter_outputfiles = arrayAppend( \@filter_outputfiles, $config{outdir} );   #Do this here as outdir may change below if temp directory specified
if ( hasval( $config{temp} ) ) {    #Use send to outdir, if temp folder specified
    $config{outdir} = $config{originalOutdir};
    if ( hasval( $config{outdir} ) ) {
        $outdir_flag = "--outdir $config{outdir}";
    } else {
        $outdir_flag = '';
    }
}
$filenames_string   = filename_string_generator_unpaired(@filter_outputfiles);     #Use this subroutine for unpaired filenames

# start timing step
print STDOUT "\nstarting HiCUP deduplicator\n";
my $start_time4 = [Time::HiRes::gettimeofday()];

if ( $config{zip} ) {                                                              #$keep not significant since this is the last step in the pipeline
    !system("$Bin/hicup_deduplicator $config{r} $datestamp_flag $outdir_flag $config{quiet} $config{samtools} $config{threads} $config{zip} $filenames_string") or die "Can't run hicup_deduplicator.\n";
} else {
    !system("$Bin/hicup_deduplicator $config{r} $datestamp_flag $outdir_flag $config{quiet} $config{samtools} $config{threads} $filenames_string") or die "Can't run hicup_deduplicator.\n";
}


#report timing results
my $diff4 = Time::HiRes::tv_interval($start_time4);
print STDOUT "\ncompleted HiCUP deduplicator\n";
print STDOUT "\n\ndeduplicator time: $diff4\n";

#Delete filtered files
unless ( $config{keep} ) {
    foreach (@filter_outputfiles) {
        my $file_to_delete = $_;
        if(-e $file_to_delete){    #Sample may contain no valid reads pairs
            unlink $file_to_delete or warn "Cannot delete '$file_to_delete'.\n";
        }
    }
}

#Change deduplication filenames to final hicup output filenames
my @deduplicator_outfiles = fileNamer( \@filter_outputfiles, \%config, 'deduplicator' );
foreach my $dedup_file (@deduplicator_outfiles) {
    if(-e $config{outdir}.$dedup_file){    #File may not be present because no valid di-tags were identified
        my $hicup_final_file = fileNamer( $dedup_file, \%config, 'hicup' );
        my $command = "mv $config{outdir}" . "$dedup_file $config{outdir}" . "$hicup_final_file";
        !system($command) or die "Could not rename file '$command' : $!";
    }
}


#Clean up the --temp output folder and move relevant files to the output directory
if ( hasval( $config{temp} ) ) {
    $config{outdir} = $config{originalOutdir};    #Use original output directory
    !system("mv $config{temp}/*.svg $config{temp}/*summary* $config{outdir}")
      or warn "Could not tidy temporary folder folder with command: 'mv $config{temp}/*.svg $config{temp}/*summary* $config{outdir}' : $!";
    remove_tree $config{temp} or warn "Could not delete temporary folder '$config{temp} with the command: 'remove_tree $config{temp}' : $!";
}

#Create the HiCUP summary reports
#Navigate to the folder where the output is written and run the script there

# start timing step
print STDOUT "\nstarting HiCUP reporter\n";
my $start_time5 = [Time::HiRes::gettimeofday()];

if($sonication_protocol){
    my $current_working_directory = cwd;
    my $hicup_reporter_command =  "$Bin/hicup_reporter .";
    if(hasval($config{outdir})){
        chdir $config{outdir};
        !system($hicup_reporter_command) or die "Can't run hicup_reporter.\n";
        chdir $current_working_directory;
    } else {
        !system($hicup_reporter_command) or die "Can't run hicup_reporter.\n";
    }
} else {
    print "Skipping summary report\n";
    print "HiCUP+ now only collates results for sonication protocol datasets\n";
}

print "HiCUP+ processing complete.\n" unless $config{quiet};

#report timing results
my $diff5 = Time::HiRes::tv_interval($start_time5);
print STDOUT "\ncompleted HiCUP reporter\n";
print STDOUT "\n\nreporter time: $diff5\n";
print STDOUT "\n\ntotal run time: $diff1\n";
exit(0);

#######################################################################################
#Subroutines                                                                          #
#######################################################################################

############################
#Subroutine "check_parameters":
#Check the user supplied parameters are ok
#Uses global variables
sub check_parameters {

    my $parameters_ok = 1;

    $parameters_ok = 0 unless ( checkAligner( \%config ) );
    $parameters_ok = 0 unless ( checkAlignerIndices( \%config ) );

    if ( hasval( $config{format} ) ) {    #Check specified FASTQ format is valid
        $parameters_ok = 0 unless ( determineAlignerFormat( $config{format} ) );
    } else {                              #If FASTQ quality format not specified, determine this automatically
        my @sequence_files = %files;
        warn "FASTQ quality format not specified, analysing file '$sequence_files[0]' to predict file format used\n";
        $config{format} = quality_checker( $sequence_files[0] );

        if ( $config{format} eq 0 ) {
            die "Unable to determine FASTQ quality format, please specify this and re-run HiCUP.\n";
        } else {
            warn "FASTQ quality set to '$config{format}'\n";
        }
    }

    unless ( $config{digest} ) {
        print "Please specify a genome digest file.\n";
        $parameters_ok = 0;
    }
    unless ( hasval( $config{index} ) ) {
        print "Please specify reference genome indices for mapping.\n";
        $parameters_ok = 0;
    }

    if ( $config{longest} ne '' ) {
        unless ( $config{longest} =~ /^\d+$/ ) {
            print "'Longest' needs to specify the insert length in base pairs (i.e. is an integer)\n";
            $parameters_ok = 0;
        }
    }

    if ( $config{shortest} ne '' ) {
        unless ( $config{shortest} =~ /^\d+$/ ) {
            print "'Shortest' needs to specify the insert length in base pairs (i.e. is an integer)\n";
            $parameters_ok = 0;
        }
    }

    if ( ( $config{longest} ne '' ) and ( $config{shortest} ne '' ) ) {
        if ( $config{longest} =~ /^\d+$/ and $config{shortest} =~ /^\d+$/ ) {
            unless ( $config{shortest} < $config{longest} ) {
                print "The shortest allowable insert size should be smaller than the longest allowable insert size.\n";
                $parameters_ok = 0;
            }
        }
    }

    if ( $config{sequence} ne '' ) {
        my @ligation_sequences = split( /,/, $config{sequence} );    #Add directly to the global variable

        #Check all the ligation sequences are in the valid format ATCG_ATCG, where before '_' denotes sequence present in the genome reference
        foreach my $sequence (@ligation_sequences) {
            unless ( $sequence =~ /^[ATCG_]+$/ ) {
                warn "The ligation sequence $sequence is not a valid DNA sequence.\n";
                $parameters_ok = 0;
            }
            unless ( ( $sequence =~ tr/\_// ) == 1 ) {
                warn "The ligation sequence should define the position where the ligation sequence no longer matches the genome reference - denoted by \'_\'.\n";
                $parameters_ok = 0;
            }
        }
    }

    if ( $config{threads} ne '' ) {
        unless ( $config{threads} =~ /^\d+$/ ) {
            print "The specified 'Threads' value needs to be an integer.\n";
            $parameters_ok = 0;
        } else {
            $config{threads} = "-threads $config{threads}";
        }
    } else {
        $config{threads} = '-threads 1';
    }

    if ( $config{quiet} ne '' ) {    #May be defined as 0 in the config file
        unless ( $config{quiet} ) {
            $config{quiet} = '';
        }
    }

    if ( $config{nofill} ne '' ) {    #May be defined as 0 in the config file
        unless ( $config{nofill} ) {
            $config{nofill} = '';
        }
    }

    #Check SAMtools is installed
    if ( !system "which samtools >/dev/null 2>&1" ) {
        $config{samtools} = `which samtools`;
        chomp $config{samtools};
    } else {
        warn "Could not find SAMtools (http://samtools.sourceforge.net), please install if you wish to compress SAM files to BAM format\n";

    }

    #Check R installed
    checkR( \%config );

    #Check the output directory exists and allows write access
    if ( $config{outdir} ne '' ) {
        unless ( -d $config{outdir} ) {
            warn "Output directory '$config{outdir}' does not exist\n";
            $parameters_ok = 0;
        }

        unless ( $config{outdir} =~ /\/$/ ) {    #Make sure that $config{outdir} ends with the forward slash character
            $config{outdir} .= '/';
        }
    }

    if ( hasval( $config{temp} ) ) {
        if ( $config{keep} ) {
            warn "Option --temp may not be specified if --keep is also specified\n";
            $parameters_ok = 0;
        } else {
            unless ( -d $config{temp} ) {
                warn "Temporary directory '$config{temp}' does not exist\n";
                $parameters_ok = 0;
            }
        }
    }

    return $parameters_ok;
}

########################################
#Subroutine "filename_string_generator":
#Receives a hash of paired filenames which it uses to creates a string of the files to
#process. It also checks if those files are empty and if so removes them from the
#pipeline.
sub filename_string_generator {

    my (%filename_pairs) = @_;
    my $filenames_string = "";

    foreach my $filename1 ( keys %filename_pairs ) {
        my $filename2 = $filename_pairs{$filename1};

        #$filename1 = $config{outdir} . $filename1;
        #$filename2 = $config{outdir} . $filename2;

        my $filename1_contains_data = 0;
        my $filename2_contains_data = 0;

        unless ( ( -e $filename1 ) and ( -e $filename2 ) ) {    #Only check filesize of files that exist
            next;
        }

        #Identify empty unzipped/zipped files
        if ( $filename1 =~ /\.gz$/ ) {
            if ( ( -s $filename1 ) > 100 ) {                    #Only perform check on small files
                $filename1_contains_data = 1;
            } else {
                open( FILENAME1, "gunzip -c $filename1 |" ) or die "Couldn't read file '$filename1' : $!";
                while (<FILENAME1>) {

                    if ( $_ =~ /\S/ ) {
                        $filename1_contains_data = 1;
                    }
                }
                close FILENAME1;
            }

            if ( ( -s $filename2 ) > 100 ) {    #    Only perform check on small files
                $filename2_contains_data = 1;
            } else {
                open( FILENAME2, "gunzip -c $filename2 |" ) or die "Couldn't read file '$filename2' : $!";
                while (<FILENAME2>) {
                    if ( $_ =~ /\S/ ) {
                        $filename2_contains_data = 1;
                    }
                }
                close FILENAME2;
            }
        } else {
            if ( -s $filename1 ) {
                $filename1_contains_data = 1;
            }
            if ( -s $filename2 ) {
                $filename2_contains_data = 1;
            }
        }
        unless ($filename1_contains_data) {
            warn "$filename1 contains no data\n";
        }
        unless ($filename2_contains_data) {
            warn "$filename2 contains no data\n";
        }
        unless ( $filename1_contains_data + $filename2_contains_data == 2 ) {
            warn "$filename1 and $filename2 shall be processed no further\n";
        } else {
            $filenames_string .= "$filename1 $filename2 ";
        }
    }

    if ( $filenames_string eq "" ) {
        die "The remaining files in the HiCUP pipeline contain no data. Processing terminated.\n";
    } else {
        return $filenames_string;
    }
}

#################################################
#Subroutine "filename_string_generator_unpaired":
#receives filenames that do not need pairing (e.g. '.pair' files) which it uses to create
#a string of the files to process. It also checks if those files are empty and, if so,
#removes them from the pipeline.
sub filename_string_generator_unpaired {

    my @input_files      = @_;
    my $filenames_string = "";

    foreach my $inputfile (@input_files) {
        $inputfile = $inputfile;

        unless ( -e $inputfile ) {    #Only check filesize of files that exist
            next;
        }

        if ( ( -s $inputfile ) < 10000 ) {    #Only check whether small files contain data

            #Open the inputfile
            if ( $inputfile =~ /\.gz$/ ) {    #Compressed SAM file
                open( INPUTFILE, "gunzip -c $inputfile |" ) or die "Couldn't read file '$inputfile' : $!";
            } elsif ( $inputfile =~ /\.bam$/ ) {    #BAM file
                open( INPUTFILE, "$config{samtools} view -h $inputfile |" ) or die "Couldn't read file '$inputfile' : $!";
            } else {                                #SAM file
                open( INPUTFILE, $inputfile ) or die "Couldn't read file '$inputfile' : $!";
            }

            #Check whether the file contains data i.e. a line not beginning with '@'
            my $data_present = 0;
            while (<INPUTFILE>) {
                if (/^@/) {
                    next;
                }
                if (/\S/) {
                    $data_present = 1;              #Data found
                    last;
                }
            }

            if ($data_present) {
                $filenames_string .= "$inputfile ";
            } else {
                warn "$inputfile contains no data\n";
            }
        } else {
            $filenames_string .= "$inputfile ";     #File > 10000 in size, so assumed to contain data
        }
    }

    if ( $filenames_string eq "" ) {
        die "All the files in the HiCUP pipeline have been removed for containing no data.\n";
    } else {
        return $filenames_string;
    }
}

__DATA__

HiCUP homepage: www.bioinformatics.babraham.ac.uk/projects/hicup

SYNOPSIS

hicup is used to run the whole HiCUP pipeline

hicup [OPTIONS]... --config [CONFIGURATION FILE]...
hicup [OPTIONS]... [FILES]

FUNCTION

The HiCUP pipeline receives FASTQ files and generates Hi-C paired read (di-tag)
files. The hicup script regulates the pipeline, passing output from one script to 
the next. 

The HiCUP pipeline uses a configuration file to set the mapping and filtering 
parameters. To produce an example configuration file:

hicup --example

Parameters may also be passed to HiCUP via the command line, overriding those 
specified in the configuration file. FASTQ file pairs should be place next to 
each other when using the command line, or on adjacent lines in the configuration 
file.

COMMAND LINE OPTIONS

--bowtie       Specify the path to Bowtie
--bowtie2      Specify the path to Bowtie 2
--dragen       Specify the path to DRAGEN
--config       Specify the configuration file
--digest       Specify the digest file listing restriction fragment co-ordinates
--example      Produce an example configuration file
--format       Specify FASTQ format
               Options: Sanger, Solexa_Illumina_1.0, Illumina_1.3, Illumina_1.5
--help         Print help message and exit
--hisat2       Specify the path to HISAT 2
--index        Path to the relevant reference genome Bowtie/Bowtie2 indices
--keep         Keep intermediate pipeline files
--longest      Maximum allowable insert size (bps)
--nofill       Hi-C protocol did NOT include a fill-in of sticky ends prior to 
               ligation step and therefore FASTQ reads shall be truncated at the 
               Hi-C restriction enzyme cut site (if present) sequence is 
               encountered
--outdir       Directory to write output files
--quiet        Suppress progress reports (except warnings)
--samtools     Specify the path to samtools
--star         Specify the path to STAR
--shortest     Minimum allowable insert size (bps)
--temp         Write intermediate files (i.e. all except summary files and files
               generated by HiCUP Deduplicator) to a specified directory
--threads      Specify the number of threads, allowing simultaneous processing 
               of multiple files
--version      Print the program version and exit
--zip          Compress output

Full instructions on running the pipeline can be found at:
www.bioinformatics.babraham.ac.uk/projects/hicup

Steven Wingett, Babraham Institute, Cambridge, UK