# HG changeset patch
# User aafc-mbb
# Date 1457986508 14400
# Node ID 57e76d405a976e34621597135c041e6d46d2dd28
Uploaded
diff -r 000000000000 -r 57e76d405a97 ITSx
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/ITSx Mon Mar 14 16:15:08 2016 -0400
@@ -0,0 +1,2695 @@
+#!/usr/bin/perl
+# ITSx ITS Extractor
+$app_title = "ITSx -- Identifies ITS sequences and extracts the ITS region";
+$app_author = "Johan Bengtsson-Palme et al., University of Gothenburg";
+$app_version = "1.0.11";
+$app_message = "";
+# ----------------------------------------------------------------- #
+
+# License information
+$license =
+" ITSx - ITS Extractor -- Identifies ITS sequences and extracts the ITS region\
+ Copyright (C) 2012-2014 Johan Bengtsson-Palme et al.\
+\
+ This program is free software: you can redistribute it and/or modify\
+ it under the terms of the GNU General Public License as published by\
+ the Free Software Foundation, either version 3 of the License, or\
+ (at your option) any later version.\
+\
+ This program is distributed in the hope that it will be useful,\
+ but WITHOUT ANY WARRANTY; without even the implied warranty of\
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the\
+ GNU General Public License for more details.\
+\
+ You should have received a copy of the GNU General Public License\
+ along with this program, in a file called 'license.txt'\
+ If not, see: http://www.gnu.org/licenses/.\
+";
+
+## BUGS:
+$bugs = "New features in this version ($app_version):\
+- None\
+\
+Fixed bugs in this version ($app_version):\
+- Fixed a bug causing newline characters to be occasionally skipped in the 'its1.full_and_partial.fasta' FASTA output file when the '--anchor' option was used\
+\
+Known bugs in this version ($app_version):\
+- None\
+";
+
+## OPTIONS:
+$options = "\
+-i {file} : DNA FASTA input file to investigate\
+-o {file} : Base for the names of output file(s)\
+-p {directory} : A path to a directory of HMM-profile collections representing ITS conserved regions, default is in the same directory as ITSx itself\
+--date {T or F} : Adds a date and time stamp to the output directory, off (F) by default\
+--reset {T or F} : Re-creates the HMM-database before ITSx is run, off (F) by default\
+
+Sequence selection options:\
+-t {character code} : Profile set to use for the search, see the User's Guide (comma-separated), default is all\
+-E {value} : Domain E-value cutoff for a sequence to be included in the output, default = 1e-5\
+-S {value} : Domain score cutoff for a sequence to be included in the output, default = 0\
+-N {value} : The minimal number of domains that must match a sequence before it is included, default = 2\
+--selection_priority {sum, domains, eval, score} : Selects what will be of highest priority when determining the origin of the sequence, default is sum\
+--search_eval {value} : The E-value cutoff used in the HMMER search, high numbers may slow down the process, cannot be used with the --search_score option, default is 0.01\
+--search_score {value} : The score cutoff used in the HMMER search, low numbers may slow down the process, cannot be used with the --search_eval option, default is to used E-value cutoff, not score\
+--allow_single_domain {e-value,score or F} : Allow inclusion of sequences that only find a single domain, given that they meet the given E-value and score thresholds, on with parameters 1e-9,0 by default\
+--allow_reorder {T or F} : Allows profiles to be in the wrong order on extracted sequences, off (F) by default\
+--complement {T or F} : Checks both DNA strands against the database, creating reverse complements, on (T) by default\
+--cpu {value} : the number of CPU threads to use, default is 1\
+--multi_thread {T or F} : Multi-thread the HMMER-search, on (T) if number of CPUs (--cpu option > 1), else off (F) by default\
+--heuristics {T or F} : Selects whether to use HMMER's heuristic filtering, off (F) by default\
+
+Output options:\
+--summary {T or F} : Summary of results output, on (T) by default\
+--graphical {T or F} : 'Graphical' output, on (T) by default\
+--fasta {T or F} : FASTA-format output of extracted ITS sequences, on (T) by default\
+--preserve {T or F} : Preserve sequence headers in input file instead of printing out ITSx headers, off (F) by default\
+--save_regions {SSU,ITS1,5.8S,ITS2,LSU,all,none} : A comma separated list of regions to output separate FASTA files for, 'ITS1,ITS2' by default\
+--anchor {integer or HMM} : Saves an additional number of bases before and after each extracted region. If set to 'HMM' all bases matching the corresponding HMM will be output, default = 0\
+--only_full {T or F} : If true, output is limited to full-length regions, off (F) by default\
+--partial {integer} : Saves additional FASTA-files for full and partial ITS sequences longer than the specified cutoff, default = 0 (off)\
+--concat {T or F} : Saves a FASTA-file with concatenated ITS sequences (with 5.8S removed), off (F) by default\
+--minlen {integer} : Minimum length the ITS regions must be to be outputted in the concatenated file (see above), default = 0\
+--positions {T or F} : Table format output containing the positions ITS sequences were found in, on (T) by default\
+--table {T or F} : Table format output of sequences containing probable ITS sequences, off (F) by default\
+--not_found {T or F} : Saves a list of non-found entries, on (T) by default\
+--detailed_results {T or F} : Saves a tab-separated list of all results, off (F) by default\
+--truncate {T or F} : Truncates the FASTA output to only contain the actual ITS sequences found, on (T) by default\
+--silent {T or F} : Supresses printing progress info to stderr, off (F) by default\
+--graph_scale {value} : Sets the scale of the graph output, if value is zero, a percentage view is shown, default = 0\
+--save_raw {T or F} : Saves all raw data for searches etc. instead of removing it on finish, off (F) by default\
+
+-h : displays this help message\
+--help : displays this help message\
+--bugs : displays the bug fixes and known bugs in this version of ITSx\
+--license : displays licensing information\
+";
+
+
+## Print title message
+print STDERR "$app_title\nby $app_author\nVersion: $app_version\n$app_message";
+print STDERR "-----------------------------------------------------------------\n";
+
+## Setup default variable values
+use List::Util qw(first max maxstr min minstr reduce shuffle sum);
+
+$bindir = $0;
+$bindir =~ s/_x//;
+$input = "";
+$output = "ITSx_out";
+$hmmscan = "";
+$profileDB = "$bindir\_db/HMMs";
+$type = "all";
+$E = 1e-5;
+$S = 0;
+$N = 2;
+$priority = "sum";
+$search_eval = 0.01;
+$search_score = "";
+$allow_single_E = 1e-9;
+$allow_single_score = 0;
+#$allow_single_E = -1; # Turns off single-domain matching by E-value
+#$allow_single_score = 0; # Turns off single-domain matching by score
+$allow_reorder = 0;
+$complement = 1;
+$cpu = 1;
+$multi_thread = "unset";
+$heuristics = 0;
+$out_sum = 1;
+$out_graph = 1;
+$out_fasta = 1;
+$out_preserve = 0;
+$out_ssu = 0;
+$out_its1 = 1;
+$out_its2 = 1;
+$out_58S = 0;
+$out_lsu = 0;
+$out_pos = 1;
+$out_table = 0;
+$out_not = 1;
+$out_date = 0;
+$out_joined = 0;
+$out_results = 0;
+$out_partial = 0;
+$out_concat = 0;
+$concat_minlen = 0;
+$truncate = 1;
+$anchor = 0;
+$only_full = 0;
+$graph_scale = 0;
+$debug = 0;
+$reset = 0;
+
+## Read command-line options
+for ($i = 0; $i <= scalar(@ARGV); $i++) { # Goes through the list of arguments
+ $arg = @ARGV[$i]; # Stores the current argument in $arg
+
+ if ($arg eq "-i") { # Read input files from -i flag
+ $i++;
+ $input = @ARGV[$i];
+ }
+ if ($arg eq "-o") { # Read output files from -o flag
+ $i++;
+ $output = @ARGV[$i];
+ }
+ if ($arg eq "-p") { # Read profile database from -p flag
+ $i++;
+ $profileDB = @ARGV[$i];
+ }
+ if ($arg eq "--hmmscan") { # Read pre-computed hmmscan output file from --hmmscan flag ('undocumented' feature)
+ $i++;
+ $hmmscan = @ARGV[$i];
+ }
+ if ($arg eq "--date") { # Determine whether or not to add a date stamp based on the --date flag
+ $i++;
+ if (substr(@ARGV[$i],0,1) =~ m/^[Ff0]/) { # Check if argument begins with "F", "f", or "0"
+ $out_date = 0;
+ } else {
+ $out_date = 1;
+ }
+ }
+ if ($arg eq "--reset") { # Reset HMM database?
+ $i++;
+ if (substr(@ARGV[$i],0,1) =~ m/^[Ff0]/) { # Check if argument begins with "F", "f", or "0"
+ $reset = 0;
+ } else {
+ $reset = 1;
+ }
+ }
+
+ if ($arg eq "-t") { # Select what types of ITSs to look for using the -t flag
+ $i++;
+ $type = @ARGV[$i];
+ }
+ if ($arg eq "-E") { # Set the E-value cutoff using the -E flag
+ $i++;
+ $E = @ARGV[$i];
+ }
+ if ($arg eq "-S") { # Set the score cutoff using the -S flag
+ $i++;
+ $S = @ARGV[$i];
+ }
+ if ($arg eq "-N") { # Set the number of found domains cutoff using the -N flag
+ $i++;
+ $N = @ARGV[$i];
+ }
+ if ($arg eq "--selection_priority") { # Set how to order the ITS types using the --selection_priority flag
+ $i++;
+ $priority = @ARGV[$i];
+ }
+ if ($arg eq "--search_eval") { # Set the E-value cutoff for the HMMER search using the --search_eval flag
+ $i++;
+ $search_eval = @ARGV[$i];
+ $search_score = ""; # Turns off score cutoff
+ }
+ if ($arg eq "--search_score") { # Set the score cutoff for the HMMER search using the --search_score flag
+ $i++;
+ $search_score = @ARGV[$i];
+ $search_eval = ""; # Turns off E-value cutoff
+ }
+ if ($arg eq "--allow_single_domain") { # Determine whether or not to allow single domain matches based on the --allow_single_domain flag
+ $i++;
+ if (substr(@ARGV[$i],0,1) =~ m/^[Ff0]/) { # Check if argument begins with "F", "f", or "0"
+ $allow_single_E = -1; # Turns off single-domain matching by E-value
+ $allow_single_score = 0; # Turns off single-domain matching by score
+ } else {
+ ($allow_single_E,$allow_single_score) = split(',',@ARGV[$i]); # Turns on single-domain matching, assigning the first given value as the E-value cutoff, and the second as score cutoff
+ }
+ }
+ if ($arg eq "--allow_reorder") { # Determine whether or not to allow the domains to be in the wrong order based on the --allow_reorder flag
+ $i++;
+ if (substr(@ARGV[$i],0,1) =~ m/^[Tt1]/) { # Check if argument begins with "T", "t", or "1"
+ $allow_reorder = 1;
+ } else {
+ $allow_reorder = 0;
+ }
+ }
+ if ($arg eq "--complement") { # Determine whether or not to scan the complementary strand of the input file based on the --complement flag
+ $i++;
+ if (substr(@ARGV[$i],0,1) =~ m/^[Tt1]/) { # Check if argument begins with "T", "t", or "1"
+ $complement = 1;
+ } else {
+ $complement = 0;
+ }
+ }
+ if ($arg eq "--cpu") { # Set the number of CPUs to use based on the --cpu flag
+ $i++;
+ $cpu = @ARGV[$i];
+ }
+ if ($arg eq "--multi_thread") { # Determine whether or not to multi-thread the HMMER step based on the --multi_thread flag
+ $i++;
+ if (substr(@ARGV[$i],0,1) =~ m/^[Tt1]/) { # Check if argument begins with "T", "t", or "1"
+ $multi_thread = 1;
+ } else {
+ $multi_thread = 0;
+ }
+ }
+ if ($arg eq "--heuristics") { # Determine whether or not to use HMMER's heuristic filtering based on the --heuristics flag
+ $i++;
+ if (substr(@ARGV[$i],0,1) =~ m/^[Tt1]/) { # Check if argument begins with "T", "t", or "1"
+ $heuristics = 1;
+ } else {
+ $heuristics = 0;
+ }
+ }
+
+ if ($arg eq "--summary") { # Determine whether or not to output a summary based on the --summary flag
+ $i++;
+ if (substr(@ARGV[$i],0,1) =~ m/^[Tt1]/) { # Check if argument begins with "T", "t", or "1"
+ $out_sum = 1;
+ } else {
+ $out_sum = 0;
+ }
+ }
+ if ($arg eq "--graphical") { # Determine whether or not to output a graphical representation of matches based on the --graphical flag
+ $i++;
+ if (substr(@ARGV[$i],0,1) =~ m/^[Tt1]/) { # Check if argument begins with "T", "t", or "1"
+ $out_graph = 1;
+ } else {
+ $out_graph = 0;
+ }
+ }
+ if ($arg eq "--detailed_results") { # Determine whether or not to output a detailed results list, based on the --detailed_results flag
+ $i++;
+ if (substr(@ARGV[$i],0,1) =~ m/^[Tt1]/) { # Check if argument begins with "T", "t", or "1"
+ $out_results = 1;
+ } else {
+ $out_results = 0;
+ }
+ }
+ if ($arg eq "--partial") { # Set the full-and-partial cutoff
+ $i++;
+ $out_partial = @ARGV[$i];
+ }
+ if ($arg eq "--anchor") { # Set the length of the sequence "anchors"
+ $i++;
+ $anchor = @ARGV[$i];
+ }
+ if ($arg eq "--only_full") { # Output only full-length regions
+ $i++;
+ if (substr(@ARGV[$i],0,1) =~ m/^[Tt1]/) { # Check if argument begins with "T", "t", or "1"
+ $only_full = 1;
+ } else {
+ $only_full = 0;
+ }
+ }
+ if ($arg eq "--save_regions") { # Determine which regions to output FASTA files for based on the --save_regions flag
+ $i++;
+ @save_regions = split(',',uc(@ARGV[$i]));
+ $out_ssu = 0;
+ $out_its1 = 0;
+ $out_its2 = 0;
+ $out_58S = 0;
+ $out_lsu = 0;
+ foreach $save_region (@save_regions) {
+ if ($save_region eq "SSU") {
+ $out_ssu = 1;
+ }
+ if ($save_region eq "ITS1") {
+ $out_its1 = 1;
+ }
+ if ($save_region eq "5.8S") {
+ $out_58S = 1;
+ }
+ if ($save_region eq "ITS2") {
+ $out_its2 = 1;
+ }
+ if ($save_region eq "LSU") {
+ $out_lsu = 1;
+ }
+ if ($save_region eq "ALL") {
+ $out_ssu = 1;
+ $out_its1 = 1;
+ $out_its2 = 1;
+ $out_58S = 1;
+ $out_lsu = 1;
+ }
+ if ($save_region eq "NONE") {
+ $out_ssu = 0;
+ $out_its1 = 0;
+ $out_its2 = 0;
+ $out_58S = 0;
+ $out_lsu = 0;
+ }
+ }
+ }
+ if ($arg eq "--positions") { # Determine whether or not to output a positions file based on the --positions flag
+ $i++;
+ if (substr(@ARGV[$i],0,1) =~ m/^[Tt1]/) { # Check if argument begins with "T", "t", or "1"
+ $out_pos = 1;
+ } else {
+ $out_pos = 0;
+ }
+ }
+ if ($arg eq "--concat") { # Determine whether or not to output a concatednated ITS1 + ITS2 file
+ $i++;
+ if (substr(@ARGV[$i],0,1) =~ m/^[Tt1]/) { # Check if argument begins with "T", "t", or "1"
+ $out_concat = 1;
+ } else {
+ $out_concat = 0;
+ }
+ }
+ if ($arg eq "--minlen") { # Set the min length of the combined ITS1 and ITS2 sequences for concatenation
+ $i++;
+ $concat_minlen = @ARGV[$i];
+ }
+ if ($arg eq "--fasta") { # Determine whether or not to output FASTA-files based on the --fasta flag
+ $i++;
+ if (substr(@ARGV[$i],0,1) =~ m/^[Tt1]/) { # Check if argument begins with "T", "t", or "1"
+ $out_fasta = 1;
+ } else {
+ $out_fasta = 0;
+ }
+ }
+ if ($arg eq "--preserve") { # Determine whether or not to preserve FASTA-headers based on the --preserve flag
+ $i++;
+ if (substr(@ARGV[$i],0,1) =~ m/^[Tt1]/) { # Check if argument begins with "T", "t", or "1"
+ $out_preserve = 1;
+ } else {
+ $out_preserve = 0;
+ }
+ }
+ if ($arg eq "--joined") { # Determine whether or not to output a FASTA-file containing ALL sorts of output sequences (for debugging)
+ $i++;
+ if (substr(@ARGV[$i],0,1) =~ m/^[Tt1]/) { # Check if argument begins with "T", "t", or "1"
+ $out_joined = 1;
+ } else {
+ $out_joined = 0;
+ }
+ }
+ if ($arg eq "--table") { # Determine whether or not to output tables of all potential matches based on the --table flag
+ $i++;
+ if (substr(@ARGV[$i],0,1) =~ m/^[Tt1]/) { # Check if argument begins with "T", "t", or "1"
+ $out_table = 1;
+ } else {
+ $out_table = 0;
+ }
+ }
+ if ($arg eq "--not_found") { # Determine whether or not to output a list of sequences that are not ITSs based on the --not_found flag
+ $i++;
+ if (substr(@ARGV[$i],0,1) =~ m/^[Tt1]/) { # Check if argument begins with "T", "t", or "1"
+ $out_not = 1;
+ } else {
+ $out_not = 0;
+ }
+ }
+ if ($arg eq "--silent") { # Determine whether or not to output anything to the screen based on the --silent flag
+ $i++;
+ if (substr(@ARGV[$i],0,1) =~ m/^[Tt1]/) { # Check if argument begins with "T", "t", or "1"
+ $silent = 1;
+ } else {
+ $silent = 0;
+ }
+ }
+ if ($arg eq "--graph_scale") { # Set the scale of the graphical output based on the --graph_scale flag
+ $i++;
+ $graph_scale = @ARGV[$i];
+ }
+ if ($arg eq "--save_raw") { # Determine whether or not to save all the raw intermediate data based on the --save_raw flag
+ $i++;
+ if (substr(@ARGV[$i],0,1) =~ m/^[Tt1]/) { # Check if argument begins with "T", "t", or "1"
+ $save_raw = 1;
+ } else {
+ $save_raw = 0;
+ }
+ }
+
+ ## If "-h" or "--help" are among the options, output usage data and options
+ if (($arg eq "-h") || ($arg eq "--help")) {
+ print "Usage: ITSx -i \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 = ) { # 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 = ) { # 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 = ) { # 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 = ) { # 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 = )) {
+ 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
diff -r 000000000000 -r 57e76d405a97 ITSx.xml
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/ITSx.xml Mon Mar 14 16:15:08 2016 -0400
@@ -0,0 +1,173 @@
+
+
+
+ ITSx -- Identifies ITS sequences and extracts the ITS region
+
+
+
+ hmmer
+
+
+
+
+
+
+ ITSx -i $input
+ --detailed_results T
+ -E $domain_eval
+ -S $domain_score
+ -N $domain_num
+ --cpu=\${GALAXY_SLOTS:-1}
+ $allow_reorder
+ $complement
+ $heuristics
+ $preserve
+ $reset
+ $truncate
+ #if $hmmer_search.search == "eval"
+ --search_eval $hmmer_search.eval_cutoff
+ #else
+ --search_score $hmmer_search.score_cutoff
+ #end if
+
+
+
+
+ Search E-value
+ Search Score
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+**Description**
+
+Identifies ITS sequences and extracts the ITS regions.
+
+ITSx is an open source software utility to extract the highly variable ITS1 and ITS2 subregions from ITS sequences, which is commonly used as a molecular barcode for e.g. fungi. As the inclusion of parts of the neighbouring, very conserved, ribosomal genes (SSU, 5S and LSU rRNA sequences) in the sequence identification process can lead to severely misleading results, ITSx identifies and extracts only the ITS regions themselves.
+For more information regarding the settings of the tool, please visit the ITSx Users Guide on http://microbiology.se/publ/itsx_users_guide.pdf
+
+-----
+
+**Inputs and Outputs**
+
+- Inputs:
+ + Accepts input in the FASTA format. It is possible to input both aligned and unaligned FASTA files, containing both DNA and RNA sequences
+
+- Outputs:
+ + Summary file of the entire run
+ + One or more detailed table containing the positions in the respective sequences where the ITS subregions were found
+ + Semi-graphical representation of hits
+ + FASTA file of all identified sequences
+ + FASTA file for ITS1 regions
+ + FASTA file for ITS2 regions
+ + A list of sequence IDs representing the entries (if they exist) that did not contain any ITS region
+ + A file containing chimeric sequences (if they are found)
+ + A file for problematic sequences
+
+-----
+
+**Tool Information**
+
+- ITSx tool v1.0.11:
+ + Release Date: December 2014
+ + URL: http://microbiology.se/software/itsx/
+ + Source URL: http://microbiology.se/sw/ITSx_1.0.11.tar.gz
+ + Copyright (C) 2012-2013 Johan Bengtsson-Palme et al.
+
+-----
+
+**Contact Information**
+
+- For more information, contact:
+ + AAFC-MBB Team
+ + E-mail: mbb@agr.gc.ca
+ + Agriculture and Agri-foods Canada, Ottawa, ON, Canada
+
+
+
+ @ARTICLE{a1,
+ author = {Bengtsson-Palme, Johan},
+ title = {ITSx: Improved software detection and extraction of ITS1 and ITS2 from ribosomal ITS sequences of fungi and other eukaryotes for use in environmental sequencing},
+ year = {2013}
+ }
+
+