fix bug on macos, h1s works with fas >= 1.4

Author: trvinh

h1s/bin/hamstr.pl

@@ -1254,7 +1254,7 @@ sub checkInput {
 
 	## 13) setting up the directories where the output files will be put into.
 	$fa_dir_neu = $outpath . '/fa_dir_' . $dbfile_short . '_' . $hmmset . '_' . $refspec[0];
-	$tmpdir = $outpath . '/' . $tmpdir;
+	# $tmpdir = $tmpdir;
 	if (!$strict) {
 		$fa_dir_neu = $outpath . '/fa_dir_' . $dbfile_short . '_' . $hmmset . '_' . $refspec[0];
 	}
@@ -1281,10 +1281,10 @@ sub checkInput {
 			`rm -rf "$fa_dir_neu"`;
 			`mkdir "$fa_dir_neu"`;
 		}
-		if (!(-e "$tmpdir")) {
+		if (!(-d "$tmpdir")) {
 			`mkdir "$tmpdir"`;
 		}
-		elsif (-e "$tmpdir" and $cleartmp) {
+		elsif (-d "$tmpdir" and $cleartmp) {
 			`rm -rf "$tmpdir"`;
 			`mkdir "$tmpdir"`;
 		}

h1s/bin/oneSeq.pl

@@ -169,6 +169,9 @@
 my $annotation_prog = "annoFAS";
 my $fas_prog = "calcFAS";
 my $hamstrFAS_prog = "hamstrFAS";
+if (check_exists_command("fdogFAS")) {
+	$hamstrFAS_prog = "fdogFAS";
+}
 
 ##### ublast Baustelle: not implemented yet
 my $runublast = 0;
@@ -607,6 +610,7 @@
 	my $final_eval_blast = $eval_blast*$eval_relaxfac;
 	my $final_eval_hmmer = $eval_hmmer*$eval_relaxfac;
 
+	$taxaPath = $genome_dir;
 	my @searchTaxa;
 	unless($groupNode) {
 		@searchTaxa = keys %taxa;
@@ -643,14 +647,16 @@
 	foreach (sort @searchTaxa) {
 		chomp(my $searchTaxon = $_);
 		my $pid = $pm->start and next;
-		my $doneTaxon = runHamstr($searchTaxon, $seqName, $finalOutput, $refSpec, $hitlimit, $representative, $strict, $coremode, $final_eval_blast, $final_eval_hmmer, $aln);
-		unless ($silent) {
-			print $doneTaxon,"\t";
-		}
-		my $doneTaxonName = getTaxonName($doneTaxon);
-		if (defined($doneTaxonName)) {
-			print $doneTaxonName, " DONE\n";
+
+		my $searchTaxonName = getTaxonName($searchTaxon);
+		if (defined($searchTaxonName)) {
+			unless ($silent) {
+				print $searchTaxon, "\t", $searchTaxonName, "\n";
+			} else {
+				print $searchTaxonName, "\n";
+			}
 		}
+		runHamstr($searchTaxon, $seqName, $finalOutput, $refSpec, $hitlimit, $representative, $strict, $coremode, $final_eval_blast, $final_eval_hmmer, $aln);
 		$pm->finish;
 	}
 	$pm->wait_all_children;
@@ -1568,14 +1574,16 @@ sub getBestOrtholog {
 			unless ($silent) {
 				print "Hamstr species: " . $key->scientific_name . " - " . @{$key->name('supplied')}[0] . "\n";
 			}
-			my $doneTaxon = runHamstr(@{$key->name('supplied')}[0], $seqName, $outputFa, $refSpec, $core_hitlimit, $core_rep, $corestrict, $coremode, $eval_blast, $eval_hmmer, $aln);
-			unless ($silent) {
-				print $doneTaxon,"\t";
-			}
-			my $doneTaxonName = getTaxonName($doneTaxon);
-			if (defined($doneTaxonName)) {
-				print $doneTaxonName, " DONE\n";
+			my $coreTaxon = @{$key->name('supplied')}[0];
+			my $coreTaxonName = getTaxonName($coreTaxon);
+			if (defined($coreTaxonName)) {
+				unless ($silent) {
+					print $coreTaxon, "\t", $coreTaxonName, "\n";
+				} else {
+					print $coreTaxonName, "\n";
+				}
 			}
+			runHamstr($coreTaxon, $seqName, $outputFa, $refSpec, $core_hitlimit, $core_rep, $corestrict, $coremode, $eval_blast, $eval_hmmer, $aln);
 			## check weather a candidate was found in the searched taxon
 			if(-e $candidatesFile) {
 
@@ -1900,11 +1908,11 @@ sub getRefTree {
 sub getTaxonName {
 	my $taxAbbr = $_[0];
 	my @tmp = split(/@/,$taxAbbr);
-	my $taxon = $db_bkp->get_taxon($tmp[1]);
+	my $taxon = $db_bkp->get_taxon(-taxonid => $tmp[1]);
 	if (defined($taxon)) {
 		return($taxon->scientific_name);
 	} else {
-		return("Unk NCBI taxon");
+		return("Unk NCBI taxon for $taxAbbr");
 	}
 }
 
@@ -2043,7 +2051,6 @@ sub runHamstr {
 	else {
 		print "No protein set available for $taxon. Failed to fetch it from database and nothing at $taxaDir. Skipping!\n";
 	}
-	return($taxon);
 }
 
 # add seed sequence to output file if not exists
@@ -2062,19 +2069,30 @@ sub addSeedSeq {
 			}
 		}
 	}
-	# get seed seq
+	# get seed sequence and add to be beginning of the fasta output
+	open(TEMP, ">$outputFa.temp") or die "Cannot create $outputFa.temp!\n";
+	my $seedFa = "";
 	if ($flag == 1) {
+		# first, write the seed to TEMP
 		my $seqio = Bio::SeqIO->new(-file => "$coreOrthologsPath/$seqName/$seqName.fa", '-format' => 'Fasta');
 		while(my $seq = $seqio->next_seq) {
 			my $id = $seq->id;
 			if ($id =~ /$refSpec/) {
-				my $seedFa = ">".$id."|1\n".$seq->seq;
-				# append to begining of outputFa
-				my $headCommand = "sed -i \'1s/^/". ">".$id."|1\\n".$seq->seq . "\\n/\' " . $outputFa;
-				system($headCommand);
+				print TEMP ">$id|1\n", $seq->seq, "\n";
+				last;
+			}
+		}
+		# then write other sequences
+		my $seqio2 = Bio::SeqIO->new(-file => "$outputFa", '-format' => 'Fasta');
+		while(my $seq = $seqio2->next_seq) {
+			my $id = $seq->id;
+			if ($id !~ /$refSpec/) {
+				print TEMP ">$id\n", $seq->seq, "\n";
 			}
 		}
 	}
+	close(TEMP);
+	system("mv $outputFa.temp $outputFa")
 }
 
 ##########################
@@ -2549,7 +2567,7 @@ sub initialCheck {
 	}
 
 	# check executable FAS
-	my $fasCheckMsg = `prepareFAS -t ./ -c 2>&1`;
+	my $fasCheckMsg = `$hamstrFAS_prog -t ./ -c 2>&1`;
 	if ($fasoff != 1 && $fasCheckMsg =~ /ERROR/) {
 		die "ERROR: greedyFAS not ready to use! Please check https://github.com/BIONF/FAS/wiki/prepareFAS\n";
 	}
@@ -2642,6 +2660,11 @@ sub checkValidFolderName {
 sub gettime { sprintf"%d.%03d",Time::HiRes::gettimeofday }
 sub roundtime { sprintf("%.2f", $_[0]); }
 
+sub check_exists_command {
+    my $check = `sh -c 'command -v $_[0]'`;
+    return $check;
+}
+
 ###########################
 sub helpMessage {
 	my $helpmessage = "

h1s/h1s.py

@@ -161,7 +161,7 @@ def h1s(args):
         sys.exit('Problem running\n%s' % (cmd))
 
 def main():
-    version = '2.2.17'
+    version = '2.2.18'
     parser = argparse.ArgumentParser(description='You are running h1s version ' + str(version) + '.')
     parser.add_argument('--version', action='version', version=str(version))
     required = parser.add_argument_group('Required arguments')

h1s/hms.py

@@ -178,7 +178,7 @@ def calcFAS (outpath, extendedFa, weightpath, cpu):
         sys.exit('Problem running\n%s' % (fasCmd))
 
 def main():
-    version = '2.2.17'
+    version = '2.2.18'
     parser = argparse.ArgumentParser(description='You are running h1s version ' + str(version) + '.')
     parser.add_argument('--version', action='version', version=str(version))
     required = parser.add_argument_group('Required arguments')

setup.py

@@ -26,7 +26,7 @@
 
 setup(
     name="h1s",
-    version="2.2.17",
+    version="2.2.18",
     python_requires='>=3.7.0',
     description="Feature-aware orthology prediction tool",
     long_description=long_description,