itsx: ITSx comparison

comparison ITSx @ 0:f82c70f54bd7 draft

Uploaded

author	okorol
date	Tue, 24 Mar 2015 12:02:48 -0400
parents
children

comparison

equal deleted inserted replaced

--1:000000000000
+:f82c70f54bd7
+#!/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

Mercurial > repos > okorol > itsx

comparison ITSx @ 0:f82c70f54bd7 draft