--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/ITSx	Tue Mar 24 12:02:48 2015 -0400
@@ -0,0 +1,2695 @@
+#!/usr/bin/perl
+# ITSx ITS Extractor
+$app_title     = "ITSx -- Identifies ITS sequences and extracts the ITS region";
+$app_author    = "Johan Bengtsson-Palme et al., University of Gothenburg";
+$app_version   = "1.0.11";
+$app_message   = "";
+# ----------------------------------------------------------------- #
+
+# License information
+$license =
+"    ITSx - ITS Extractor -- Identifies ITS sequences and extracts the ITS region\
+    Copyright (C) 2012-2014 Johan Bengtsson-Palme et al.\
+\
+    This program 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.\
+\
+    This program 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 this program, in a file called 'license.txt'\
+    If not, see: http://www.gnu.org/licenses/.\
+";
+
+## BUGS:
+$bugs = "New features in this version ($app_version):\
+- None\
+\
+Fixed bugs in this version ($app_version):\
+- Fixed a bug causing newline characters to be occasionally skipped in the 'its1.full_and_partial.fasta' FASTA output file when the '--anchor' option was used\
+\
+Known bugs in this version ($app_version):\
+- None\
+";
+
+## OPTIONS:
+$options = "\
+-i {file} : DNA FASTA input file to investigate\
+-o {file} : Base for the names of output file(s)\
+-p {directory} : A path to a directory of HMM-profile collections representing ITS conserved regions, default is in the same directory as ITSx itself\
+--date {T or F} : Adds a date and time stamp to the output directory, off (F) by default\
+--reset {T or F} : Re-creates the HMM-database before ITSx is run, off (F) by default\
+
+Sequence selection options:\
+-t {character code} : Profile set to use for the search, see the User's Guide (comma-separated), default is all\
+-E {value} : Domain E-value cutoff for a sequence to be included in the output, default = 1e-5\
+-S {value} : Domain score cutoff for a sequence to be included in the output, default = 0\
+-N {value} : The minimal number of domains that must match a sequence before it is included, default = 2\
+--selection_priority {sum, domains, eval, score} : Selects what will be of highest priority when determining the origin of the sequence, default is sum\
+--search_eval {value} : The E-value cutoff used in the HMMER search, high numbers may slow down the process, cannot be used with the --search_score option, default is 0.01\
+--search_score {value} : The score cutoff used in the HMMER search, low numbers may slow down the process, cannot be used with the --search_eval option, default is to used E-value cutoff, not score\
+--allow_single_domain {e-value,score or F} : Allow inclusion of sequences that only find a single domain, given that they meet the given E-value and score thresholds, on with parameters 1e-9,0 by default\
+--allow_reorder {T or F} : Allows profiles to be in the wrong order on extracted sequences, off (F) by default\
+--complement {T or F} : Checks both DNA strands against the database, creating reverse complements, on (T) by default\
+--cpu {value} : the number of CPU threads to use, default is 1\
+--multi_thread {T or F} : Multi-thread the HMMER-search, on (T) if number of CPUs (--cpu option > 1), else off (F) by default\
+--heuristics {T or F} : Selects whether to use HMMER's heuristic filtering, off (F) by default\
+
+Output options:\
+--summary {T or F} : Summary of results output, on (T) by default\
+--graphical {T or F} : 'Graphical' output, on (T) by default\
+--fasta {T or F} : FASTA-format output of extracted ITS sequences, on (T) by default\
+--preserve {T or F} : Preserve sequence headers in input file instead of printing out ITSx headers, off (F) by default\
+--save_regions {SSU,ITS1,5.8S,ITS2,LSU,all,none} : A comma separated list of regions to output separate FASTA files for, 'ITS1,ITS2' by default\
+--anchor {integer or HMM} : Saves an additional number of bases before and after each extracted region. If set to 'HMM' all bases matching the corresponding HMM will be output, default = 0\
+--only_full {T or F} : If true, output is limited to full-length regions, off (F) by default\
+--partial {integer} : Saves additional FASTA-files for full and partial ITS sequences longer than the specified cutoff, default = 0 (off)\
+--concat {T or F} : Saves a FASTA-file with concatenated ITS sequences (with 5.8S removed), off (F) by default\
+--minlen {integer} : Minimum length the ITS regions must be to be outputted in the concatenated file (see above), default = 0\
+--positions {T or F} : Table format output containing the positions ITS sequences were found in, on (T) by default\
+--table {T or F} : Table format output of sequences containing probable ITS sequences, off (F) by default\
+--not_found {T or F} : Saves a list of non-found entries, on (T) by default\
+--detailed_results {T or F} : Saves a tab-separated list of all results, off (F) by default\
+--truncate {T or F} : Truncates the FASTA output to only contain the actual ITS sequences found, on (T) by default\
+--silent {T or F} : Supresses printing progress info to stderr, off (F) by default\
+--graph_scale {value} : Sets the scale of the graph output, if value is zero, a percentage view is shown, default = 0\
+--save_raw {T or F} : Saves all raw data for searches etc. instead of removing it on finish, off (F) by default\
+
+-h : displays this help message\
+--help : displays this help message\
+--bugs : displays the bug fixes and known bugs in this version of ITSx\
+--license : displays licensing information\
+";
+
+
+## Print title message
+print STDERR "$app_title\nby $app_author\nVersion: $app_version\n$app_message";
+print STDERR "-----------------------------------------------------------------\n";
+
+## Setup default variable values
+use List::Util qw(first max maxstr min minstr reduce shuffle sum);
+
+$bindir = $0;
+$bindir =~ s/_x//;
+$input = "";
+$output = "ITSx_out";
+$hmmscan = "";
+$profileDB = "$bindir\_db/HMMs";
+$type = "all";
+$E = 1e-5;
+$S = 0;
+$N = 2;
+$priority = "sum";
+$search_eval = 0.01;
+$search_score = "";
+$allow_single_E = 1e-9;
+$allow_single_score = 0;
+#$allow_single_E = -1;   # Turns off single-domain matching by E-value
+#$allow_single_score = 0;   # Turns off single-domain matching by score
+$allow_reorder = 0;
+$complement = 1;
+$cpu = 1;
+$multi_thread = "unset";
+$heuristics = 0;
+$out_sum = 1;
+$out_graph = 1;
+$out_fasta = 1;
+$out_preserve = 0;
+$out_ssu = 0;
+$out_its1 = 1;
+$out_its2 = 1;
+$out_58S = 0;
+$out_lsu = 0;
+$out_pos = 1;
+$out_table = 0;
+$out_not = 1;
+$out_date = 0;
+$out_joined = 0;
+$out_results = 0;
+$out_partial = 0;
+$out_concat = 0;
+$concat_minlen = 0;
+$truncate = 1;
+$anchor = 0;
+$only_full = 0;
+$graph_scale = 0;
+$debug = 0;
+$reset = 0;
+
+## Read command-line options
+for ($i = 0; $i <= scalar(@ARGV); $i++) {   # Goes through the list of arguments
+  $arg = @ARGV[$i];   # Stores the current argument in $arg
+
+  if ($arg eq "-i") {   # Read input files from -i flag
+    $i++;
+    $input = @ARGV[$i];
+  }
+  if ($arg eq "-o") {   # Read output files from -o flag
+    $i++;
+    $output = @ARGV[$i];
+  }
+  if ($arg eq "-p") {   # Read profile database from -p flag
+    $i++;
+    $profileDB = @ARGV[$i];
+  }
+  if ($arg eq "--hmmscan") {   # Read pre-computed hmmscan output file from --hmmscan flag ('undocumented' feature)
+    $i++;
+    $hmmscan = @ARGV[$i];
+  }
+  if ($arg eq "--date") {   # Determine whether or not to add a date stamp based on the --date flag
+    $i++;
+    if (substr(@ARGV[$i],0,1) =~ m/^[Ff0]/) {   # Check if argument begins with "F", "f", or "0"
+      $out_date = 0;
+    } else {
+      $out_date = 1;
+    }
+  }
+  if ($arg eq "--reset") {   # Reset HMM database?
+    $i++;
+    if (substr(@ARGV[$i],0,1) =~ m/^[Ff0]/) {   # Check if argument begins with "F", "f", or "0"
+      $reset = 0;
+    } else {
+      $reset = 1;
+    }
+  }
+
+  if ($arg eq "-t") {   # Select what types of ITSs to look for using the -t flag
+    $i++;
+    $type = @ARGV[$i];
+  }
+  if ($arg eq "-E") {   # Set the E-value cutoff using the -E flag
+    $i++;
+    $E = @ARGV[$i];
+  }
+  if ($arg eq "-S") {   # Set the score cutoff using the -S flag
+    $i++;
+    $S = @ARGV[$i];
+  }
+  if ($arg eq "-N") {   # Set the number of found domains cutoff using the -N flag
+    $i++;
+    $N = @ARGV[$i];
+  }
+  if ($arg eq "--selection_priority") {   # Set how to order the ITS types using the --selection_priority flag
+    $i++;
+    $priority = @ARGV[$i];
+  }
+  if ($arg eq "--search_eval") {   # Set the E-value cutoff for the HMMER search using the --search_eval flag
+    $i++;
+    $search_eval = @ARGV[$i];
+    $search_score = "";   #  Turns off score cutoff
+  }
+  if ($arg eq "--search_score") {   # Set the score cutoff for the HMMER search using the --search_score flag
+    $i++;
+    $search_score = @ARGV[$i];
+    $search_eval = "";   #  Turns off E-value cutoff
+  }
+  if ($arg eq "--allow_single_domain") {   # Determine whether or not to allow single domain matches based on the --allow_single_domain flag
+    $i++;
+    if (substr(@ARGV[$i],0,1) =~ m/^[Ff0]/) {   # Check if argument begins with "F", "f", or "0"
+      $allow_single_E = -1;   # Turns off single-domain matching by E-value
+      $allow_single_score = 0;   # Turns off single-domain matching by score
+    } else {
+      ($allow_single_E,$allow_single_score) = split(',',@ARGV[$i]);   # Turns on single-domain matching, assigning the first given value as the E-value cutoff, and the second as score cutoff
+    }
+  }
+  if ($arg eq "--allow_reorder") {   # Determine whether or not to allow the domains to be in the wrong order based on the --allow_reorder flag
+    $i++;
+    if (substr(@ARGV[$i],0,1) =~ m/^[Tt1]/) {   # Check if argument begins with "T", "t", or "1"
+      $allow_reorder = 1;
+    } else {
+      $allow_reorder = 0;
+    }
+  }
+  if ($arg eq "--complement") {   # Determine whether or not to scan the complementary strand of the input file based on the --complement flag
+    $i++;
+    if (substr(@ARGV[$i],0,1) =~ m/^[Tt1]/) {   # Check if argument begins with "T", "t", or "1"
+      $complement = 1;
+    } else {
+      $complement = 0;
+    }
+  }
+  if ($arg eq "--cpu") {   # Set the number of CPUs to use based on the --cpu flag
+    $i++;
+    $cpu = @ARGV[$i];
+  }
+  if ($arg eq "--multi_thread") {   # Determine whether or not to multi-thread the HMMER step based on the --multi_thread flag
+    $i++;
+    if (substr(@ARGV[$i],0,1) =~ m/^[Tt1]/) {   # Check if argument begins with "T", "t", or "1"
+      $multi_thread = 1;
+    } else {
+      $multi_thread = 0;
+    }
+  }
+  if ($arg eq "--heuristics") {   # Determine whether or not to use HMMER's heuristic filtering based on the --heuristics flag
+    $i++;
+    if (substr(@ARGV[$i],0,1) =~ m/^[Tt1]/) {   # Check if argument begins with "T", "t", or "1"
+      $heuristics = 1;
+    } else {
+      $heuristics = 0;
+    }
+  }
+
+  if ($arg eq "--summary") {   # Determine whether or not to output a summary based on the --summary flag
+    $i++;
+    if (substr(@ARGV[$i],0,1) =~ m/^[Tt1]/) {   # Check if argument begins with "T", "t", or "1"
+      $out_sum = 1;
+    } else {
+      $out_sum = 0;
+    }
+  }
+  if ($arg eq "--graphical") {   # Determine whether or not to output a graphical representation of matches based on the --graphical flag
+    $i++;
+    if (substr(@ARGV[$i],0,1) =~ m/^[Tt1]/) {   # Check if argument begins with "T", "t", or "1"
+      $out_graph = 1;
+    } else {
+      $out_graph = 0;
+    }
+  }
+  if ($arg eq "--detailed_results") {   # Determine whether or not to output a detailed results list, based on the --detailed_results flag
+    $i++;
+    if (substr(@ARGV[$i],0,1) =~ m/^[Tt1]/) {   # Check if argument begins with "T", "t", or "1"
+      $out_results = 1;
+    } else {
+      $out_results = 0;
+    }
+  }
+  if ($arg eq "--partial") {   # Set the full-and-partial cutoff
+    $i++;
+    $out_partial = @ARGV[$i];
+  }
+  if ($arg eq "--anchor") {   # Set the length of the sequence "anchors"
+    $i++;
+    $anchor = @ARGV[$i];
+  }
+  if ($arg eq "--only_full") {   # Output only full-length regions
+    $i++;
+    if (substr(@ARGV[$i],0,1) =~ m/^[Tt1]/) {   # Check if argument begins with "T", "t", or "1"
+      $only_full = 1;
+    } else {
+      $only_full = 0;
+    }
+  }
+  if ($arg eq "--save_regions") {   # Determine which regions to output FASTA files for based on the --save_regions flag
+    $i++;
+    @save_regions = split(',',uc(@ARGV[$i]));
+    $out_ssu = 0;
+    $out_its1 = 0;
+    $out_its2 = 0;
+    $out_58S = 0;
+    $out_lsu = 0;
+    foreach $save_region (@save_regions) {
+      if ($save_region eq "SSU") {
+	$out_ssu = 1;
+      }
+      if ($save_region eq "ITS1") {
+	$out_its1 = 1;
+      }
+      if ($save_region eq "5.8S") {
+	$out_58S = 1;
+      }
+      if ($save_region eq "ITS2") {
+	$out_its2 = 1;
+      }
+      if ($save_region eq "LSU") {
+	$out_lsu = 1;
+      }
+      if ($save_region eq "ALL") {
+	$out_ssu = 1;
+	$out_its1 = 1;
+	$out_its2 = 1;
+	$out_58S = 1;
+	$out_lsu = 1;
+      }
+      if ($save_region eq "NONE") {
+	$out_ssu = 0;
+	$out_its1 = 0;
+	$out_its2 = 0;
+	$out_58S = 0;
+	$out_lsu = 0;
+      }
+    }
+  }
+  if ($arg eq "--positions") {   # Determine whether or not to output a positions file based on the --positions flag
+    $i++;
+    if (substr(@ARGV[$i],0,1) =~ m/^[Tt1]/) {   # Check if argument begins with "T", "t", or "1"
+      $out_pos = 1;
+    } else {
+      $out_pos = 0;
+    }
+  }
+  if ($arg eq "--concat") {   # Determine whether or not to output a concatednated ITS1 + ITS2 file
+    $i++;
+    if (substr(@ARGV[$i],0,1) =~ m/^[Tt1]/) {   # Check if argument begins with "T", "t", or "1"
+      $out_concat = 1;
+    } else {
+      $out_concat = 0;
+    }
+  }
+  if ($arg eq "--minlen") {   # Set the min length of the combined ITS1 and ITS2 sequences for concatenation
+    $i++;
+    $concat_minlen = @ARGV[$i];
+  }
+  if ($arg eq "--fasta") {   # Determine whether or not to output FASTA-files based on the --fasta flag
+    $i++;
+    if (substr(@ARGV[$i],0,1) =~ m/^[Tt1]/) {   # Check if argument begins with "T", "t", or "1"
+      $out_fasta = 1;
+    } else {
+      $out_fasta = 0;
+    }
+  }
+  if ($arg eq "--preserve") {   # Determine whether or not to preserve FASTA-headers based on the --preserve flag
+    $i++;
+    if (substr(@ARGV[$i],0,1) =~ m/^[Tt1]/) {   # Check if argument begins with "T", "t", or "1"
+      $out_preserve = 1;
+    } else {
+      $out_preserve = 0;
+    }
+  }
+  if ($arg eq "--joined") {   # Determine whether or not to output a FASTA-file containing ALL sorts of output sequences (for debugging)
+    $i++;
+    if (substr(@ARGV[$i],0,1) =~ m/^[Tt1]/) {   # Check if argument begins with "T", "t", or "1"
+      $out_joined = 1;
+    } else {
+      $out_joined = 0;
+    }
+  }
+  if ($arg eq "--table") {   # Determine whether or not to output tables of all potential matches based on the --table flag
+    $i++;
+    if (substr(@ARGV[$i],0,1) =~ m/^[Tt1]/) {   # Check if argument begins with "T", "t", or "1"
+      $out_table = 1;
+    } else {
+      $out_table = 0;
+    }
+  }
+  if ($arg eq "--not_found") {   # Determine whether or not to output a list of sequences that are not ITSs based on the --not_found flag
+    $i++;
+    if (substr(@ARGV[$i],0,1) =~ m/^[Tt1]/) {   # Check if argument begins with "T", "t", or "1"
+      $out_not = 1;
+    } else {
+      $out_not = 0;
+    }
+  }
+  if ($arg eq "--silent") {   # Determine whether or not to output anything to the screen based on the --silent flag
+    $i++;
+    if (substr(@ARGV[$i],0,1) =~ m/^[Tt1]/) {   # Check if argument begins with "T", "t", or "1"
+      $silent = 1;
+    } else {
+      $silent = 0;
+    }
+  }
+  if ($arg eq "--graph_scale") {   # Set the scale of the graphical output based on the --graph_scale flag
+    $i++;
+    $graph_scale = @ARGV[$i];
+  }
+  if ($arg eq "--save_raw") {   # Determine whether or not to save all the raw intermediate data based on the --save_raw flag
+    $i++;
+    if (substr(@ARGV[$i],0,1) =~ m/^[Tt1]/) {   # Check if argument begins with "T", "t", or "1"
+      $save_raw = 1;
+    } else {
+      $save_raw = 0;
+    }
+  }
+
+  ## If "-h" or "--help" are among the options, output usage data and options
+  if (($arg eq "-h") || ($arg eq "--help")) {
+    print "Usage: ITSx -i <input file> -o <output file>\nOptions:$options";
+    print "-----------------------------------------------------------------\n";
+    exit;   # Exit ITSx
+  }
+
+  ## If "--bugs" is among the options, output bugs and features information
+  if ($arg eq "--bugs") {
+    print "$bugs\n";
+    exit;   # Exit ITSx
+  }
+
+  ## If "--license" is among the options, output license information
+  if ($arg eq "--license") {
+    print "$license\n";
+    exit;   # Exit ITSx
+  }
+
+  if ($arg eq "--debug") {   # Run ITSx in debug mode
+    $debug = 1;
+  }
+  if ($arg eq "--pipeline") {   # Run ITSx in pipeline mode
+    $pipeline = 1;
+  }
+}
+
+## Setup some variables dependent on input
+
+if ($multi_thread eq "unset") {   # If the multi-thread option is not set
+  if ($cpu > 1) {   # Then if the number of CPUs used > 1, then multi-thread HMMER searches
+    $multi_thread = 1;
+  } else {   # Else, run HMMER searches sequentially on one CPU
+    $multi_thread = 0;
+  }
+}
+
+if ($hmmscan ne "") {   # If a pre-computed hmmscan output is supplied
+  $output = $hmmscan;   # Then set the base of the output directory name to be the same as that hmmscan output file
+}
+
+## Check for binaries
+
+chomp($path = `which hmmpress`);   # Get the path for hmmpress
+if ($path eq "") {   # If the path is empty, then show an error message and exit ITSx
+  print STDERR "FATAL ERROR :: Could not locate HMMER binaries! It seems that hmmpress is not installed properly.\
+Consult the manual for installation instructions. Note that HMMER3 is required. Previous HMMER-versions will not work.\
+This error is fatal, and ITSx will now abort.\n";
+  print STDERR "-----------------------------------------------------------------\n";
+  exit;
+}
+
+chomp($path = `which hmmscan`);   # Get the path for hmmscan
+if ($path eq "") {   # If the path is empty, then show an error message and exit ITSx
+  print STDERR "FATAL ERROR :: Could not locate HMMER binaries! It seems that hmmscan is not installed properly.\
+Consult the manual for installation instructions. Note that HMMER3 is required. Previous HMMER-versions will not work.\
+This error is fatal, and ITSx will now abort.\n";
+  print STDERR "-----------------------------------------------------------------\n";
+  exit;
+}
+
+
+## Check for database
+chomp($errormsg = `ls $profileDB* 2>&1 1>/dev/null`);   # Get the error msg when looking for the profile database
+if (substr($errormsg,0,4) eq "ls: ") {   # If the error message begins with "ls: ", then show an error message and exit ITSx
+  print STDERR "FATAL ERROR :: The specified profile database could not be found.\
+Consult the manual for installation instructions.\
+This error is fatal, and ITSx will now abort.\n";
+  print STDERR "-----------------------------------------------------------------\n";
+  exit;
+}
+
+if ($pipeline == 0) {   # If ITSx is not run in pipeline mode (i.e. from ITSx)
+  if ($out_date == 1) {   # If a date and time stamp should be supplied
+    ($sec,$min,$hour,$mday,$mon,$year,$wday,$yday,$isdst) = localtime(time);   # Get the date and time
+    $year = $year + 1900;   # Format the year
+    $mon = $mon + 1;   # Format the month
+    if ($mon < 10) {   # Add a zero to the month, if needed
+      $mon = "0" . $mon;
+    }
+    if ($mday < 10) {   # Add a zero to the day, if needed
+      $mday = "0" . $mday;
+    }
+    if ($hour < 10) {   # Add a zero to the hour, if needed
+      $hour = "0" . $hour;
+    }
+    if ($min < 10) {   # Add a zero to the minute, if needed
+      $min = "0" . $min;
+    }
+    $outputDate = ".$year\-$mon\-$mday\_$hour.$min";   # Create a date and time stamp
+    $outputDate =~ s./.-.g;   # Remove any potential slashes in the output name (as this will confuse ITSx's file naming)
+    $output = $output . $outputDate;   # Add the date and time stamp top the output base name
+  }
+}
+
+$tempDir = "ITSx_temp_directory__$output";   # Setup a temporary directory variable
+$tempDir =~ s./.-.g;   # Remove any potential slashes in the output name (as this will confuse ITSx's file naming)
+
+if ($pipeline == 0) {   # If not running in pipeline mode
+  ## Create a summary file
+  if ($out_sum == 1) {   # If summary output is on
+    $now = localtime;   # Get the current time
+    open (SUMMARY, ">$output.summary.txt");   # Create the summary file
+    print SUMMARY "ITSx run started at $now.\n";   # Output the starting time for the analysis
+    print SUMMARY "-----------------------------------------------------------------\n";
+    close (SUMMARY);   # Close summary file
+  }
+}
+
+## Create a temporary directory for ITSx
+if ($pipeline == 0) {   # If ITSx is not run in pipeline mode (i.e. from ITSx)
+  `mkdir $tempDir 2> /dev/null`;   # Create a temporary directory
+}
+
+
+## Prepare profile database
+## Get the current time and output info message
+$now = localtime;
+if ($silent == 0) {
+  print STDERR "$now : Preparing HMM database (should be quick)...\n";
+}
+
+## Setup profile index
+%profileIndex = {};
+
+# A = Alveolata
+# B = Bryophyta
+# C = Bacillariophyta
+# D = Amoebozoa
+# E = Euglenozoa
+# F = Fungi
+# G = Chlorophyta (green algae)
+# H = Rhodophyta (red algae)
+# I = Phaeophyceae (brown algae)
+# L = Marchantiophyta (liverworts)
+# M = Metazoa
+# N = Microsporidia
+# O = Oomycota
+# P = Haptophyceae (prymnesiophytes)
+# Q = Raphidophyceae
+# R = Rhizaria
+# S = Synurophyceae
+# T = Tracheophyta (higher plants)
+# U = Eustigmatophyceae
+# X = Apusozoa
+# Y = Parabasalia
+
+$profileIndex{"A"} = "alveolates";
+$profileIndex{"B"} = "bryophyta";
+$profileIndex{"C"} = "bacillariophyta";
+$profileIndex{"D"} = "amoebozoa";
+$profileIndex{"E"} = "euglenozoa";
+$profileIndex{"F"} = "fungi";
+$profileIndex{"G"} = "chlorophyta";
+$profileIndex{"H"} = "rhodophyta";
+$profileIndex{"I"} = "phaeophyceae";
+$profileIndex{"J"} = "undefined";
+$profileIndex{"K"} = "undefined";
+$profileIndex{"L"} = "marchantiophyta";
+$profileIndex{"M"} = "metazoa";
+$profileIndex{"N"} = "microsporidia";
+$profileIndex{"O"} = "oomycota";
+$profileIndex{"P"} = "haptophyceae";
+$profileIndex{"Q"} = "raphidophyceae";
+$profileIndex{"R"} = "rhizaria";
+$profileIndex{"S"} = "synurophyceae";
+$profileIndex{"T"} = "tracheophyta";
+$profileIndex{"U"} = "eustigmatophyceae";
+$profileIndex{"V"} = "undefined";
+$profileIndex{"W"} = "undefined";
+$profileIndex{"X"} = "apusozoa";
+$profileIndex{"Y"} = "parabasalia";
+$profileIndex{"Z"} = "undefined";
+
+
+@profileList = split(',',uc($type));   # Get the list of profile types
+foreach $entry (@profileList) {   # Go through the entered types
+  if (($entry eq "ALL") || ($entry eq ".")) {   # If "all" among the entries
+    push(@profileSet,"A");   # Add the alveolates profiles to the investigation set
+    push(@profileSet,"B");   # Add the bryophytes profiles to the investigation set
+    push(@profileSet,"C");   # Add the bacillariophyta profiles to the investigation set
+    push(@profileSet,"D");   # Add the amoebozoa profiles to the investigation set
+    push(@profileSet,"E");   # Add the euglenozoa profiles to the investigation set
+    push(@profileSet,"F");   # Add the fungi profiles to the investigation set
+    push(@profileSet,"G");   # Add the green algae profiles to the investigation set
+    push(@profileSet,"H");   # Add the red algae profiles to the investigation set
+    push(@profileSet,"I");   # Add the brown algae profiles to the investigation set
+    push(@profileSet,"L");   # Add the liverworts profiles to the investigation set
+    push(@profileSet,"M");   # Add the metazoa profiles to the investigation set
+    #push(@profileSet,"N");   # Add the microsporidia profiles to the investigation set
+    push(@profileSet,"O");   # Add the oomycetes profiles to the investigation set
+    push(@profileSet,"P");   # Add the prymnesiophytes profiles to the investigation set
+    push(@profileSet,"Q");   # Add the raphidophytes profiles to the investigation set
+    push(@profileSet,"R");   # Add the rhizaria profiles to the investigation set
+    push(@profileSet,"S");   # Add the synurophyceae profiles to the investigation set
+    push(@profileSet,"T");   # Add the tracheophyta (higher plants) profiles to the investigation set
+    push(@profileSet,"U");   # Add the eustigmatophytes profiles to the investigation set
+    #push(@profileSet,"X");   # Add the apusozoa profiles to the investigation set
+    #push(@profileSet,"Y");   # Add the parabasalia profiles to the investigation set
+  } else {
+    if (length($entry) == 1) {  # If the name has only one character
+      if ($entry =~ m/[ABCDERFHILMNOPQRSTUXY]/) {  # If the selected set exists
+	push(@profileSet,$entry);   # Add the selected profiles to the investigation set
+      }
+    } else {
+      if ($entry =~ m/ALVEOL/) {
+	push(@profileSet,"A");   # Add the alveolates profiles to the investigation set
+      }
+      if ($entry =~ m/BRYO/) {
+	push(@profileSet,"B");   # Add the bryophytes profiles to the investigation set
+      }
+      if ($entry =~ m/MOSS/) {
+	push(@profileSet,"B");   # Add the bryophytes profiles to the investigation set
+      }
+      if ($entry =~ m/BACILL/) {
+	push(@profileSet,"C");   # Add the bacillariophyta profiles to the investigation set
+      }
+      if ($entry =~ m/DIATOM/) {
+	push(@profileSet,"C");   # Add the bacillariophyta profiles to the investigation set
+      }
+      if ($entry =~ m/AMOEB/) {
+	push(@profileSet,"D");   # Add the amoebozoa profiles to the investigation set
+      }
+      if ($entry =~ m/EUGLE/) {
+	push(@profileSet,"E");   # Add the euglenozoa profiles to the investigation set
+      }
+      if ($entry =~ m/FUNG/) {
+	push(@profileSet,"F");   # Add the fungi profiles to the investigation set
+      }
+      if ($entry =~ m/GREEN/) {
+	push(@profileSet,"G");   # Add the green algae profiles to the investigation set
+      }
+      if ($entry =~ m/CHLORO/) {
+	push(@profileSet,"G");   # Add the green algae profiles to the investigation set
+      }
+      if ($entry =~ m/RED-AL/) {
+	push(@profileSet,"H");   # Add the red algae profiles to the investigation set
+      }
+      if ($entry =~ m/RHODO/) {
+	push(@profileSet,"H");   # Add the red algae profiles to the investigation set
+      }
+      if ($entry =~ m/BROWN/) {
+	push(@profileSet,"I");   # Add the brown algae profiles to the investigation set
+      }
+      if ($entry =~ m/PHAEOP/) {
+	push(@profileSet,"I");   # Add the brown algae profiles to the investigation set
+      }
+      if ($entry =~ m/LIVER/) {
+	push(@profileSet,"L");   # Add the liverworts profiles to the investigation set
+      }
+      if ($entry =~ m/MARCH/) {
+	push(@profileSet,"L");   # Add the liverworts profiles to the investigation set
+      }
+      if ($entry =~ m/METAZ/) {
+	push(@profileSet,"M");   # Add the metazoa profiles to the investigation set
+      }
+      if ($entry =~ m/ANIMAL/) {
+	push(@profileSet,"M");   # Add the metazoa profiles to the investigation set
+      }
+      if ($entry =~ m/MICROSPOR/) {
+	push(@profileSet,"N");   # Add the microsporidia profiles to the investigation set
+      }
+      if ($entry =~ m/OOMYC/) {
+	push(@profileSet,"O");   # Add the oomycetes profiles to the investigation set
+      }
+      if ($entry =~ m/PRYMN/) {
+	push(@profileSet,"P");   # Add the prymnesiophytes profiles to the investigation set
+      }
+      if ($entry =~ m/HAPTO/) {
+	push(@profileSet,"P");   # Add the prymnesiophytes profiles to the investigation set
+      }
+      if ($entry =~ m/RAPHID/) {
+	push(@profileSet,"Q");   # Add the raphidophytes profiles to the investigation set
+      }
+      if ($entry =~ m/RHIZA/) {
+	push(@profileSet,"R");   # Add the rhizaria profiles to the investigation set
+      }
+      if ($entry =~ m/SYNUR/) {
+	push(@profileSet,"S");   # Add the synurophyceae profiles to the investigation set
+      }
+      if ($entry =~ m/TRACHE/) {
+	push(@profileSet,"T");   # Add the tracheophyta profiles to the investigation set
+      }
+      if ($entry =~ m/PLANTS/) {
+	push(@profileSet,"T");   # Add the tracheophyta profiles to the investigation set
+      }
+      if ($entry =~ m/EUSTIG/) {
+	push(@profileSet,"U");   # Add the eustigmatophytes profiles to the investigation set
+      }
+      if ($entry =~ m/APUSO/) {
+	push(@profileSet,"X");   # Add the apusozoa profiles to the investigation set
+      }
+      if ($entry =~ m/PARAB/) {
+	push(@profileSet,"Y");   # Add the parabasalia profiles to the investigation set
+      }
+    }
+  }
+}
+
+foreach $set (@profileSet) {   # For each set of profiles in the the full profile set for investigation
+  $hmmPath = $profileDB . "/" . $set . ".hmm";   # Determine the path to the HMM-file
+  chomp($modelCount = `grep -c "//" $hmmPath`);   # Count the number of models in the HMM-file
+  push(@modelCount,$modelCount);   # Add the number of models in this HMM-file to the list of model counts
+  if ($reset == 1) {
+    `rm -f $hmmPath.h3* 2> /dev/null`;   # Delete old HMM-files
+    `hmmpress $hmmPath 2> /dev/null`;   # Prepare the HMM-file for searching
+    ## Redirecting stderr is a quick and dirty solution to get rid of the messages... Could be made more elegant
+  }
+}
+
+## Clean-up input files and create complementary strand if needed
+## Get the current time and output an info message
+$now = localtime;
+if ($silent == 0) {
+  print STDERR "$now : Checking and handling input sequence data (should not take long)...\n";
+}
+
+## Open the summary file for writing
+if ($out_sum == 1) {   # If summary output is on
+  open (SUMMARY, ">>$output.summary.txt");   # Append to the summary file
+}
+
+if ($input ne "") {   # If an input file is given
+  ## Read from file
+  open (SEQUENCES, $input);   # Open the input file for reading
+  open (MAIN, ">$tempDir/main.fasta");   # Create a temporary file for storing the cleaned sequences representing the main strand
+  open (COMPLEMENT, ">$tempDir/complement.fasta");   # Create a temporary file for storing the cleaned sequences representing the complementary strand
+  $inputSequenceCount = 0;   # Reset input sequence counter
+  while ($sequence = <SEQUENCES>) {   # Repeat for every line in the input file
+    chomp($sequence);   # Truncate any potential line feeds
+    $sequence =~ tr/\r\n//d;  # Remove all carriage return and new line characters
+    if (substr($sequence,0,1) eq ">") {   # If a new FASTA entry is found in the input file
+      $inputSequenceCount++;   # Add one to the input sequence counter
+      print MAIN $mainSeq . "\n";   # Write the previous main DNA sequence to the main sequence file
+      print MAIN $sequence . " main\n";   # Write the definition line of the new sequence to the main sequence file
+      $sequenceDB{"$sequenceID"} = $mainSeq;  # Add sequence to sequence database
+      $headers{"$sequenceID"} = $header;  # Add the header to the header database
+      push(@sequenceOrder,$sequenceID);  # Add this sequence ID to the ordered list of sequences
+      ($sequenceID) = split(" ",substr($sequence,1));
+      $header = $sequence;  # Save the sequence header
+      $mainSeq = "";   # Empty the main sequence entry
+      if ($complement == 1) {   # If the complementary file should be written
+	$complementSeq = reverse($complementSeq);   # Reverse the complementary DNA sequence
+	print COMPLEMENT $complementSeq . "\n";   # Write the previous complementary DNA sequence to the complementary sequence file
+	print COMPLEMENT $sequence . " complement\n";   # Write the definition line of the new sequence to the complement sequence file
+	$complementSeq = "";   # Empty the complementary sequence entry
+      }
+    } else {    # If this line is just a continuation of the current DNA sequence
+      $mseq = $sequence;   # Store this part of the DNA sequence in the intermediate varaible $mseq
+      $mseq =~ s/[ .-]//g;   # Remove any gap characters present in this sequence part (good if input was an alignment)
+      $mseq =~ tr/[a-z]/[A-Z]/;   # Make all letters uppercase
+      $mseq =~ s/[^A-Z]//g;  # Remove all non-alphabetic characters
+      $mseq =~ tr/U/T/;   # Exchanges U:s for T:s (Uracil to Thymine, good if input was RNA sequence)
+      $mainSeq = $mainSeq . $mseq;   # Add the intermediate DNA sequence to the end of the main DNA sequence entry
+      
+      if ($complement == 1) {   # If the complementary file should be written
+	$cseq = $sequence;   # Store this part of the DNA sequence in the intermediate varaible $cseq
+	$cseq =~ s/[ .-]//g;   # Remove any gap characters present in this sequence part (good if input was an alignment)
+	$cseq =~ tr/[a-z]/[A-Z]/;   # Make all letters uppercase
+	$cseq =~ s/[^A-Z]//g;  # Remove all non-alphabetic characters
+	$cseq =~ tr/ACGTURYSWKMBDHVN/TGCAAYRWSMKVHDBN/;   # Replace all characters with its complementary base
+	$complementSeq = $complementSeq . $cseq;   # Add the intermediate DNA sequence to the end of the complementary DNA sequence entry
+      }
+    }
+  }
+  ## When the input file's end is reached
+  print MAIN $mainSeq . "\n";   # Write the last main DNA sequence to the main sequence file
+  if ($complement == 1) {   # If the complementary file should be written
+    $complementSeq = reverse($complementSeq);   # Reverse the complementary DNA sequence
+    print COMPLEMENT $complementSeq . "\n";   # Write the last complementary DNA sequence to the complementary sequence file
+  }
+  $sequenceDB{"$sequenceID"} = $mainSeq;  # Add sequence to sequence database
+  $headers{"$sequenceID"} = $header;  # Add the header to the header database
+  push(@sequenceOrder,$sequenceID);  # Add this sequence ID to the ordered list of sequences
+  $mainSeq = "";   # Empty the main sequence entry
+  $complementSeq = "";   # Empty the complementary sequence entry
+  close (SEQUENCES);   # Close the sequence input file
+  close (COMPLEMENT);   # Close the complementary output file
+  close (MAIN);   # Close the main output file
+} else {   # If no input file is supplied, then read from stdin instead
+  $input = "$tempDir/main.fasta";   # Set up a temporary input file path
+  open (MAIN, ">$tempDir/main.fasta");   # Create a temporary file for storing the cleaned sequences representing the main strand
+  open (COMPLEMENT, ">$tempDir/complement.fasta");   # Create a temporary file for storing the cleaned sequences representing the complementary strand
+  $inputSequenceCount = 0;   # Reset input sequence counter
+  while ($sequence = <STDIN>) {   # Repeat for every line in the standard input
+    chomp($sequence);   # Truncate any potential line feeds
+    if (substr($sequence,0,1) eq ">") {   # If a new FASTA entry is found in the input
+      $inputSequenceCount++;   # Add one to the input sequence counter
+      print MAIN $mainSeq . "\n";   # Write the previous main DNA sequence to the main sequence file
+      print MAIN $sequence . " main\n";   # Write the definition line of the new sequence to the main sequence file
+      $sequenceDB{"$sequenceID"} = $mainSeq;  # Add sequence to sequence database
+      $headers{"$sequenceID"} = $header;  # Add the header to the header database
+      push(@sequenceOrder,$sequenceID);  # Add this sequence ID to the ordered list of sequences
+      $sequenceID = split(" ",substr($sequence,1));
+      $header = $sequence;  # Save the sequence header
+      $mainSeq = "";   # Empty the main sequence entry
+      if ($complement == 1) {   # If the complementary file should be written
+	$complementSeq = reverse($complementSeq);   # Reverse the complementary DNA sequence
+	print COMPLEMENT $complementSeq . "\n";   # Write the previous complementary DNA sequence to the complementary sequence file
+	print COMPLEMENT $sequence . " complement\n";   # Write the definition line of the new sequence to the complement sequence file
+	$complementSeq = "";   # Empty the complementary sequence entry
+      }
+    } else {    # If this line is just a continuation of the current DNA sequence
+      $mseq = $sequence;   # Store this part of the DNA sequence in the intermediate varaible $mseq
+      $mseq =~ s/[ .-]//g;   # Remove any gap characters present in this sequence part (good if input was an alignment)
+      $mseq =~ tr/[a-z]/[A-Z]/;   # Make all letters uppercase
+      $mseq =~ s/[^A-Z]//g;  # Remove all non-alphabetic characters
+      $mseq =~ tr/U/T/;   # Exchanges U:s for T:s (Uracil to Thymine, good if input was RNA sequence)
+      $mainSeq = $mainSeq . $mseq;   # Add the intermediate DNA sequence to the end of the main DNA sequence entry
+      
+      if ($complement == 1) {   # If the complementary file should be written
+	$cseq = $sequence;   # Store this part of the DNA sequence in the intermediate varaible $cseq
+	$cseq =~ s/[ .-]//g;   # Remove any gap characters present in this sequence part (good if input was an alignment)
+	$cseq =~ tr/[a-z]/[A-Z]/;   # Make all letters uppercase
+	$cseq =~ s/[^A-Z]//g;  # Remove all non-alphabetic characters
+	$cseq =~ tr/ACGTURYSWKMBDHVN/TGCAAYRWSMKVHDBN/;   # Replace all characters with its complementary base
+	$complementSeq = $complementSeq . $cseq;   # Add the intermediate DNA sequence to the end of the complementary DNA sequence entry
+      }
+    }
+  }
+  ## When the input file's end is reached
+  print MAIN $mainSeq . "\n";   # Write the last main DNA sequence to the main sequence file
+  if ($complement == 1) {   # If the complementary file should be written
+    $complementSeq = reverse($complementSeq);   # Reverse the complementary DNA sequence
+    print COMPLEMENT $complementSeq . "\n";   # Write the last complementary DNA sequence to the complementary sequence file
+  }
+  $sequenceDB{"$sequenceID"} = $mainSeq;  # Add sequence to sequence database
+  $headers{"$sequenceID"} = $header;  # Add the header to the header database
+  push(@sequenceOrder,$sequenceID);  # Add this sequence ID to the ordered list of sequences
+  $mainSeq = "";   # Empty the main sequence entry
+  $complementSeq = "";   # Empty the complementary sequence entry
+  close (COMPLEMENT);   # Close the complementary output file
+  close (MAIN);   # Close the main output file
+}
+
+if ($out_sum == 1) {   # If summary output should be written
+  print SUMMARY "Number of sequences in input file:       \t$inputSequenceCount\n";   # Write info on the number of input sequences to the summary file
+}
+
+
+## Perform HMM-scan
+if ($hmmscan eq "") {   # If a pre-computed hmmscan output file is not supplied
+  if ($heuristics == 0) {  # If HMMER's heuristic filtering should not be used
+    $heurMax = "--max";  # Set the heurMax to "--max" (indicating that HMMER should turn off filtering)
+  } else {
+    $heurMax = "";  # Set the heurMax to empty (indicating that HMMER should turn on filtering)
+  }
+  if ($multi_thread == 0) {   # If multi-threading is off
+    ## Get the current time and output info message
+    $now = localtime;
+    if ($silent == 0) {
+      print STDERR "$now : Comparing sequences to HMM database (this may take a long while)...\n";
+    }
+    foreach $set (@profileSet) {   # Go sequentially through all profile sets to search for
+      $hmmPath = $profileDB . "/" . $set . ".hmm";   # Set the path to the HMM-file of the current set
+      if ($search_eval ne "") {   # If E-value cutoff is use for the search
+	hmmerSearch("hmmscan --cpu $cpu $heurMax -E $search_eval $hmmPath $tempDir/main.fasta","$tempDir/main.$set.hmmscan","M",$set);   # Call HMMER with E-value cutoff
+      } else {   # If score cutoff is use for the search
+	hmmerSearch("hmmscan --cpu $cpu $heurMax -T $search_score $hmmPath $tempDir/main.fasta","$tempDir/main.$set.hmmscan","M",$set);   # Call HMMER with score cutoff
+      }
+      if ($complement == 1) {   # If the complementary file should be scanned
+	if ($search_eval ne "") {   # If E-value cutoff is use for the search
+	  hmmerSearch("hmmscan --cpu $cpu $heurMax -E $search_eval $hmmPath $tempDir/complement.fasta","$tempDir/complement.$set.hmmscan","C",$set);   # Call HMMER with E-value cutoff
+	} else {   # If score cutoff is use for the search
+	  hmmerSearch("hmmscan --cpu $cpu $heurMax -T $search_score $hmmPath $tempDir/complement.fasta","$tempDir/complement.$set.hmmscan","C",$set);   # Call HMMER with score cutoff
+	}
+      }
+    }
+  } else {   # If multi-threading is on
+    ## Get the current time and output info message
+    $now = localtime;
+    if ($silent == 0) {
+      print STDERR "$now : Doing paralellised comparison to HMM database (this may take a long while)...\n";
+    }
+
+    ## Determining number of cpus per thread
+    if ($complement == 1) {   # If the complementary file should be scanned
+      $hmmcpu = int(0.5 * $cpu / scalar(@profileSet));   # Assign X CPUs to each thread, X = 0.5 * (TOTAL_CPUs_USED) / (TOTAL_NUMBER_OF_PROFILE_SETS)
+    } else {   # If the complementary file should not be scanned
+      $hmmcpu = int($cpu / scalar(@profileSet));   # Assign X CPUs to each thread, X = (TOTAL_CPUs_USED) / (TOTAL_NUMBER_OF_PROFILE_SETS)
+    }
+    if ($hmmcpu < 1) {   # If the number of CPUs per thread is smaller than 1
+      $hmmcpu = 1;   # Give each thread at least one CPU to work on
+    }
+    ## Main strand searches...
+    $cpuCount = 0;
+    foreach $set (@profileSet) {   # Go through each profile set to investigate
+      if ($cpuCount < $cpu) {
+	$cpuCount++;
+	$pid = fork();   # Fork off a copy of this process for this set
+      } else {
+	$deceasedPID = wait();   # Wait until a PID is finished, and gather its number
+	$pid = fork();   # Fork off a copy of this process for this set
+      }
+      if ($pid != 0) {   # If this is the parent process
+	push(@pids,$pid);   # Add the new process ID to the list of active process IDs
+      } else {   # If this is the new child process
+	$hmmPath = $profileDB . "/" . $set . ".hmm";   # Set the path to the HMM-file of the current set
+	if ($search_eval ne "") {   # If E-value cutoff is use for the search
+	  hmmerSearch("hmmscan --cpu $hmmcpu $heurMax -E $search_eval $hmmPath $tempDir/main.fasta","$tempDir/main.$set.hmmscan","M",$set);   # Call HMMER with E-value cutoff
+	} else {   # If score cutoff is use for the search
+	  hmmerSearch("hmmscan --cpu $hmmcpu $heurMax -T $search_score $hmmPath $tempDir/main.fasta","$tempDir/main.$set.hmmscan","M",$set);   # Call HMMER with score cutoff
+	}
+	## Stop child process...
+	exit;   # Exits the child process
+      }
+    }
+    ## Revese strand searches...
+    if ($complement == 1) {   # If the complementary file should be scanned
+      foreach $set (@profileSet) {   # Go through each profile set to investigate
+	if ($cpuCount < $cpu) {
+	  $cpuCount++;
+	  $pid = fork();   # Fork off a copy of this process for this set
+	} else {
+	  $deceasedPID = wait();   # Wait until a PID is finished, and gather its number
+	  $pid = fork();   # Fork off a copy of this process for this set
+	}
+	if ($pid != 0) {   # If this is the parent process
+	  push(@pids,$pid);   # Add the new process ID to the list of active process IDs
+	} else {   # If this is the new child process
+	  $hmmPath = $profileDB . "/" . $set . ".hmm";   # Set the path to the HMM-file of the current set
+	  if ($search_eval ne "") {   # If E-value cutoff is use for the search
+	    hmmerSearch("hmmscan --cpu $hmmcpu $heurMax -E $search_eval $hmmPath $tempDir/complement.fasta","$tempDir/complement.$set.hmmscan","C",$set);   # Call HMMER with E-value cutoff
+	  } else {   # If score cutoff is use for the search
+	    hmmerSearch("hmmscan --cpu $hmmcpu $heurMax -T $search_score $hmmPath $tempDir/complement.fasta","$tempDir/complement.$set.hmmscan","C",$set);   # Call HMMER with score cutoff
+	  }
+	  ## Stop child process...
+	  exit;   # Exits the child process
+	}
+      }
+    }
+    ## Get the current time and output the active process IDs
+    $now = localtime;
+    #print STDERR "    $now : Active PIDs: ";
+    #foreach $p (@pids) {   # Go through the list of PIDs
+    #  print STDERR "$p ";   # Print the PID
+    #}
+    #print STDERR "\n";   # Print a new line
+    do {   # Loop until all child PIDs have finished.
+      $deceasedPID = wait();   # Wait until a PID is finished, and gather its number
+      $now = localtime;   # Get the current time
+      if ($deceasedPID > -1) {   # If the PID that finished wasn't the last active one
+#	print STDERR "    $now : PID $deceasedPID finished.\n";   # Print finished PID
+      } else {   # If PID that finished was the last
+	print STDERR "    $now : All processes finished.\n";   # Print that all PIDs have finished
+      }
+    } until (wait() == -1);   # Do this loop until all PIDs have finished
+    $now = localtime;   # Get current time
+    print STDERR "$now : Parallel HMM-scan finished.\n";   # Print informative finishing message
+  }
+} else {    # If a pre-computed hmmscan file is supplied then
+## Get the current time and output that the hmmscan step is skipped
+  $now = localtime;
+  if ($silent == 0) {
+    print STDERR "$now : Skipping hmmscan! Using $hmmscan as input for the analysis instead.\n";
+  }
+}
+
+## Analyse HMM-scan output
+## Get the current time and output info
+$now = localtime;
+if ($silent == 0) {
+  print STDERR "$now : Analysing results of HMM-scan (this might take quite some time)...\n";
+}
+
+## Set up output files
+if ($out_table == 1) {   # If table output is on
+  open (TABLE, ">$output.hmmer.table");   # Create a table output file
+}
+if ($out_graph == 1) {   # If graphical output is on
+  open (GRAPH, ">$output.graph");   # Create a graph output file
+}
+if ($out_not == 1) {   # If not-found output is on
+  open (NOTFOUND, ">$tempDir/$output\_hmmer_no_detections.txt");   # Create a HMMER not-found output file
+}
+
+$setI = 0;   # Set the profile set indicator to zero
+foreach $set (@profileSet) {   # Go through all the profile sets to be investigated
+  for ($co = 0; $co <= 1; $co++) {   # Do main (and complementary) strand analysis in order
+    if ($co > 0) {   # If main strand analysis is finished
+      if ($complement == 1) {   # If complementary strand should be analysed
+	open (HMMOUTPUT, "$tempDir/complement.$set.hmmscan");   # Open hmmscan output for reading
+	open (SEQUENCES, "$tempDir/complement.fasta");   # Open complementary sequence file for reading
+	if ($out_table == 1) {   # If table output is on, write a header for this set
+	  print TABLE "***********************************************************\n";
+	  print TABLE "$set matches on complementary strand:\n";
+	}
+	if ($out_graph == 1) {   # If graphical output is on, write a header for this set
+	  print GRAPH "***********************************************************\n";
+	  print GRAPH "$set matches on complementary strand:\n";
+	}
+      } else {   # If complementary strand should not be analysed
+	last;   # Exit this loop
+      }
+    }
+    if ($co == 0) {   # If main strand analysis is not finished
+      open (HMMOUTPUT, "$tempDir/main.$set.hmmscan");   # Open hmmscan output for reading
+      open (SEQUENCES, "$tempDir/main.fasta");   # Open main sequence file for reading
+      if ($out_table == 1) {   # If table output is on, write a header for this set
+	print TABLE "***********************************************************\n";
+	print TABLE "$set matches on main strand:\n";
+      }
+      if ($out_graph == 1) {   # If graphical output is on, write a header for this set
+	print GRAPH "***********************************************************\n";
+	print GRAPH "$set matches on main strand:\n";
+      }
+    }
+  ## Read and analyse hmmscan output file
+  while ($line = <HMMOUTPUT>) {   # Read in the hmmscan output file, line by line
+    chomp($line);   # Remove any potential line feeds
+    if (substr($line,0,13) eq "## New query:") {   # If this line begin with "## New query:", then this is a new entry
+      undef %hits;   # Empty the hits hash
+      undef %evals;   # Empty the e-value hash
+      undef %scores;   # Empty the score hash
+      $querytemp = substr($line,14);   # Extract everything from this line, except for the start ("## New query:")
+      ($query,$length) = split('\t',$querytemp);  # Split query name, length and DNA sequence
+      if ($co == 0) {   # If main strand analysis
+	$DNA = $sequenceDB{"$query"};
+      } else {  # If complementary strand
+	$cseq = $sequenceDB{"$query"};   # Store the DNA sequence in the intermediate varaible $cseq
+	$cseq =~ tr/[a-z]/[A-Z]/;   # Make all letters uppercase
+	$cseq =~ s/[^A-Z]//g;  # Remove all non-alphabetic characters
+	$cseq =~ tr/ACGTURYSWKMBDHVN/TGCAAYRWSMKVHDBN/;   # Replace all characters with its complementary base
+	$DNA = reverse($cseq);
+      }
+    } else {
+      
+      if ($line ne "//") {
+	## Find domain annotations...
+	($query,$matchProfile,$length,$domNo,$sign,$score,$bias,$cE,$iE,$hmmFrom,$hmmTo,$hmmends,$queryFrom,$queryTo,$queryends,$envFrom,$envTo,$envends,$acc) = split('\t',$line);   # Split the line into a collection of stat variables
+
+	$useQueryFrom = $queryFrom;
+	if ($hmmFrom > 1) {  # If the HMM-profile is not matched from the beginning
+	  $hmmDiff = $hmmFrom - 1;
+	  if ($useQueryFrom > $hmmDiff) {
+	    $useQueryFrom = $queryFrom - $hmmDiff;
+	  } else {
+	    $useQueryFrom = 1;
+	  }
+	}
+
+	if (uc($anchor) eq "HMM") {
+	  $anchorLen = $hmmTo - $hmmFrom + 1;
+	} else {
+	  $anchorLen = $anchor;
+	}
+	  
+
+	$query_profile_match = $query . ":" . $matchProfile;
+	$profileExists = 0;   # Assume that the newly found match profile ($matchProfile) is not already found for this sequence
+	if (exists($hits{$query_profile_match})) {   # If a profile from the list is the same as the match profile
+	  ($hitFrom,$hitTo,$hitProfile,$hitScore,$hitE) = split('\t',$hits{$query_profile_match});   # Split the entry in list into stat variables
+	  if ($iE < $hitE) {   # If the new match profile has a smaller E-value than the one from the list
+	    $hits{$query_profile_match} = "$useQueryFrom\t$envTo\t$matchProfile\t$score\t$iE\t$anchorLen";   # Replace the data in the hit list with the data for the newly found match profile
+	    $evals{$query_profile_match} = $iE;   # Replace the E-value in the hit list with the E-value for the newly found match profile
+	    $scores{$query_profile_match} = $score;   # Replace the score in the hit list with the score for the newly found match profile
+	  }
+	  $profileExists = 1;   # Indicate that this match profile was found in the hit list
+	}
+	if ($profileExists == 0) {   # If the match profile was not found in the hit list
+	  if (($iE <= $E) && ($score >= $S)) {  # If this hits lives up to the minimal score and E-value cutoffs
+	    $hits{$query_profile_match} = "$useQueryFrom\t$envTo\t$matchProfile\t$score\t$iE\t$anchorLen";   # Add the data for the newly found match profile
+	    $evals{$query_profile_match} = $iE;   # Add the E-value for the newly found match profile
+	    $scores{$query_profile_match} = $score;   # Add the score for the newly found match profile
+	  }
+	}
+      } else {   # If the line only contains "//", the end of this sequence's hmmscan entry is reached
+      ## Save analysis results
+	@sortedKeys = sort {$hits{$a} <=> $hits{$b}} keys(%hits);   # Sort the the list of hits numerically ascending (smallest first)
+	undef @sortedHits;
+	undef @scores;
+	undef @evals;
+	foreach $key (@sortedKeys) {
+	  push(@sortedHits, $hits{$key});   # Add the hit to the list of hits numerically ascending (smallest first)
+	  push(@scores, $scores{$key});  # Add the score to the scores array
+	  push(@evals, $evals{$key});  # Add the E-value to the evals array
+	}
+
+	## If the number of hits > N, the min eval < E and the max score > S then include query sequence
+	## OR if a single domain satisfies the thresholds and this is allowed, include it!
+	if ( ((scalar(@sortedHits) >= $N) && (min(@evals) <= $E) && (max(@scores) >= $S)) ||
+	     ((scalar(@sortedHits) > 0) && ($allow_single_E >= 0) && (min(@evals) <= $allow_single_E) && (max(@scores) >= $allow_single_score)) ) {
+	  if ($debug == 1) {   # If debugging mode is on
+	    print STDERR $query . " :\t" . scalar(@sortedHits) . "\t" . min(@evals) . "\t" . max(@scores) . "\n";   # Print some top hit statistics
+	  }
+	  
+	  ## Save some total stats to be able to determine origin of ITS sequence
+	  if (scalar(@evals) > 0) {   # If there are any E-values stored
+	    $averageE = sum(@evals) / scalar(@evals);   # Calculate the average E-value for this profile set
+	    $averageScore = sum(@scores) / scalar(@scores);   # Calculate the average score for this profile set
+	    $numberOfDomains = scalar(@sortedHits);   # Calculate the number of domains matched on this sequence
+	    
+	    #$scoreSum = sum(@scores) / @modelCount[$setI];   # Calculate score sum as: sum / (no. of profiles of this given type)
+	    $scoreSum = sum(@scores) / 4;   # Calculate score sum as: sum / (no. of profiles of this given type)
+	    
+	    $saveThis = "$query\t$set\t$co\t$numberOfDomains\t$averageE\t$averageScore\t$scoreSum\t$DNA\t";   # Collect the variables to save for this sequence and this profile set
+	    foreach $hit (@sortedHits) {   # Go through the list of hits and add specific information to save from each hit
+	      ($hitFrom,$hitTo,$hitProfile,$hitScore,$hitE,$hitanchorlen) = split('\t',$hit);   # Extract information from this hit
+	      $saveThis = $saveThis . "$hitFrom;$hitTo;$hitProfile;$hitScore;$hitE;$hitanchorlen\t";   # Add information to the list of variables to save
+	    }
+	    #push(@allHits, $saveThis);   # Add this information to the collection of all hits for this sequence, across all profile sets
+	    if (exists($allHits{$query})) {
+	      $allHits{$query} = $allHits{$query} . "\n" . $saveThis;
+	    } else {
+	      $allHits{$query} = $saveThis;
+	    }
+	  }
+	  
+	  if ($out_table == 1) {   # If table output is on
+	    print TABLE $query . "\t" . $length . "\t";   # Print query and length information to table
+	    foreach $hit (@sortedHits) {   # Go through each hit in the hit list
+	      ($hitFrom,$hitTo,$hitProfile,$hitScore,$hitE,$hitanchorlen) = split('\t',$hit);   # Extract data corresponding to this hit
+	      print TABLE "$hitFrom - $hitTo: $hitProfile ($hitScore, $hitE)\t";   # Print hit information to table
+	    }
+	    print TABLE "\n";   # Print new line
+	  }
+	  
+	  if ($out_graph == 1) {   # If graphical output is on
+	    print GRAPH ">> " . $query . "\t" . $length . " bp\n";   # Print a sequence header
+	    $insertPoint = 0;   # Set the domain insert point to beginning of line
+	    $hi = 0;   # Set hit number to zero
+	    foreach $hit (@sortedHits) {   # Go through the hit list
+	      ($hitFrom,$hitTo,$hitProfile,$hitScore,$hitE,$anchorLen) = split('\t',$hit);   # Split the hit into stat variables
+	      if ($graph_scale == 0) {   # If the graph scale is scaled individually to 100% for each sequence
+		$pFrom = $hitFrom / $length * 100;   # Set the profile start on graph relative to its position in the sequence
+		$pTo = $hitTo / $length * 100;   # Set the profile end on graph relative to its position in the sequence
+		$pEnd = 100;   # Set the end of the sequence graph to be at 100
+	      } else {   # If the scale is the same for all sequences
+		$pFrom = $hitFrom * $graph_scale;   # Set the profile start on graph scaled to the parameter given
+		$pTo = $hitTo * $graph_scale;   # Set the profile end on graph scaled to the parameter given
+		$pEnd = $length * $graph_scale;   # Set the end of the sequence graph to be at the end of the sequence scaled to the parameter given
+	      }
+	      for ($insertPoint = $insertPoint; $insertPoint <= $pFrom; $insertPoint++) {   # Go forward through the sequence, moving the insert point one step at a time until the beginning of the next profile is reached
+		print GRAPH "-";   # Print a "-"
+	      }
+	      print GRAPH substr($hitProfile,2,3);   # When the profile is reached, print its name
+	      $insertPoint = $insertPoint + 3;   # Move the insert point three steps forward, to account for the inserted name
+	      ($nextHitStart,$nextHitEnd,$nextProfile) = split('\t',@sortedHits[$hi + 1]);   # Check where the next hit in the list is located
+	      if (($nextHitStart <= $hitTo) && ($nextHitStart > 0)) {   # If the next hit in the list overlaps with this profile
+		if ($graph_scale == 0) {   # If the scale is relative
+		  $pTo = $nextHitStart / $length * 100 - 1;   # Change the profile end on the graph to be where this next profile starts
+		} else {   # If the scale is the same for all sequences
+		  $pTo = $nextHitStart * $graph_scale - 1;   # Change the profile end on the graph to be where this next profile starts
+		}
+	      }
+	      for ($insertPoint = $insertPoint; $insertPoint <= $pTo; $insertPoint++) {   # Go forward through the sequence, moving the insert point one step at a time until the end of the current profile is reached
+		print GRAPH "=";   # Print a "="
+	      }
+	      if (($nextHitStart <= $hitTo) && ($nextHitStart > 0)) {   # If the next hit in the list overlaps with this profile
+		print GRAPH ">";   # Print a ">" to indicate the profile overlap
+		$insertPoint++;   # Move the insert point one additional step forward to account for the ">" inserted
+	      }
+	      $hi++;   # Increase the hit number by one
+	    }
+	    for ($insertPoint = $insertPoint; $insertPoint <= $pEnd; $insertPoint++) {   # If there is no more profile matches to sequence, go forward through the sequence, moving the insert point one step at a time until the end of the sequence is reached
+	      print GRAPH "-";   # Print a "-"
+	    }
+	    print GRAPH "\n";   # Print a new line, indicating the end of this sequence entry
+	  }
+	} else {   # If this sequence didn't find any good-enough profile matches
+	  if ($out_not == 1) {   # If not-found output is on
+	    print NOTFOUND $query . "\n";   # Print the name of this query to the not-found list
+	  }
+	}
+      }
+    }
+  }
+  close (SEQUENCES);   # Close the input sequence file
+  close (HMMOUTPUT);   # Close the hmmscan output file
+  }
+  $setI++;   # Add one to the profile set indicator
+}
+
+## Close output files
+if ($out_table == 1) {   # If table output is on, close the table file
+  close (TABLE);
+}
+if ($out_graph == 1) {   # If graphical output is on, close the graph file
+  close (GRAPH);
+}
+if ($out_not == 1) {   # If not-found output is on, close the not-found file
+  close (NOTFOUND);
+  #$profileCount = scalar(@profileSet);   # Count the number of profile sets
+  #if ($complement == 1) {   # If complementary strand was scanned
+  #  $profileCount = $profileCount * 2;   # Double the number of profile sets that was investigated (and thus the number of not-founds that could at max be found)
+  #}
+  #`sort $tempDir/hmmer_no_detections.txt | uniq -c | grep "  *$profileCount " | sed "s/  *$profileCount //" > $output\_no_detections.txt`;   # Sort the not-found list, count the number of profile sets having no matches for each query. Save those that have only non-matches to the hmmer-not-found file
+}
+
+
+## Create total collected output and FASTA output
+
+if ($out_results == 1) {
+  open (RESULTS, ">$output.extraction.results");   # Create a results file
+}
+open (RAWOUT, ">$tempDir/ITSx_output.raw");   # Create a raw output file for ALL data
+open (PROBLEM, ">$output.problematic.txt");   # Create a file for problematic entries
+$foundProblem = 0;
+if ($out_pos == 1) {
+  open (POS, ">$output.positions.txt");   # Create a positions file
+}
+if ($out_fasta == 1) {   # If FASTA output should be written
+  open (FASTA, ">$output.full.fasta");   # Create a FASTA output file for found sequences
+  if ($allow_reorder == 0) {   # If reordering of domains is not allowed
+    open (CHIMERA, ">$output.chimeric.fasta");   # Create a FASTA file for potential chimera sequences with profile matches in the wrong order
+    $foundChimera = 0;
+  }
+  if ($out_partial > 0) {
+    open (FULLPARTIAL, ">$output.full_and_partial.fasta");
+  }
+}
+
+if ($out_joined == 1) {   # If SSU FASTA output should be written
+  open (JOINED, ">$output.joined.fasta");
+}
+if ($out_ssu == 1) {   # If SSU FASTA output should be written
+  open (SSU, ">$output.SSU.fasta");
+}
+if ($out_lsu == 1) {   # If LSU FASTA output should be written
+  open (LSU, ">$output.LSU.fasta");
+}
+if ($out_58S == 1) {   # If 5.8S FASTA output should be written
+  open (MID, ">$output.5_8S.fasta");
+}
+if ($out_its1 == 1) {   # If ITS1 FASTA output should be written
+  open (ITS1, ">$output.ITS1.fasta");
+  if ($out_partial > 0) {
+    open (ITS1PARTIAL, ">$output.ITS1.full_and_partial.fasta");
+  }
+}
+if ($out_its2 == 1) {   # If ITS2 FASTA output should be written
+  open (ITS2, ">$output.ITS2.fasta");
+  if ($out_partial > 0) {
+    open (ITS2PARTIAL, ">$output.ITS2.full_and_partial.fasta");
+  }
+}
+if ($out_concat == 1) {  # If concatenated output should be written
+  open (CONCAT, ">$output.concat.fasta");
+}
+
+undef @sortedHits;   # Empty the array of sorted hits
+
+# @sortedHits = sort @allHits;   # Sort the full list of hits in alphabetical order (to be able to analyse all sequences with same ID at once)
+if ($out_not == 1) {   # If not-found output is on
+  open (NOTFOUND, ">$output\_no_detections.txt");   # Create a not-found output file
+  $noDetect = 0;
+}
+foreach $sequenceID (@sequenceOrder) {   # Sort the full list of hits in their original order (to be able to analyse all sequences with same ID at once)
+  if ($sequenceID ne "") {
+    $countsInList = 0;
+    if (exists($allHits{$sequenceID})) {
+      @allHits = split('\n',$allHits{$sequenceID});
+    } else {
+      undef @allHits;
+    }
+    foreach $line (@allHits) {
+      @item = split('\t',$line);   # Split the line into an array
+      if (@item[0] eq $sequenceID) {  # If this item corresponds to the current sequence ID
+	push(@sortedHits,$line);  # Add it to the sorted list of hits
+	$countsInList++;
+      }
+    }
+    if ($countsInList == 0) {  # If no matches were found
+      if ($out_not == 1) {   # If not-found output is on
+	print NOTFOUND "$sequenceID\n";  # Output the sequence ID
+	$noDetect++;
+      }
+    }
+  }
+}
+if ($out_not == 1) {   # If not-found output is on
+  close (NOTFOUND);
+  if ($noDetect == 0) {
+    `rm $output\_no_detections.txt 2> /dev/null`;
+  }
+}
+
+## Set all counts for different ITS types to zero
+undef @itsCounts;
+$itsChimeric = 0;
+$itsMain = 0;
+$itsCompl = 0;
+
+push(@sortedHits,"--END--");   # Add a last item to the sorted list, so that all items are securely saved
+$lc = 1;   # Set the line count to one
+foreach $line (@sortedHits) {   # Go through the list of found hits
+  print RAWOUT "$line\n";   # Write the raw data associated with this hit to the raw data output file
+  @item = split('\t',$line);   # Split the line into an array
+  ## If this sequence ID is the same as the saved ones, then add it, else empty the array sequence ID and save
+  if ((@seqID[0] ne @item[0]) && (@item[0] ne "") || ($lc > scalar(@sortedHits))) {   # If this sequence ID is not the same as the last one and is non-empty, or if the end of the list has been reached
+    ## Save profile-type which is most likely...
+    if ($priority eq "sum") {   # If the sum-of-scores algorithm should be used to determine the most likely profile-type
+      ## Reset variables to unrealisticly high or low values
+      $best = 0;
+      $bestCount = 0;
+      $bestEval = 1000;
+      $bestScore = -1000;
+      $bestSum = -1000;
+      
+      for ($i = 0; $i < scalar(@seqScoreSum); $i++) {   # Go through all sum-of-scores entries
+	if (@seqScoreSum[$i] > $bestSum) {   # If the current value is larger than the previous top value
+	  $best = $i;   # Set the best value to be the current value
+	  ## Set all other best variables to those corresponding to the current value
+	  $bestCount = @seqDomCounts[$i];
+	  $bestEval = @seqAvgE[$i];
+	  $bestScore = @seqAvgScore[$i];
+	  $bestSum = @seqScoreSum[$i];
+	}
+	if (@seqScoreSum[$i] == $bestSum) {   # If the current value is equal to the previous top value
+	  if (@seqDomCounts[$i] > $bestCount) {   # If the current domain count is larger than the previous top domain count
+	    $best = $i;   # Set the best value to be the current value
+	    ## Set all other best variables to those corresponding to the current value
+	    $bestCount = @seqDomCounts[$i];
+	    $bestEval = @seqAvgE[$i];
+	    $bestScore = @seqAvgScore[$i];
+	    $bestSum = @seqScoreSum[$i];
+	  }
+	  if (@seqDomCounts[$i] == $bestCount) {   # If the current domain count is equal to the previous top domain count
+	    if (@seqAvgE[$i] < $bestEval) {   # If the current E-value is smaller than the previous top E-value
+	      $best = $i;   # Set the best value to be the current value
+	      ## Set all other best variables to those corresponding to the current value
+	      $bestCount = @seqDomCounts[$i];
+	      $bestEval = @seqAvgE[$i];
+	      $bestScore = @seqAvgScore[$i];
+	      $bestSum = @seqScoreSum[$i];
+	    }
+	  }
+	}
+      }
+    }
+    if ($priority eq "domains") {   # If the number of found domains should be used to determine the most likely profile-type
+      ## Reset variables to unrealisticly high or low values
+      $best = 0;
+      $bestCount = 0;
+      $bestEval = 1000;
+      $bestScore = -1000;
+
+      for ($i = 0; $i < scalar(@seqDomCounts); $i++) {   # Go through all domain count entries
+	if (@seqDomCounts[$i] > $bestCount) {   # If the current domain count is larger than the previous top domain count
+	  $best = $i;   # Set the best value to be the current value
+	  ## Set all other best variables to those corresponding to the current value
+	  $bestCount = @seqDomCounts[$i];
+	  $bestEval = @seqAvgE[$i];
+	  $bestScore = @seqAvgScore[$i];
+	}
+	if (@seqDomCounts[$i] == $bestCount) {   # If the current domain count is equal to the previous top domain count
+	  if (@seqAvgE[$i] < $bestEval) {   # If the current E-value is smaller than the previous top E-value
+	    $best = $i;   # Set the best value to be the current value
+	    ## Set all other best variables to those corresponding to the current value
+	    $bestCount = @seqDomCounts[$i];
+	    $bestEval = @seqAvgE[$i];
+	    $bestScore = @seqAvgScore[$i];
+	  }
+	  if (@seqAvgE[$i] == $bestEval) {   # If the current E-value is equal to the previous top E-value
+	    if (@seqAvgScore[$i] > $bestScore) {   # If the current average score is larger than the previous top average score
+	      $best = $i;   # Set the best value to be the current value
+	      ## Set all other best variables to those corresponding to the current value
+	      $bestCount = @seqDomCounts[$i];
+	      $bestEval = @seqAvgE[$i];
+	      $bestScore = @seqAvgScore[$i];
+	    }
+	  }
+	}
+      }
+    }
+    if ($priority eq "eval") {   # If the average E-value should be used to determine the most likely profile-type
+      ## Reset variables to unrealisticly high or low values
+      $best = 0;
+      $bestCount = 0;
+      $bestEval = 1000;
+      $bestScore = -1000;
+
+      for ($i = 0; $i < scalar(@seqDomCounts); $i++) {   # Go through all domain counts entries
+	if (@seqAvgE[$i] < $bestEval) {   # If the current E-value is smaller than the previous top E-value
+	  $best = $i;   # Set the best value to be the current value
+	  ## Set all other best variables to those corresponding to the current value
+	  $bestCount = @seqDomCounts[$i];
+	  $bestEval = @seqAvgE[$i];
+	  $bestScore = @seqAvgScore[$i];
+	}
+	if (@seqAvgE[$i] == $bestEval) {   # If the current E-value is equal to the previous top E-value
+	  if (@seqAvgScore[$i] > $bestScore) {   # If the current average score is larger than the previous top average score
+	    $best = $i;   # Set the best value to be the current value
+	    ## Set all other best variables to those corresponding to the current value
+	    $bestCount = @seqDomCounts[$i];
+	    $bestEval = @seqAvgE[$i];
+	    $bestScore = @seqAvgScore[$i];
+	  }
+	  if (@seqAvgScore[$i] == $bestScore) {   # If the current average score is equal to the previous top average score
+	    if (@seqDomCounts[$i] > $bestCount) {   # If the current number of domains is larger than the previous top number of domains
+	      $best = $i;   # Set the best value to be the current value
+	      ## Set all other best variables to those corresponding to the current value
+	      $bestCount = @seqDomCounts[$i];
+	      $bestEval = @seqAvgE[$i];
+	      $bestScore = @seqAvgScore[$i];
+	    }
+	  }
+	}
+      }
+    }
+    if ($priority eq "score") {   # If the average score should be used to determine the most likely profile-type
+      ## Reset variables to unrealisticly high or low values
+      $best = 0;
+      $bestCount = 0;
+      $bestEval = 1000;
+      $bestScore = -1000;
+
+      for ($i = 0; $i < scalar(@seqDomCounts); $i++) {   # Go through all domain counts entries
+	if (@seqAvgScore[$i] > $bestScore) {   # If the current average score is larger than the previous top average score
+	  $best = $i;   # Set the best value to be the current value
+	  ## Set all other best variables to those corresponding to the current value
+	  $bestCount = @seqDomCounts[$i];
+	  $bestEval = @seqAvgE[$i];
+	  $bestScore = @seqAvgScore[$i];
+	}
+	if (@seqAvgScore[$i] == $bestScore) {   # If the current average score is equal to the previous top average score
+	  if (@seqAvgE[$i] < $bestEval) {   # If the current E-value is smaller than the previous top E-value
+	    $best = $i;   # Set the best value to be the current value
+	    ## Set all other best variables to those corresponding to the current value
+	    $bestCount = @seqDomCounts[$i];
+	    $bestEval = @seqAvgE[$i];
+	    $bestScore = @seqAvgScore[$i];
+	  }
+	  if (@seqAvgE[$i] == $bestEval) {   # If the current E-value is equal to the previous top E-value
+	    if (@seqDomCounts[$i] > $bestCount) {   # If the current number of domains is larger than the previous top number of domains
+	      $best = $i;   # Set the best value to be the current value
+	      ## Set all other best variables to those corresponding to the current value
+	      $bestCount = @seqDomCounts[$i];
+	      $bestEval = @seqAvgE[$i];
+	      $bestScore = @seqAvgScore[$i];
+	    }
+	  }
+	}
+      }
+    }
+    
+    if (@seqID[$best] ne "") {   # If the sequence ID of the most likely profile is not empty
+
+      $allanchorLens = @anchorLens[$best];  # Get the best anchor lengths
+      @allanchorLens = split(',', $allanchorLens);  # Split the anchor lens into an array
+
+      $chimeric = 0;   # Assume the sequence is not chimeric
+      if ($allow_reorder == 0) {   # If re-order of domain is not allowed
+	$domain_order = @allSeqDomains[$best];   # Gather the order the domains are found in
+	@domain_order = split(' ',$domain_order);   # Split the list into an array
+	@sorted_domain_order = sort {$a cmp $b} @domain_order;   # Sort the array alphabetically
+	for ($di = 0; $di <= scalar(@domain_order); $di++) {   # Go through the sorted array
+	  if ((@domain_order[$di] ne @sorted_domain_order[$di]) || (@problemCode[$best] =~ m/C/)) {    # Check if the order of the arrays differ at any poiny
+	    $chimeric = 1;   # If they do differ, mark the sequence as chimeric
+	  }
+	}
+      }
+
+      $seqDNALength = length(@seqDNA[$best]);   # Get the length of the DNA sequence
+
+      ## Print sequence and match data...
+      ## Order of columns in the output file:
+      ## ID Length Type Main/Compl Domains Avg.Eval Avg.Score Start End Start_domain End_domain Chimeric
+      if ($out_results == 1) {
+	print RESULTS @seqID[$best] . "\t" . $seqDNALength . "\t" . @seqITSType[$best] . "\t" . @seqCompl[$best] . "\t" . @seqDomCounts[$best] . "\t" . @seqAvgE[$best] . "\t" . @seqAvgScore[$best] . "\t" . @seqScoreSum[$best] . "\t" . @dnaStart[$best] . "\t" . @dnaEnd[$best] . "\t" . @startDomain[$best] . "\t" . @endDomain[$best] . "\t";   # Print sequence and match data to the results file
+	if ($chimeric == 1) {   # If the sequence was regarded chimeric
+	  print RESULTS "Chimeric\t";   # Add a chimeric tag to the entry
+	} else {   # If not chimeric
+	  print RESULTS "\t";   # Add an empty column
+	}
+	$allDomains = @allSeqDomains[$best];   # Get the domain order of the entry
+	$allDomains =~ tr/ /,/;   # Replace spaces with commas in the domain order string
+	$allDomains = substr($allDomains,0,length($allDomains) - 1);   # Remove the last character (a comma)
+	print RESULTS $allDomains;   # Write the domain order to the results file
+	print RESULTS "\t";   # Write a tab to the results file
+      }
+
+      if ($out_pos == 1) {
+	$out_all_pos = 1;
+	if ($out_all_pos == 1) {  # Output positions of all domains
+	  $seqPartLen = @dnaEnd[$best] - @dnaStart[$best] + 1;
+	  if ($seqPartLen < 0) {
+	    $seqPartLen = $seqPartLen * -1;
+	  }
+	  ## Print the positions of all identified domains to the position file
+	  print POS @seqID[$best] . "\t" . $seqDNALength . " bp." . "\t";
+	  if (@problemCode[$best] !~ m/S/) {
+	    print POS "SSU: " . @ssuStart[$best] . "-" . @ssuEnd[$best] . "\t";
+	  } else {
+	    print POS "SSU: Not found\t";
+	  }
+	  if ((@problemCode[$best] =~ m/X/) || ((@problemCode[$best] =~ m/[15]/) && (@problemCode[$best] =~ m/S/))) {
+	    print POS "ITS1: Not found\t";
+	  } else {
+	    if (@problemCode[$best] =~ m/O/) {
+	      print POS "ITS1: " . (@ssuEnd[$best] + 1) . "-" . (@ssuEnd[$best] + 1) . "\t";
+	    } else {
+	      print POS "ITS1: " . @its1Start[$best] . "-" . @its1End[$best] . "\t";
+	    }
+	  }
+	  if (@problemCode[$best] !~ m/[125OP]/) {
+	    print POS "5.8S: " . @midStart[$best] . "-" . @midEnd[$best] . "\t";
+	  } else {
+	    if (@problemCode[$best] =~ m/5/) {
+	      print POS "5.8S: Not found\t";
+	    } else {
+	      if (@problemCode[$best] =~ m/1/) {
+		print POS "5.8S: No start\t";
+	      } else {
+		if (@problemCode[$best] =~ m/O/) {
+		  print POS "5.8S: Overlap SSU\t";
+		} else {
+		  if (@problemCode[$best] =~ m/P/) {
+		    print POS "5.8S: Overlap LSU\t";
+		  } else {
+		    print POS "5.8S: No end\t";
+		  }
+		}
+	      }
+	    }
+	  }
+	  if ((@problemCode[$best] =~ m/Y/) ||((@problemCode[$best] =~ m/[25]/) && (@problemCode[$best] =~ m/L/))) {
+	    print POS "ITS2: Not found\t";
+	  } else {
+	    if (@problemCode[$best] =~ m/O/) {
+	      print POS "ITS2: " . (@lsuStart[$best] - 1) . "-" . (@lsuStart[$best] - 1) . "\t";
+	    } else {
+	      print POS "ITS2: " . @its2Start[$best] . "-" . @its2End[$best] . "\t";
+	    }
+	  }
+	  if (@problemCode[$best] !~ m/L/) {
+	    print POS "LSU: " . @lsuStart[$best] . "-" . @lsuEnd[$best] . "\t";
+	  } else {
+	    print POS "LSU: Not found" . "\t";
+	  }
+
+	  if (@problemCode[$best] =~ m/5/) {
+	    print POS "Broken or partial sequence, no 5.8S! ";
+	  } else {
+	    if (@problemCode[$best] =~ m/[12]/) {
+	      print POS "Broken or partial sequence, only partial 5.8S! ";
+	    }
+	  }
+	  if (@problemCode[$best] =~ m/B/) {
+	    print POS "ITS region too long! ";
+	  }
+	  if (@problemCode[$best] =~ m/O/) {
+	    print POS "5.8S seem to overlap with SSU! ";
+	  }
+	  if (@problemCode[$best] =~ m/P/) {
+	    print POS "5.8S seem to overlap with LSU! ";
+	  }
+	  if (@problemCode[$best] =~ m/C/) {
+	    print POS "Chimeric! ";
+	  }
+	  print POS "\n";
+
+	} else {  # Output only ITS positions
+	  $seqPartLen = @dnaEnd[$best] - @dnaStart[$best] + 1;
+	  print POS @seqID[$best] . "\t" . $seqDNALength . " bp." . "\t" . "ITS1: " . @its1Start[$best] . "-" . @its1End[$best] . "\t" . "ITS2: " . @its2Start[$best] . "-" . @its2End[$best] . "\n";  # Print the positions of the ITS sequences to the position file
+	}
+      }
+
+      if (@problem[$best] ne "") {
+	$foundProblem++;
+	print PROBLEM @seqID[$best] . "\t" . @problem[$best] . "\n";
+      }
+      
+      ## Set extended type string (the string going into the definition line of the FASTA file)
+      $extendedType = $profileIndex{@seqITSType[$best]} . " ITS sequence";
+      @itsCounts[ord(@seqITSType[$best])]++;   # Add one ITS to the appropriate counter
+
+      if (@seqCompl[$best] == 1) {   # If domains were found on complementary strand
+	$extendedStrand = "complementary strand";   # Set the strand string to complementary
+	$itsCompl++;   # Add one to the complementary strand counter
+      } else {   # If domains were found on the main strand
+	$extendedStrand = "main strand";   # Set the strand string to main
+	$itsMain++;   # Add one to the main strand counter
+      }
+      
+      ## Print (extracted) ITS sequence...
+      if ($out_fasta == 1) {   # If FASTA-output is on
+	if ($truncate == 0) {   # If the whole sequence should be kept in output file
+	  if ($chimeric == 0) {   # If the sequence is not chimeric
+	    if (@problemCode[$best] !~ m/[SL]/) {
+	      if ($out_preserve == 1) {  # If sequence headers should be preserved
+		print FASTA $headers{@seqID[$best]} . "\n";   # Write FASTA definition line
+	      } else {
+		print FASTA ">" . @seqID[$best] . "|" . @seqITSType[$best] ." " . $extendedType . " (" . $seqDNALength . " bp)  on " . $extendedStrand . "\n";   # Write FASTA definition line
+	      }
+	      print FASTA  @seqDNA[$best] . "\n";   # Write DNA sequence
+	    }
+	    if ($out_partial > 0) {
+	      if (@problemCode[$best] !~ m/[SL125]/) {
+	      if ($out_preserve == 1) {  # If sequence headers should be preserved
+		print FULLPARTIAL $headers{@seqID[$best]} . "\n";   # Write FASTA definition line
+	      } else {
+		print FULLPARTIAL ">" . @seqID[$best] . "|" . @seqITSType[$best] ." " . $extendedType . " (" . $seqDNALength . " bp) Full ITS region on " . $extendedStrand . "\n";   # Write FASTA defline
+	      }
+		print FULLPARTIAL  @seqDNA[$best] . "\n";   # Write DNA sequence
+	      } else {
+		$its1PartLen = @its1End[$best] - @its1Start[$best] + 1;
+		$its2PartLen = @its2End[$best] - @its2Start[$best] + 1;
+		if (($out_partial < $its1PartLen) && ($out_partial < $its2PartLen)) {
+		  if ($out_preserve == 1) {  # If sequence headers should be preserved
+		    print FULLPARTIAL $headers{@seqID[$best]} . "\n";   # Write FASTA definition line
+		  } else {
+		    print FULLPARTIAL ">" . @seqID[$best] . "|" . @seqITSType[$best] ." " . $extendedType . " (" . $seqDNALength . " bp) Partial ITS region on " . $extendedStrand . "\n";   # Write FASTA defline
+		  }
+		  print FULLPARTIAL  @seqDNA[$best] . "\n";   # Write DNA sequence
+		}
+	      }
+	    }
+	  } else {   # If sequence is regarded chimeric
+	    $itsChimeric++;   # Add one to the chimeric counter
+	    if ($out_preserve == 1) {  # If sequence headers should be preserved
+	      print CHIMERA $headers{@seqID[$best]} . "\n";   # Write FASTA definition line
+	    } else {
+	      print CHIMERA ">" . @seqID[$best] . "|" . @seqITSType[$best] ." Chimeric " . $extendedType . " (" . length(@seqDNA[$best]) . " bp)  on " . $extendedStrand . "\n";   # Write FASTA definition line
+	    }
+	    print CHIMERA  @seqDNA[$best] . "\n";   # Write DNA sequence
+	    $foundChimera++;
+	  }
+	} else {    # If only the ITS part of the sequence should be saved to output file
+	  $fastaStartPoint = @dnaStart[$best] - 1;   # Start extraction at the start of the first domain
+	  $fastaEndPoint = @dnaEnd[$best]+10;   # End extraction 10 bp after the last domain
+
+	  if (@seqDomCounts[$best] > 1) {   # If more than one domain was found
+	    if (substr(@startDomain[$best],0,5) eq "1_SSU")  {   # If the first domain was SSU
+	      $fastaStartPoint = @ssuEnd[$best];   # Set the start point of the extraction to the end of the SSU domain
+	    }
+	    if (substr(@startDomain[$best],0,5) eq "4_LSU")  {   # If the first domain was LSU
+	      $fastaStartPoint = @lsuEnd[$best];   # Set the start point of the extraction to the end of the LSU domain
+	    }
+	    if (substr(@endDomain[$best],0,5) eq "1_SSU")  {   # If the last domain was SSU
+	      $fastaEndPoint = @ssuStart[$best] - 1;   # Set the end point of the extraction to the start of the SSU domain
+	    }
+	    if (substr(@endDomain[$best],0,5) eq "4_LSU")  {   # If the last domain was LSU
+	      $fastaEndPoint = @lsuStart[$best] - 1;   # Set the end point of the extraction to the start of the LSU domain
+	    }
+	  }
+	  
+	  if ($fastaStartPoint < 0) {   # If the start point is smaller than zero, set the start point to zero
+	    $fastaStartPoint = 0;
+	  }
+	  if ($fastaEndPoint > length(@seqDNA[$best])) {   # If the end point is larger than the length of the sequence, set the end point to the sequence end
+	    $fastaEndPoint = length(@seqDNA[$best]);
+	  }
+
+	  $fastaLength = $fastaEndPoint - $fastaStartPoint + 1;
+
+	  if ($chimeric == 0) {   # If the sequence is not chimeric
+	    if (@problemCode[$best] !~ m/[SL]/) {
+	      if ($out_preserve == 1) {  # If sequence headers should be preserved
+		print FASTA $headers{@seqID[$best]} . "\n";   # Write FASTA definition line
+	      } else {
+		print FASTA ">" . @seqID[$best] . "|" . @seqITSType[$best] . " " . $extendedType . " (" . $fastaLength . " bp) on " . $extendedStrand . "\n";   # Write FASTA definition line
+	      }
+	      if ($anchorLen > 0) {
+		if ($fastaStartPoint - @allanchorLens[0] > 0) {
+		  print FASTA substr(@seqDNA[$best],$fastaStartPoint - @allanchorLens[0],$fastaLength + @allanchorLens[0] + @allanchorLens[3]) . "\n";   # Write DNA sequence
+		} else {
+		  print FASTA substr(@seqDNA[$best],0,$fastaLength + $fastaStartPoint + @allanchorLens[3]) . "\n";   # Write DNA sequence
+		}
+	      } else {
+		print FASTA substr(@seqDNA[$best],$fastaStartPoint,$fastaLength) . "\n";   # Write DNA sequence
+	      }
+	    }
+	    if ($out_partial > 0) {
+	      if (@problemCode[$best] !~ m/[SL125]/) {
+		if ($out_preserve == 1) {  # If sequence headers should be preserved
+		  print FULLPARTIAL $headers{@seqID[$best]} . "\n";   # Write FASTA definition line
+		} else {
+		  print FULLPARTIAL ">" . @seqID[$best] . "|" . @seqITSType[$best] ." " . $extendedType . " (" . $seqDNALength . " bp) Full ITS region on " . $extendedStrand . "\n";   # Write FASTA defline
+		}
+		print FULLPARTIAL substr(@seqDNA[$best],$fastaStartPoint,$fastaLength) . "\n";   # Write DNA sequence
+	      } else {
+		$its1PartLen = @its1End[$best] - @its1Start[$best] + 1;
+		$its2PartLen = @its2End[$best] - @its2Start[$best] + 1;
+		if (($out_partial < $its1PartLen) && ($out_partial < $its2PartLen)) {
+		  if ($out_preserve == 1) {  # If sequence headers should be preserved
+		    print FULLPARTIAL $headers{@seqID[$best]} . "\n";   # Write FASTA definition line
+		  } else {
+		    print FULLPARTIAL ">" . @seqID[$best] . "|" . @seqITSType[$best] ." " . $extendedType . " (" . $seqDNALength . " bp) Partial ITS region on " . $extendedStrand . "\n";   # Write FASTA defline
+		  }
+		  print FULLPARTIAL substr(@seqDNA[$best],$fastaStartPoint,$fastaLength) . "\n";   # Write DNA sequence
+		}
+	      }
+	    }
+	  } else {   # If sequence is regarded chimeric
+	    $itsChimeric++;   # Add one to the chimeric counter
+	    if ($out_preserve == 1) {  # If sequence headers should be preserved
+	      print CHIMERA $headers{@seqID[$best]} . "\n";   # Write FASTA definition line
+	    } else {
+	      print CHIMERA ">" . @seqID[$best] . "|" . @seqITSType[$best] . " Chimeric " . $extendedType . " (" . $fastaLength . " bp) on " . $extendedStrand . "\n";   # Write FASTA definition line
+	    }
+	    print CHIMERA @seqDNA[$best] . "\n";   # Write DNA sequence
+	    $foundChimera++;
+	  }
+	}
+      }
+
+      if ($out_joined == 1) {
+	if ($chimeric == 0) {   # If the sequence is not chimeric
+	  if ($out_preserve == 1) {  # If sequence headers should be preserved
+	    print JOINED $headers{@seqID[$best]} . "\n";   # Write FASTA definition line
+	  } else {
+	    print JOINED ">" . @seqID[$best] . "|" . @seqITSType[$best] . " " . $extendedType . " (" . $fastaLength . " bp) From domain " . @startDomain[$best] . " to " . @endDomain[$best] . " on " . $extendedStrand . " Found domains: ";   # Write FASTA definition line
+	    print JOINED substr(@allSeqDomains[$best],0,length(@allSeqDomains[$best]) - 1) . "\n";   # Write domain order
+	  }
+	  print JOINED substr(@seqDNA[$best],$fastaStartPoint,$fastaLength) . "\n";   # Write DNA sequence
+	} else {   # If sequence is regarded chimeric
+	  if ($out_preserve == 1) {  # If sequence headers should be preserved
+	    print JOINED $headers{@seqID[$best]} . "\n";   # Write FASTA definition line
+	  } else {
+	    print JOINED ">" . @seqID[$best] . "|" . @seqITSType[$best] . " Chimeric " . $extendedType . " (" . $fastaLength . " bp) From domain " . @startDomain[$best] . " to " . @endDomain[$best] . " on " . $extendedStrand . " Found domains: ";   # Write FASTA definition line
+	    print JOINED substr(@allSeqDomains[$best],0,length(@allSeqDomains[$best]) - 1) . "\n";   # Write domain order
+	  }
+	  print JOINED substr(@seqDNA[$best],$fastaStartPoint,$fastaLength) . "\n";   # Write DNA sequence
+	}
+      }
+
+      ## Write SSU sequence to file
+      if ($out_ssu == 1) {  # If SSU output is on
+	if (@problemCode[$best] !~ m/S/) {
+	  if ($only_full == 0) {
+	    $seqPartLen =  @ssuEnd[$best] - @ssuStart[$best] + 1;
+	    if ($out_preserve == 1) {  # If sequence headers should be preserved
+	      print SSU $headers{@seqID[$best]} . "\n";   # Write FASTA definition line
+	    } else {
+	      print SSU ">" . @seqID[$best] . "|" . @seqITSType[$best] . "|SSU " . "Extracted SSU sequence (" . $seqPartLen . " bp)\n";   # Write FASTA definition line
+	    }
+	    if (uc($anchor) ne "HMM") {
+	      if (@ssuStart[$best] - 1 - $anchorLen > 0) {
+		print SSU substr(@seqDNA[$best], @ssuStart[$best] - 1 - $anchorLen, $anchorLen);
+	      } else {
+		print SSU substr(@seqDNA[$best], 0, @ssuStart[$best] - 1);
+	      }
+	      print SSU substr(@seqDNA[$best], @ssuStart[$best] - 1, @ssuEnd[$best] - @ssuStart[$best] + 1);
+	      if (@ssuEnd[$best] - @ssuStart[$best] + 1 - $anchorLen > 0) {
+		print SSU substr(@seqDNA[$best], @ssuEnd[$best], $anchorLen);
+	      } else {
+		print SSU substr(@seqDNA[$best], @ssuEnd[$best]);
+	      }
+	    } else {
+	      print SSU substr(@seqDNA[$best], @ssuStart[$best] - 1, @ssuEnd[$best] - @ssuStart[$best] + 1);
+	    }
+	    print SSU "\n";   # Write DNA sequence
+	  }
+	}
+      }
+
+      ## Write LSU sequence to file
+      if ($out_lsu == 1) {  # If LSU output is on
+	if ($only_full == 0) {
+	  if (@problemCode[$best] !~ m/L/) {
+	    $seqPartLen = @lsuEnd[$best] - @lsuStart[$best] + 1;
+	    if ($out_preserve == 1) {  # If sequence headers should be preserved
+	      print LSU $headers{@seqID[$best]} . "\n";   # Write FASTA definition line
+	    } else {
+	      print LSU ">" . @seqID[$best] . "|" . @seqITSType[$best] . "|LSU " . "Extracted LSU sequence (" . $seqPartLen . " bp)\n";   # Write FASTA definition line
+	    }
+	    if (uc($anchor) ne "HMM") {
+	      if (@lsuStart[$best] - 1 - $anchorLen > 0) {
+		print LSU substr(@seqDNA[$best], @lsuStart[$best] - 1 - $anchorLen, $anchorLen);
+	      } else {
+		print LSU substr(@seqDNA[$best], 0, @lsuStart[$best] - 1);
+	      }
+	      print LSU substr(@seqDNA[$best], @lsuStart[$best] - 1, @lsuEnd[$best] - @lsuStart[$best] + 1);
+	      #if (@lsuEnd[$best] - @lsuStart[$best] + 1 - $anchorLen > 0) {
+	      #  print LSU substr(@seqDNA[$best], @lsuEnd[$best] - @lsuStart[$best] + 1, $anchorLen);
+	      #} else {
+	      #  print LSU substr(@seqDNA[$best], @lsuEnd[$best] - @lsuStart[$best] + 1);
+	      #}
+	    } else {
+	      print LSU substr(@seqDNA[$best], @lsuStart[$best] - 1, @lsuEnd[$best] - @lsuStart[$best] + 1);
+	    }
+	    print LSU "\n";   # Write DNA sequence
+	  }
+	}
+      }
+
+      ## Write 5.8S sequence to file
+      if ($out_58S == 1) {  # If 5.8S output is on
+	if (@problemCode[$best] !~ m/[125]/) {
+	  $seqPartLen = @midEnd[$best] - @midStart[$best] + 1;
+	  if ($out_preserve == 1) {  # If sequence headers should be preserved
+	    print MID $headers{@seqID[$best]} . "\n";   # Write FASTA definition line
+	  } else {
+	    print MID ">" . @seqID[$best] . "|" . @seqITSType[$best] . "|5.8S " . "Extracted 5.8S sequence (" . $seqPartLen . " bp)\n";   # Write FASTA definition line
+	  }
+	  if (uc($anchor) ne "HMM") {
+	    if (@midStart[$best] - 1 - $anchorLen > 0) {
+	      print MID substr(@seqDNA[$best], @midStart[$best] - 1 - $anchorLen, $anchorLen);
+	    } else {
+	      print MID substr(@seqDNA[$best], 0, @midStart[$best] - 1);
+	    }
+	    print MID substr(@seqDNA[$best], @midStart[$best] - 1, @midEnd[$best] - @midStart[$best] + 1);
+	    if (@midEnd[$best] - @midStart[$best] + 1 - $anchorLen > 0) {
+	      print MID substr(@seqDNA[$best], @midEnd[$best], $anchorLen);
+	    } else {
+	      print MID substr(@seqDNA[$best], @midEnd[$best]);
+	    }
+	  } else {
+	    print MID substr(@seqDNA[$best], @midStart[$best] - 1, @midEnd[$best] - @midStart[$best] + 1);
+	  }
+	  print MID "\n";   # Write DNA sequence
+	}
+      }
+
+      ## Write ITS1 sequence to file
+      if ($out_its1 == 1) {  # If ITS1 output is on
+	$seqPartLen = @its1End[$best] - @its1Start[$best] + 1;
+	if ($seqPartLen > 1) {
+	  if ( (($only_full == 0) && (@problemCode[$best] !~ m/[15]/)) || (($only_full == 1) && (@problemCode[$best] !~ m/[S15]/))) {
+	    if ($out_preserve == 1) {  # If sequence headers should be preserved
+	      print ITS1 $headers{@seqID[$best]} . "\n";   # Write FASTA definition line
+	    } else {
+	      print ITS1 ">" . @seqID[$best] . "|" . @seqITSType[$best] . "|ITS1 " . "Extracted ITS1 sequence " . @its1Start[$best] . "-" . @its1End[$best] . " (" . $seqPartLen . " bp)\n";   # Write FASTA definition line
+	    }
+	    if (uc($anchor) ne "HMM") {
+	      if (@its1Start[$best] - 1 - $anchorLen > 0) {
+		print ITS1 substr(@seqDNA[$best], @its1Start[$best] - 1 - $anchorLen, $anchorLen);
+	      } else {
+		print ITS1 substr(@seqDNA[$best], 0, @its1Start[$best] - 1);
+	      }
+	      print ITS1 substr(@seqDNA[$best], @its1Start[$best] - 1, @its1End[$best] - @its1Start[$best] + 1);
+	      if (length(@seqDNA[$best]) - @its1End[$best] - $anchorLen > 0) {
+		print ITS1 substr(@seqDNA[$best], @its1End[$best], $anchorLen);
+	      } else {
+		print ITS1 substr(@seqDNA[$best], @its1End[$best]);
+	      }
+	    } else {
+	      if (@its1Start[$best] - 1 - @allanchorLens[0] > 0) {
+		print ITS1 substr(@seqDNA[$best], @its1Start[$best] - 1 - @allanchorLens[0], @allanchorLens[0]);
+	      } else {
+		print ITS1 substr(@seqDNA[$best], 0, @its1Start[$best] - 1);
+	      }
+	      print ITS1 substr(@seqDNA[$best], @its1Start[$best] - 1, @its1End[$best] - @its1Start[$best] + 1);
+	      if (length(@seqDNA[$best]) - @its1End[$best] - @allanchorLens[1] > 0) {
+		print ITS1 substr(@seqDNA[$best], @its1End[$best], @allanchorLens[1]);
+	      } else {
+		print ITS1 substr(@seqDNA[$best], @its1End[$best]);
+	      }
+	      
+	    }
+	    print ITS1 "\n";   # Write DNA sequence
+	  }
+	}
+	if ($out_partial > 0) {
+	  if (@problemCode[$best] !~ m/[S15]/) {
+	    if ($out_preserve == 1) {  # If sequence headers should be preserved
+	      print ITS1PARTIAL $headers{@seqID[$best]} . "\n";   # Write FASTA definition line
+	    } else {
+	      print ITS1PARTIAL ">" . @seqID[$best] . "|" . @seqITSType[$best] . "|ITS1 " . "Extracted Full ITS1 sequence " . @its1Start[$best] . "-" . @its1End[$best] . " (" . $seqPartLen . " bp)\n";   # Write FASTA definition line
+	    }
+	    if (uc($anchor) ne "HMM") {
+	      if (@its1Start[$best] - 1 - $anchorLen > 0) {
+		print ITS1PARTIAL substr(@seqDNA[$best], @its1Start[$best] - 1 - $anchorLen, $anchorLen);
+	      } else {
+		print ITS1PARTIAL substr(@seqDNA[$best], 0, @its1Start[$best] - 1);
+	      }
+	      print ITS1PARTIAL substr(@seqDNA[$best], @its1Start[$best] - 1, @its1End[$best] - @its1Start[$best] + 1);
+	      if (length(@seqDNA[$best]) - @its1End[$best] - $anchorLen > 0) {
+		print ITS1PARTIAL substr(@seqDNA[$best], @its1End[$best], $anchorLen);
+	      } else {
+		print ITS1PARTIAL substr(@seqDNA[$best], @its1End[$best]);
+	      }
+	    } else {
+	      if (@its1Start[$best] - 1 - @allanchorLens[0] > 0) {
+		print ITS1PARTIAL substr(@seqDNA[$best], @its1Start[$best] - 1 - @allanchorLens[0], @allanchorLens[0]);
+	      } else {
+		print ITS1PARTIAL substr(@seqDNA[$best], 0, @its1Start[$best] - 1);
+	      }
+	      print ITS1PARTIAL substr(@seqDNA[$best], @its1Start[$best] - 1, @its1End[$best] - @its1Start[$best] + 1);
+	      if (length(@seqDNA[$best]) - @its1End[$best] - @allanchorLens[1] > 0) {
+		print ITS1PARTIAL substr(@seqDNA[$best], @its1End[$best], @allanchorLens[1]);
+	      } else {
+		print ITS1PARTIAL substr(@seqDNA[$best], @its1End[$best]);
+	      }
+	      
+	    }
+	    print ITS1PARTIAL "\n";   # Write DNA sequence
+	  } else {
+	    if ($out_partial < $seqPartLen) {
+	      if ($out_preserve == 1) {  # If sequence headers should be preserved
+		print ITS1PARTIAL $headers{@seqID[$best]} . "\n";   # Write FASTA definition line
+	      } else {
+		print ITS1PARTIAL ">" . @seqID[$best] . "|" . @seqITSType[$best] . "|ITS1 " . "Extracted Partial ITS1 sequence " . @its1Start[$best] . "-" . @its1End[$best] . " (" . $seqPartLen . " bp)\n";   # Write FASTA definition line
+	      }
+	      if (uc($anchor) ne "HMM") {
+		if (@its1Start[$best] - 1 - $anchorLen > 0) {
+		  print ITS1PARTIAL substr(@seqDNA[$best], @its1Start[$best] - 1 - $anchorLen, $anchorLen);
+		} else {
+		  print ITS1PARTIAL substr(@seqDNA[$best], 0, @its1Start[$best] - 1);
+		}
+		print ITS1PARTIAL substr(@seqDNA[$best], @its1Start[$best] - 1, @its1End[$best] - @its1Start[$best] + 1);
+		if (length(@seqDNA[$best]) - @its1End[$best] - $anchorLen > 0) {
+		  print ITS1PARTIAL substr(@seqDNA[$best], @its1End[$best], $anchorLen);
+		} else {
+		  print ITS1PARTIAL substr(@seqDNA[$best], @its1End[$best]);
+		}
+	      } else {
+		if (@its1Start[$best] - 1 - @allanchorLens[0] > 0) {
+		  print ITS1PARTIAL substr(@seqDNA[$best], @its1Start[$best] - 1 - @allanchorLens[0], @allanchorLens[0]);
+		} else {
+		  print ITS1PARTIAL substr(@seqDNA[$best], 0, @its1Start[$best] - 1);
+		}
+		print ITS1PARTIAL substr(@seqDNA[$best], @its1Start[$best] - 1, @its1End[$best] - @its1Start[$best] + 1);
+		if (length(@seqDNA[$best]) - @its1End[$best] - @allanchorLens[1] > 0) {
+		  print ITS1PARTIAL substr(@seqDNA[$best], @its1End[$best], @allanchorLens[1]);
+		} else {
+		  print ITS1PARTIAL substr(@seqDNA[$best], @its1End[$best]);
+		}
+	      }
+	      print ITS1PARTIAL "\n";   # Write DNA sequence
+	    }
+	  }
+	}
+      }
+      
+      ## Write ITS2 sequence to file
+      if ($out_its2 == 1) {  # If ITS2 output is on
+	$seqPartLen = @its2End[$best] - @its2Start[$best] + 1;
+	if ($seqPartLen > 1) {
+	  if ( (($only_full == 0) && (@problemCode[$best] !~ m/[25]/)) || (($only_full == 1) && (@problemCode[$best] !~ m/[L25]/))) {
+	    if ($out_preserve == 1) {  # If sequence headers should be preserved
+	      print ITS2 $headers{@seqID[$best]} . "\n";   # Write FASTA definition line
+	    } else {
+	      print ITS2 ">" . @seqID[$best] . "|" . @seqITSType[$best] . "|ITS2 " . "Extracted ITS2 sequence " . @its2Start[$best] . "-" . @its2End[$best] . " (" . $seqPartLen . " bp)\n";   # Write FASTA definition line
+	    }
+	    if (uc($anchor) ne "HMM") {
+	      if (@its2Start[$best] - 1 - $anchorLen > 0) {
+		print ITS2 substr(@seqDNA[$best], @its2Start[$best] - 1 - $anchorLen, $anchorLen);
+	      } else {
+		print ITS2 substr(@seqDNA[$best], 0, @its2Start[$best] - 1);
+	      }
+	      print ITS2 substr(@seqDNA[$best], @its2Start[$best] - 1, @its2End[$best] - @its2Start[$best] + 1);
+	      if (length(@seqDNA[$best]) - @its2End[$best] - $anchorLen > 0) {
+		print ITS2 substr(@seqDNA[$best], @its2End[$best], $anchorLen);
+	      } else {
+		print ITS2 substr(@seqDNA[$best], @its2End[$best]);
+	      }
+	    } else {
+	      if (@its2Start[$best] - 1 - @allanchorLens[2] > 0) {
+		print ITS2 substr(@seqDNA[$best], @its2Start[$best] - 1 - @allanchorLens[2], @allanchorLens[2]);
+	      } else {
+		print ITS2 substr(@seqDNA[$best], 0, @its2Start[$best] - 1);
+	      }
+	      print ITS2 substr(@seqDNA[$best], @its2Start[$best] - 1, @its2End[$best] - @its2Start[$best] + 1);
+	      if (length(@seqDNA[$best]) - @its2End[$best] - @allanchorLens[3] > 0) {
+		print ITS2 substr(@seqDNA[$best], @its2End[$best], @allanchorLens[3]);
+	      } else {
+		print ITS2 substr(@seqDNA[$best], @its2End[$best]);
+	      }
+	    }
+	    print ITS2 "\n";   # Write DNA sequence
+	  }
+	  if ($out_partial > 0) {
+	    if (@problemCode[$best] !~ m/[L25]/) {
+	      if ($out_preserve == 1) {  # If sequence headers should be preserved
+		print ITS2PARTIAL $headers{@seqID[$best]} . "\n";   # Write FASTA definition line
+	      } else {
+		print ITS2PARTIAL ">" . @seqID[$best] . "|" . @seqITSType[$best] . "|ITS2 " . "Extracted Full ITS2 sequence " . @its2Start[$best] . "-" . @its2End[$best] . " (" . $seqPartLen . " bp)\n";   # Write FASTA definition line
+	      }
+	      if (uc($anchor) ne "HMM") {
+		if (@its2Start[$best] - 1 - $anchorLen > 0) {
+		  print ITS2PARTIAL substr(@seqDNA[$best], @its2Start[$best] - 1 - $anchorLen, $anchorLen);
+		} else {
+		  print ITS2PARTIAL substr(@seqDNA[$best], 0, @its2Start[$best] - 1);
+		}
+		print ITS2PARTIAL substr(@seqDNA[$best], @its2Start[$best] - 1, @its2End[$best] - @its2Start[$best] + 1);
+		if (length(@seqDNA[$best]) - @its2End[$best] - $anchorLen > 0) {
+		  print ITS2PARTIAL substr(@seqDNA[$best], @its2End[$best], $anchorLen);
+		} else {
+		  print ITS2PARTIAL substr(@seqDNA[$best], @its2End[$best]);
+		}
+	      } else {
+		if (@its2Start[$best] - 1 - @allanchorLens[2] > 0) {
+		  print ITS2PARTIAL substr(@seqDNA[$best], @its2Start[$best] - 1 - @allanchorLens[2], @allanchorLens[2]);
+		} else {
+		  print ITS2PARTIAL substr(@seqDNA[$best], 0, @its2Start[$best] - 1);
+		}
+		print ITS2PARTIAL substr(@seqDNA[$best], @its2Start[$best] - 1, @its2End[$best] - @its2Start[$best] + 1);
+		if (length(@seqDNA[$best]) - @its2End[$best] - @allanchorLens[3] > 0) {
+		  print ITS2PARTIAL substr(@seqDNA[$best], @its2End[$best], @allanchorLens[3]);
+		} else {
+		  print ITS2PARTIAL substr(@seqDNA[$best], @its2End[$best]);
+		}
+	      }
+	      print ITS2PARTIAL "\n";   # Write DNA sequence
+	    } else {
+	      if ($out_partial < $seqPartLen) {
+		if ($out_preserve == 1) {  # If sequence headers should be preserved
+		  print ITS2PARTIAL $headers{@seqID[$best]} . "\n";   # Write FASTA definition line
+		} else {
+		  print ITS2PARTIAL ">" . @seqID[$best] . "|" . @seqITSType[$best] . "|ITS2 " . "Extracted Partial ITS2 sequence " . @its2Start[$best] . "-" . @its2End[$best] . " (" . $seqPartLen . " bp)\n";   # Write FASTA definition line
+		}
+		if (uc($anchor) ne "HMM") {
+		  if (@its2Start[$best] - 1 - $anchorLen > 0) {
+		    print ITS2PARTIAL substr(@seqDNA[$best], @its2Start[$best] - 1 - $anchorLen, $anchorLen);
+		  } else {
+		    print ITS2PARTIAL substr(@seqDNA[$best], 0, @its2Start[$best] - 1);
+		  }
+		  print ITS2PARTIAL substr(@seqDNA[$best], @its2Start[$best] - 1, @its2End[$best] - @its2Start[$best] + 1);
+		  if (length(@seqDNA[$best]) - @its2End[$best] - $anchorLen > 0) {
+		    print ITS2PARTIAL substr(@seqDNA[$best], @its2End[$best], $anchorLen);
+		  } else {
+		    print ITS2PARTIAL substr(@seqDNA[$best], @its2End[$best]);
+		  }
+		} else {
+		  if (@its2Start[$best] - 1 - @allanchorLens[2] > 0) {
+		    print ITS2PARTIAL substr(@seqDNA[$best], @its2Start[$best] - 1 - @allanchorLens[2], @allanchorLens[2]);
+		  } else {
+		    print ITS2PARTIAL substr(@seqDNA[$best], 0, @its2Start[$best] - 1);
+		  }
+		  print ITS2PARTIAL substr(@seqDNA[$best], @its2Start[$best] - 1, @its2End[$best] - @its2Start[$best] + 1);
+		  if (length(@seqDNA[$best]) - @its2End[$best] - @allanchorLens[3] > 0) {
+		    print ITS2PARTIAL substr(@seqDNA[$best], @its2End[$best], @allanchorLens[3]);
+		  } else {
+		    print ITS2PARTIAL substr(@seqDNA[$best], @its2End[$best]);
+		  }
+		}
+		print ITS2PARTIAL "\n";   # Write DNA sequence
+	      }
+	    }
+	  }
+	}
+      }
+
+      ## Output concatenated ITS1 + ITS2 sequences
+      if ($out_concat == 1) {
+	$seqPartLen1 = @its1End[$best] - @its1Start[$best] + 1;
+	$seqPartLen2 = @its2End[$best] - @its2Start[$best] + 1;
+	if (($seqPartLen1 >= $concat_minlen) && ($seqPartLen2 >= $concat_minlen)) {
+	  if ($out_preserve == 1) {  # If sequence headers should be preserved
+	    print CONCAT $headers{@seqID[$best]} . "\n";   # Write FASTA definition line
+	  } else {
+	    print CONCAT ">" . @seqID[$best] . "|" . @seqITSType[$best] . "|ITS1+2 " . "Concatenated ITS1 and ITS2 sequences (" . ($seqPartLen1 + $seqPartLen2) . " bp)\n";   # Write FASTA definition line
+	  }
+	  print CONCAT substr(@seqDNA[$best], @its1Start[$best] - 1, @its1End[$best] - @its1Start[$best] + 1);   # Write ITS1 DNA sequence
+	  print CONCAT "-----";  # Write spacer
+	  print CONCAT substr(@seqDNA[$best], @its2Start[$best] - 1, @its2End[$best] - @its2Start[$best] + 1);   # Write ITS2 DNA sequence
+	  print CONCAT "\n";  # Write newline
+	} else {
+	  if ($seqPartLen1 >= $concat_minlen) {
+	    if ($out_preserve == 1) {  # If sequence headers should be preserved
+	      print CONCAT $headers{@seqID[$best]} . "\n";   # Write FASTA definition line
+	    } else {
+	      print CONCAT ">" . @seqID[$best] . "|" . @seqITSType[$best] . "|ITS1 " . "ITS1 sequence (ITS2 too short) (" . $seqPartLen1 . " bp)\n";   # Write FASTA definition line
+	    }
+	    print CONCAT substr(@seqDNA[$best], @its1Start[$best] - 1, @its1End[$best] - @its1Start[$best] + 1);   # Write ITS1 DNA sequence
+	    print CONCAT "-----";  # Write spacer
+	    print CONCAT "\n";  # Write newline
+	  }
+	  if ($seqPartLen2 >= $concat_minlen) {
+	    if ($out_preserve == 1) {  # If sequence headers should be preserved
+	      print CONCAT $headers{@seqID[$best]} . "\n";   # Write FASTA definition line
+	    } else {
+	      print CONCAT ">" . @seqID[$best] . "|" . @seqITSType[$best] . "|ITS2 " . "ITS2 sequence (ITS1 too short) (" . $seqPartLen2 . " bp)\n";   # Write FASTA definition line
+	    }
+	    print CONCAT "-----";  # Write spacer
+	    print CONCAT substr(@seqDNA[$best], @its2Start[$best] - 1, @its2End[$best] - @its2Start[$best] + 1);   # Write ITS2 DNA sequence
+	    print CONCAT "\n";  # Write newline
+	  }
+	}
+      }
+
+      ## Print all sequences to the joined file for debugging
+      if ($out_joined == 1) {
+	$seqPartLen =  @ssuEnd[$best] - @ssuStart[$best] + 1;
+	if ($out_preserve == 1) {  # If sequence headers should be preserved
+	  print JOINED $headers{@seqID[$best]} . "\n";   # Write FASTA definition line
+	} else {
+	  print JOINED ">" . @seqID[$best] . "|" . @seqITSType[$best] . "|SSU " . "Extracted SSU sequence (" . $seqPartLen . " bp)\n";   # Write FASTA definition line
+	}
+	print JOINED substr(@seqDNA[$best], @ssuStart[$best] - 1, @ssuEnd[$best] - @ssuStart[$best] + 1) . "\n";   # Write DNA sequence
+	
+	$seqPartLen = @its1End[$best] - @its1Start[$best] + 1;
+	if ($out_preserve == 1) {  # If sequence headers should be preserved
+	  print JOINED $headers{@seqID[$best]} . "\n";   # Write FASTA definition line
+	} else {
+	  print JOINED ">" . @seqID[$best] . "|" . @seqITSType[$best] . "|ITS1 " . "Extracted ITS1 sequence " . @its1Start[$best] . "-" . @its1End[$best] . " (" . $seqPartLen . " bp)\n";   # Write FASTA definition line
+	}
+	print JOINED substr(@seqDNA[$best], @its1Start[$best] - 1, @its1End[$best] - @its1Start[$best] + 1) . "\n";   # Write DNA sequence
+	
+	$seqPartLen = @midEnd[$best] - @midStart[$best] + 1;
+	if ($out_preserve == 1) {  # If sequence headers should be preserved
+	  print JOINED $headers{@seqID[$best]} . "\n";   # Write FASTA definition line
+	} else {
+	  print JOINED ">" . @seqID[$best] . "|" . @seqITSType[$best] . "|5.8S " . "Extracted 5.8S sequence (" . $seqPartLen . " bp)\n";   # Write FASTA definition line
+	}
+	print JOINED substr(@seqDNA[$best], @midStart[$best] - 1, @midEnd[$best] - @midStart[$best] + 1) . "\n";   # Write DNA sequence
+	
+	$seqPartLen = @its2End[$best] - @its2Start[$best] + 1;
+	if ($out_preserve == 1) {  # If sequence headers should be preserved
+	  print JOINED $headers{@seqID[$best]} . "\n";   # Write FASTA definition line
+	} else {
+	  print JOINED ">" . @seqID[$best] . "|" . @seqITSType[$best] . "|ITS2 " . "Extracted ITS2 sequence (" . $seqPartLen . " bp)\n";   # Write FASTA definition line
+	}
+	print JOINED substr(@seqDNA[$best], @its2Start[$best] - 1, @its2End[$best] - @its2Start[$best] + 1) . "\n";   # Write DNA sequence
+	
+	$seqPartLen = @lsuEnd[$best] - @lsuStart[$best] + 1;
+	if ($out_preserve == 1) {  # If sequence headers should be preserved
+	  print JOINED $headers{@seqID[$best]} . "\n";   # Write FASTA definition line
+	} else {
+	  print JOINED ">" . @seqID[$best] . "|" . @seqITSType[$best] . "|LSU " . "Extracted LSU sequence (" . $seqPartLen . " bp)\n";   # Write FASTA definition line
+	}
+	print JOINED substr(@seqDNA[$best], @lsuStart[$best] - 1, @lsuEnd[$best] - @lsuStart[$best] + 1) . "\n";   # Write DNA sequence
+      }
+      
+      
+      if ($out_results == 1) {
+	## Print info on all matches, also not top ones to the results file...
+	for ($i = 0; $i < scalar(@seqITSType); $i++) {   # Go through all the possible ITS types...
+	  print RESULTS @seqITSType[$i] . ": " . @seqDomCounts[$i] . " " . @seqAvgE[$i] . " " . @seqAvgScore[$i];   # Write some info on this type (Type, Domain count, Average E-value, Average score)
+	  if ($i < scalar(@seqITSType) - 1) {   # If this is no the last domain type
+	    print RESULTS ", ";   # Write a comma
+	  }
+	}
+	print RESULTS "\n";   # Write end of line
+      }
+    }
+    
+    ## Undefine all used arrays for the next round...
+    undef @seqID;
+    undef @seqITSType;
+    undef @seqCompl;
+    undef @seqDomCounts;
+    undef @seqAvgE;
+    undef @seqAvgScore;
+    undef @seqScoreSum;
+    undef @seqDNA;
+    undef @allSeqDomains;
+    undef @dnaStart;
+    undef @dnaEnd;
+    undef @ssuStart;
+    undef @ssuEnd;
+    undef @lsuStart;
+    undef @lsuEnd;
+    undef @midStart;
+    undef @midEnd;
+    undef @its1Start;
+    undef @its1End;
+    undef @its2Start;
+    undef @its2End;
+    undef @startDomain;
+    undef @endDomain;
+    undef @domain_order;
+    undef @sorted_domain_order;
+    undef @problem;
+    undef @problemCode;
+    undef @anchorLens;
+  }
+  if (@item[0] ne "") {   # Add this entry to the set (regardless if the entry has the same ID as entries already in the set), as long as it is non-empty
+    push(@seqID, @item[0]);   # Add sequence ID
+    push(@seqITSType, @item[1]);   # Add ITS type
+    push(@seqCompl, @item[2]);   # Add main/complementary strand info
+    push(@seqDomCounts, @item[3]);   # Add domain count
+    push(@seqAvgE, @item[4]);   # Add average E-value
+    push(@seqAvgScore, @item[5]);   # Add average score
+    push(@seqScoreSum, @item[6]);   # Add sum-of-scores
+    push(@seqDNA, @item[7]);   # Add DNA sequence
+    
+    ## Determine first and last domains, and their positions
+    ## Set variables to unrealistic values
+    $dnaEnd = length(@item[7]);
+    $startDomain = "***";
+    $endDomain = "***";
+    $dnaStart = 1;
+    $allDomains = "";
+    $ssuStart = 1;
+    $ssuEnd = 0;
+    $midStart = 1;
+    $midEnd = 0;
+    $lsuStart = 1;
+    $lsuEnd = 0;
+    $its1Start = 1;
+    $its1End = length(@item[7]);
+    $its2Start = 1;
+    $its2End = length(@item[7]);
+    $problem = "";
+    $problemCode = "";
+
+    $ssuFound = 0;
+    $lsuFound = 0;
+    $midSFound = 0;
+    $midEFound = 0;
+    $order = "";
+    undef @hitanchorlens;
+    
+    for ($i = 8; $i < scalar(@item); $i++) {   # Go through the list of found domains in this sequence
+      ($hitStart,$hitEnd,$hitProfile,$hitScore,$hitEval,$hitanchorlen) = split(';',@item[$i]);   # Separate the hit stats into variables
+      $allDomains = $allDomains . $hitProfile . " ";   # Add found domain to the list of all domains
+
+#      if ($hitStart < $dnaStart) {   # If this domain is the first one so far
+#	$dnaStart = $hitStart;   # Set the start of the ITS sequence to this domain's start
+#	$startDomain = $hitProfile;   # Set this domain as the starting domain
+#      }
+#      if ($hitEnd > $dnaEnd) {   # If this domain is the last one so far
+#	$dnaEnd = $hitEnd;   # Set the end of the ITS sequence to this domain's end
+#	$endDomain = $hitProfile;   # Set this domain as the ending domain
+#      }
+
+      if (substr($hitProfile,0,5) eq "1_SSU" ) {   # If this domain is the SSU's end
+	$dnaStart = $hitStart;   # Set the start of the ITS sequence to this domain's start
+	$startDomain = $hitProfile;   # Set this domain as the starting domain
+	$ssuStart = 1;   # Set the start of the SSU sequence to the start of the sequence
+	$ssuEnd = $hitEnd;   # Set the end of the SSU sequence to this domain's end
+	$its1Start = $hitEnd + 1;   # Set the start of the ITS1 sequence to right after this domain's end
+	if ($midEFound == 0) {
+	  $its2Start = $hitEnd + 1;   # Set the end of the ITS2 sequence to right after this domain's end
+	}
+	$ssuFound = 1;
+	$order = $order . "1";
+	@hitanchorlens[1] = $hitanchorlen;
+      }
+
+      if (substr($hitProfile,0,5) eq "2_5.8" ) {   # If this domain is the 5.8S's start
+	if ($startDomain eq "***") {
+	  $dnaStart = 1;   # Set the start of the ITS sequence to the start of the sequence
+	  $startDomain = $hitProfile;   # Set this domain as the starting domain
+	}
+	$its1End = $hitStart - 1;   # Set the end of the ITS1 sequence to right before this domain's start
+	if ($midEFound == 0) {
+	  $its2Start = $hitEnd + 1;   # Set the start of the ITS2 sequence to right after this domain's end
+	}
+	$midStart = $hitStart;   # Set the start of the 5.8S sequence to this domain's start
+	$midSFound = 1;
+	$order = $order . "2";
+	@hitanchorlens[2] = $hitanchorlen;
+      }
+
+      if (substr($hitProfile,0,5) eq "3_End" ) {   # If this domain is the 5.8S's end
+	if ($startDomain eq "***") {
+	  $dnaStart = 1;   # Set the start of the ITS sequence to the start of the sequence
+	  $startDomain = $hitProfile;   # Set this domain as the starting domain
+	}
+	if ($midSFound == 0) {
+	  $its1End = $hitStart - 1;   # Set the end of the ITS1 sequence to right before this domain's start
+	}
+	$its2Start = $hitEnd + 1;   # Set the end of the ITS2 sequence to right after this domain's end
+	$midEnd = $hitEnd;   # Set the end of the 5.8S sequence to this domain's end
+	$midEFound = 1;
+	$order = $order . "3";
+	@hitanchorlens[3] = $hitanchorlen;
+      }
+
+      if (substr($hitProfile,0,5) eq "4_LSU" ) {   # If this domain is the LSU's start
+	if ($startDomain eq "***") {
+	  $dnaStart = 1;   # Set the start of the ITS sequence to the start of the sequence
+	  $startDomain = $hitProfile;   # Set this domain as the starting domain
+	}
+	$dnaEnd = $hitEnd;   # Set the end of the ITS sequence to this domain's end
+	$endDomain = $hitProfile;   # Set this domain as the ending domain
+	$lsuStart = $hitStart;   # Set the start of the LSU sequence to the start of the this domain
+	$lsuEnd = length(@item[7]);   # Set the end of the LSU sequence to the end of the sequence
+	$its2End = $hitStart - 1;   # Set the end of the ITS2 sequence to right before this domain's start
+	if ($midSFound == 0) {
+	  $its1End = $hitStart - 1;   # Set the end of the ITS1 sequence to right before this domain's start
+	}
+	$lsuFound = 1;
+	$order = $order . "4";
+	@hitanchorlens[4] = $hitanchorlen;
+      }
+    }
+
+    if ($ssuFound == 0) {
+      if ($problem ne "") {
+	$problem = $problem . "; ";
+      }
+      $problem = $problem . "End of SSU sequence not found";
+      $problemCode = $problemCode . "S";
+      $order = "1" .  $order;
+    }
+
+    if ($lsuFound == 0) {
+      if ($problem ne "") {
+	$problem = $problem . "; ";
+      }
+      $problem = $problem . "Start of LSU sequence not found";
+      $problemCode = $problemCode . "L";
+      $order = $order . "4";
+    }
+
+    if (($midSFound == 0) && ($midEFound == 0)) {
+      if ($problem ne "") {
+	$problem = $problem . "; ";
+      }
+      $problem = $problem . "The 5.8S sequence was not found at all";
+      $problemCode = $problemCode . "5";
+    } else {
+      if ($midSFound == 0) {
+	if ($problem ne "") {
+	  $problem = $problem . "; ";
+	}
+	$problem = $problem . "Start of 5.8S sequence not found";
+	$problemCode = $problemCode . "1";
+      }
+      if ($midEFound == 0) {
+	if ($problem ne "") {
+	  $problem = $problem . "; ";
+	}
+	$problem = $problem . "End of 5.8S sequence not found";
+	$problemCode = $problemCode . "2";
+      }
+    }
+
+    if ($dnaEnd - $dnaStart > 1500) {
+      if ($problem ne "") {
+	$problem = $problem . "; ";
+      }
+      $problem = $problem . "ITS region is suspiciously long (> 1500 bp)";
+      $problemCode = $problemCode . "B";
+    }
+
+    if (length($order) == 4) {
+      if ($order ne "1234") {
+	if ($problem ne "") {
+	  $problem = $problem . "; ";
+	}
+	$problem = $problem . "Domains found in wrong order, sequence may be chimeric";
+	$problemCode = $problemCode . "C";
+      }
+    } else {
+      if ((substr($order,0,1) ne "1") || (substr($order,-1) ne "4")) {
+	if ($problem ne "") {
+	  $problem = $problem . "; ";
+	}
+	$problem = $problem . "Domains found in wrong order, sequence may be chimeric";
+	$problemCode = $problemCode . "C";
+      } else {
+	if (($midSFound == 1) && ($midEFound == 1)) {
+	  if ($order !~ m/23/) {
+	    if ($problem ne "") {
+	      $problem = $problem . "; ";
+	    }
+	    $problem = $problem . "Domains found in wrong order, sequence may be chimeric";
+	    $problemCode = $problemCode . "C";
+	  }
+	}
+      }
+    }
+
+    if ($its1End - $its1Start < 0) {  # 5.8S overlaps SSU
+      if (($midSFound == 1) && ($ssuFound == 1)) {
+	$problemCode = $problemCode . "OC";
+	$its1Start = 0;
+	$its1End = 0;
+	if ($problem ne "") {
+	  $problem = $problem . "; ";
+	}
+	$problem = $problem . "SSU seems to overlap 5.8S, sequence may be chimeric";
+      } else {
+	$problemCode = $problemCode . "X";
+	$its1Start = 0;
+	$its1End = 0;
+	if ($problem ne "") {
+	  $problem = $problem . "; ";
+	}
+	$problem = $problem . "No ITS1 sequence";
+      }
+    }
+    if ($its2End - $its2Start < 0) {  # 5.8S overlaps LSU
+      if (($midEFound == 1) && ($lsuFound == 1)) {
+	$problemCode = $problemCode . "PC";
+	$its2Start = 0;
+	$its2End = 0;
+	if ($problem ne "") {
+	  $problem = $problem . "; ";
+	}
+	$problem = $problem . "LSU seems to overlap 5.8S, sequence may be chimeric";
+       } else {
+	$problemCode = $problemCode . "Y";
+	$its2Start = 0;
+	$its2End = 0;
+	if ($problem ne "") {
+	  $problem = $problem . "; ";
+	}
+	$problem = $problem . "No ITS2 sequence";
+      }
+    }
+
+
+    if (($problemCode =~ m/[15]/) && ($problemCode =~ m/[S]/)) {  # Sequence lack all indications of ITS1
+      $its1Start = 0;
+      $its1End = 0;
+    }
+    if (($problemCode =~ m/[25]/) && ($problemCode =~ m/[L]/)) {  # Sequence lack all indications of ITS2
+      $its2Start = 0;
+      $its2End = 0;
+    }
+
+
+    $allhitanchorlens = @hitanchorlens[1] . "," . @hitanchorlens[2] . "," . @hitanchorlens[3] . "," . @hitanchorlens[4];
+
+    push(@allSeqDomains,$allDomains);   # Add list of all domains
+    push(@dnaStart, $dnaStart);   # Add start of the ITS sequence
+    push(@startDomain, $startDomain);   # Add start domain
+    push(@dnaEnd, $dnaEnd);   # Add end of ITS sequence
+    push(@endDomain, $endDomain);   # Add end domain
+    push(@ssuStart, $ssuStart);   # Add start of the SSU sequence
+    push(@ssuEnd, $ssuEnd);   # Add end of the SSU sequence
+    push(@midStart, $midStart);   # Add start of the 5.8S sequence
+    push(@midEnd, $midEnd);   # Add end of the 5.8S sequence
+    push(@lsuStart, $lsuStart);   # Add start of the LSU sequence
+    push(@lsuEnd, $lsuEnd);   # Add end of the LSU sequence
+    push(@its1Start, $its1Start);   # Add start of the ITS1 sequence
+    push(@its1End, $its1End);   # Add end of the ITS1 sequence
+    push(@its2Start, $its2Start);   # Add start of the ITS2 sequence
+    push(@its2End, $its2End);   # Add end of the ITS2 sequence
+    push(@problem, $problem);  # Add potential problem info
+    push(@problemCode, $problemCode);  # Add potential problem info in code form
+    push(@anchorLens, $allhitanchorlens);  # Add the anchor lengths
+  }
+  $lc++;   # Increase the line count by one
+}
+
+## Save results to the summary file
+if ($out_sum == 1) {   # If summary file should be written
+  $itsTotal = 0;   # Reset the total ITS sum
+  foreach $typeCount (@itsCounts) {   # Add ITSs from all different origins
+    $itsTotal += $typeCount;   # Add the number of ITSs for this origin
+  }
+  ## Write info on the found ITS sequences to the summary file
+  print SUMMARY "Sequences detected as ITS by ITSx:\t$itsTotal\n";
+  print SUMMARY "  On main strand:         \t$itsMain\n";
+  print SUMMARY "  On complementary strand:\t$itsCompl\n";
+  if ($allow_reorder == 0) {   # If re-ordering of domains is not allowed
+    print SUMMARY "Sequences detected as chimeric by ITSx:\t$itsChimeric\n";   # Write the number of reported chimeric sequences to the summary file
+  }
+  ## Write info on the found ITS sequence types to the summary file
+  print SUMMARY "ITS sequences by preliminary origin:\n";
+  print SUMMARY "  Alveolates:      \t" . int(@itsCounts[ord("A")]) . "\n";
+  print SUMMARY "  Amoebozoa:       \t" . int(@itsCounts[ord("D")]) . "\n";
+  print SUMMARY "  Bacillariophyta: \t" . int(@itsCounts[ord("C")]) . "\n";
+  print SUMMARY "  Brown algae:     \t" . int(@itsCounts[ord("I")]) . "\n";
+  print SUMMARY "  Bryophytes:      \t" . int(@itsCounts[ord("B")]) . "\n";
+  print SUMMARY "  Euglenozoa:      \t" . int(@itsCounts[ord("E")]) . "\n";
+  print SUMMARY "  Eustigmatophytes:\t" . int(@itsCounts[ord("U")]) . "\n";
+  print SUMMARY "  Fungi:           \t" . int(@itsCounts[ord("F")]) . "\n";
+  print SUMMARY "  Green algae:     \t" . int(@itsCounts[ord("G")]) . "\n";
+  print SUMMARY "  Liverworts:      \t" . int(@itsCounts[ord("L")]) . "\n";
+  print SUMMARY "  Metazoa:         \t" . int(@itsCounts[ord("M")]) . "\n";
+  print SUMMARY "  Microsporidia:   \t" . int(@itsCounts[ord("N")]) . "\n";
+  print SUMMARY "  Oomycetes:       \t" . int(@itsCounts[ord("O")]) . "\n";
+  print SUMMARY "  Prymnesiophytes: \t" . int(@itsCounts[ord("P")]) . "\n";
+  print SUMMARY "  Raphidophytes:   \t" . int(@itsCounts[ord("Q")]) . "\n";
+  print SUMMARY "  Red algae:       \t" . int(@itsCounts[ord("H")]) . "\n";
+  print SUMMARY "  Rhizaria:        \t" . int(@itsCounts[ord("R")]) . "\n";
+  print SUMMARY "  Synurophyceae:   \t" . int(@itsCounts[ord("S")]) . "\n";
+  print SUMMARY "  Tracheophyta:    \t" . int(@itsCounts[ord("T")]) . "\n";
+  print SUMMARY "-----------------------------------------------------------------\n";
+  close (SUMMARY);   # Close the summary file
+}
+
+if ($out_results == 1) {
+  close (RESULTS);   # Close the results file
+}
+close (RAWOUT);   # Close the raw output file
+close (PROBLEM);   # Close the file for problematic entries
+if ($foundProblem == 0) {
+  `rm $output.problematic.txt 2> /dev/null`;
+}
+if ($out_pos == 1) {
+  close (POS);   # Close the positions file
+}
+if ($out_fasta == 1) {   # If FASTA output is on
+  close (FASTA);   # Close the FASTA output file
+  if ($allow_reorder == 0) {   # If re-ordering of domains is not allowed
+    close (CHIMERA);   # Close the chimera file
+    if ($foundChimera == 0) {
+      `rm $output.chimeric.fasta 2> /dev/null`;
+    }
+  }
+  if ($out_partial > 0) {
+    close FULLPARTIAL;
+  }
+}
+if ($out_ssu == 1) {
+  close (SSU);   # Close the SSU file
+}
+if ($out_lsu == 1) {
+  close (LSU);   # Close the LSU file
+}
+if ($out_58S == 1) {
+  close (MID);   # Close the 58S file
+}
+if ($out_its1 == 1) {
+  close (ITS1);   # Close the ITS1 file
+  if ($out_partial > 0) {
+    close ITS1PARTIAL;
+  }
+}
+if ($out_its2 == 1) {
+  close (ITS2);   # Close the ITS2 file
+  if ($out_partial > 0) {
+    close ITS2PARTIAL;
+  }
+}
+if ($out_concat == 1) {
+  close CONCAT;
+}
+
+if ($out_not == 1) {   # If not-found output is on
+  open (NOTFOUND, "$output\_no_detections.txt");   # Open the not-found output file
+  while ($line = <NOTFOUND>) {  # Read all entries from file
+    chomp($line);  # Remove newline char
+    push(@nodetectionlist,$line);  # Add to non-detection list
+  }
+  close NOTFOUND;
+
+  open (NOTFOUND, ">$output\_no_detections.fasta");   # Create a not-found FASTA output file
+  foreach $seqID (@nodetectionlist) {  # For all non-detections
+    $seq = $sequenceDB{"$seqID"};  # Get sequence from sequence database
+    print NOTFOUND ">$seqID\n";  # Print not found sequence ID
+    print NOTFOUND $seq . "\n";  # Print not found sequence
+  }
+  close NOTFOUND;
+}
+
+## Clean up and finish
+
+if ($pipeline == 0) {   # If ITSx is not called from the pipeline mode (i.e. from ITSx)
+  if ($save_raw == 1) {   # If raw data should be saved
+    `mv $tempDir $output\_ITSx_raw_output`;   # Change the name of the temporary directory to ..._ITSx_raw_output
+  } else {   # Else, discard the raw data
+    `rm -rf $tempDir`;   # Remove the temporary directory
+  }
+}
+
+## Get the current time and output a finished message
+$now = localtime;
+if ($silent == 0) {
+  print STDERR "$now : Extraction finished!\n";
+  print STDERR "-----------------------------------------------------------------\n";
+  print STDERR "Thank you for using ITSx!\n";
+  print STDERR "Please report bugs or unsupported lineages to itsx\@microbiology.se\n";
+  print STDERR "\n";
+}
+
+## Write end time a summary file
+if ($pipeline == 0) {    # If not running in pipeline mode
+  if ($out_sum == 1) {   # If summary output is on
+    open (SUMMARY, ">>$output.summary.txt");   # Append to the summary file
+    print SUMMARY "ITSx run finished at $now.\n";   # Write ending time for the analysis
+    close (SUMMARY);   # Close summary file
+  }
+}
+
+sub hmmerSearch {
+  $hmmerCommand = $_[0];
+  $outputFile = $_[1];
+  $strand = $_[2];
+  $profileSet = $_[3];
+  open (HMMEROUTPUT, ">$outputFile");
+  open (HMMER, "$hmmerCommand |");
+  $totalHitCount = 0;
+  $hitCount = 0;
+  $SSUCount = 0;
+  $LSUCount = 0;
+  $startCount = 0;
+  $endCount = 0;
+  $maxCount = 1;
+  while (chomp($line = <HMMER>)) {
+    if (substr($line,0,6) eq "Query:") {
+      $hitCount = 0;
+      $SSUCount = 0;
+      $LSUCount = 0;
+      $startCount = 0;
+      $endCount = 0;
+      undef @bestScore;
+      undef @bestEntry;
+      
+      $query = substr($line,7);
+      $queryLength = $query;
+      $queryLength =~ s/.*  *//;
+      $queryLength =~ s/[^0-9]//g;
+      $query =~ s/  *//;
+      $query =~ s/  *.*//;
+      
+      print HMMEROUTPUT "## New query:\t$query\t$queryLength\n";
+
+    }
+    if (substr($line,0,12) eq "Description:") {
+      $desc = $line;
+    }
+    if (substr($line,0,3) eq ">> ") {
+      ($tempshit,$hmmerSubjectName) = split(' ',$line);
+    }
+    if ($line =~ m/[0-9] ! /) {
+      $stats = $line;
+      $stats =~ s/  */\t/g;
+      ($empty,$no,$excl,$score) = split('\t',$stats);
+      $hitCount++;
+      $totalHitCount++;
+      if ($maxCount == 0) {
+	print HMMEROUTPUT "$query\t$hmmerSubjectName\t$queryLength$stats\n";
+      } else {
+	if (substr($hmmerSubjectName,0,5) eq "1_SSU") {
+	  $SSUCount++;
+	  if ($SSUCount <= $maxCount) {
+	    @bestScore[1] = $score;
+	    @bestEntry[1] = "$query\t$hmmerSubjectName\t$queryLength$stats\n";
+	  } else {
+	    if ($score > @bestScore[1]) {
+	      @bestScore[1] = $score;
+	      @bestEntry[1] = "$query\t$hmmerSubjectName\t$queryLength$stats\n";
+	    }
+	  }
+	}
+	if (substr($hmmerSubjectName,0,5) eq "4_LSU") {
+	  $LSUCount++;
+	  if ($LSUCount <= $maxCount) {
+	    @bestScore[4] = $score;
+	    @bestEntry[4] = "$query\t$hmmerSubjectName\t$queryLength$stats\n";
+	  } else { 
+	    if ($score > @bestScore[4]) {
+	      @bestScore[4] = $score;
+	      @bestEntry[4] = "$query\t$hmmerSubjectName\t$queryLength$stats\n";
+	    }
+	  }
+	}
+	if (substr($hmmerSubjectName,0,5) eq "2_5.8") {
+	  $startCount++;
+	  if ($startCount <= $maxCount) {
+	    @bestScore[2] = $score;
+	    @bestEntry[2] = "$query\t$hmmerSubjectName\t$queryLength$stats\n";
+	  } else {
+	    if ($score > @bestScore[2]) {
+	      @bestScore[2] = $score;
+	      @bestEntry[2] = "$query\t$hmmerSubjectName\t$queryLength$stats\n";
+	    }
+	  }
+	}
+	if (substr($hmmerSubjectName,0,5) eq "3_End") {
+	  $endCount++;
+	  if ($endCount <= $maxCount) {
+	    @bestScore[3] = $score;
+	    @bestEntry[3] = "$query\t$hmmerSubjectName\t$queryLength$stats\n";
+	  } else {
+	    if ($score > @bestScore[3]) {
+	      @bestScore[3] = $score;
+	      @bestEntry[3] = "$query\t$hmmerSubjectName\t$queryLength$stats\n"
+	    }
+	  }
+	}
+      }
+    }
+    if (substr($line,0,2) eq "//") {
+      if ($maxCount > 0) {
+	if ($SSUCount > 0) {
+	  print HMMEROUTPUT @bestEntry[1];
+	}
+	if ($startCount > 0) {
+	  print HMMEROUTPUT @bestEntry[2];
+	}
+	if ($endCount > 0) {
+	  print HMMEROUTPUT @bestEntry[3];
+	}
+	if ($LSUCount > 0) {
+	  print HMMEROUTPUT @bestEntry[4];
+	}
+      }
+      print HMMEROUTPUT "//\n";
+    }
+  }
+  close(HMMEROUTPUT);
+  $now = localtime;
+  if ($silent == 0) {
+    if ($strand eq "M") {
+      print STDERR "    $now : " . ucfirst($profileIndex{$profileSet}) . " analysis of main strand finished.\n";   # Print finished type
+    } else {
+      print STDERR "    $now : " . ucfirst($profileIndex{$profileSet}) . " analysis of complementary strand finished.\n";   # Print finished type
+    }
+  }
+}
+
+
+## Please send beers, pizzas, cakes, fruit pies, job positions and other types of feedback to:
+## johan.bengtsson [at] microbiology.se
+## Looking forward to hearing from you.... visit my website: www.microbiology.se for info on my research
+## //Johan Bengtsson, 2012-2014