Add hashing to plot-tree

Author: reslp

bin/libphylociraptor/usageinfo.py

@@ -251,10 +251,13 @@ def hi():
 	-i, --intrees		Relative paths to input trees in results folder, seperated by commas. (Default: all)
 	-o, --outprefix		Output file name prefix.
 	-n, --nquartets		Total number of quartets to be calculated. Use either this or --stopby.
+	OR
 	-b, --stopby		Stoping criterion. There are two options:
 					-b conflicts=100 stops as soon as 100 conflicts have been found.
 					--stopby tipcoverage=100 stops as soon as every tip is in at least 100 quartets.
+
 Optional Arguments:
+	--config-file           Relative custom config-file path. (Default: data/config.yaml)
 	-s, --seed		Random seed number for reproducibility. (Default: random)
 	-l, --lineagefile	Lineagefile created with phylociraptor util get-lineage.
 				Mandatory when using-a/--selecttaxa with a specific taxonomic level.
@@ -278,6 +281,7 @@ def hi():
 				Use 'all' to plot all available trees.
 	
 Optional Arguments:
+	--config-file           Relative custom config-file path. (Default: data/config.yaml)
 	-o, --outprefix		Output file name prefix.
 	-l, --lineagefile	Lineagefile created with phylociraptor util get-lineage.
 	-e, --level		Taxonomic level in lineage file which should be plotted.

bin/plot-tree.R

@@ -283,11 +283,13 @@ if (level == "none" || lineage_file == "none") {
   for (ntree in 1:length(treenames)) {
     #extract filename information:
     treename <- treenames[ntree]
-    bs_cutoff <- strsplit(strsplit(treename,"/")[[1]][2],"-")[[1]][2]
+
+    bs_cutoff <- strsplit(strsplit(treename,"cutoff-")[[1]][2],"/")[[1]][1]
     algorithm <- strsplit(treename,"/")[[1]][3]
-    alitrim <- strsplit(treename,"/")[[1]][4]
+    alitrim <- strsplit(strsplit(treename,".",fixed=T)[[1]][1], "/")[[1]][5]
+    hash <- strsplit(strsplit(treename,".",fixed=T)[[1]][2], "/")[[1]][1]
     prefix <- paste( algorithm, alitrim, bs_cutoff, sep="-")
-    
+
     #reroot tree
     tree <- read.tree(treename)
     if (outgroup != "none") { #reroot tree in case an outgroup was specified
@@ -308,7 +310,7 @@ if (level == "none" || lineage_file == "none") {
     # extract legend then remove it
     cat(paste0("    Write PDF: ",prefix,"-",level,"-tree.pdf\n"))
     pdf(file = paste0(prefix,"-",level,"-tree.pdf"), width = 10, height=get_pdf_height(tree))
-    print(t2 + theme(legend.position="none") +  plot_annotation(title = prefix, caption=paste0("Taxonomy: ", level,". Random seed:", seed,".\nOutgroup: ", outgroups), theme=theme(plot.title=element_text(hjust=0.5, size=16))))#+ plot_layout(guides = 'none')# & theme(legend.position='bottom')
+    print(t2 + theme(legend.position="none") +  plot_annotation(title = prefix, caption=paste0("Taxonomy: ", level,". Random seed: ", seed,".\nOutgroup: ", outgroups, "\nHash: ", hash), theme=theme(plot.title=element_text(hjust=0.5, size=16))))#+ plot_layout(guides = 'none')# & theme(legend.position='bottom')
     garbage <- dev.off()
     if (singlet == TRUE) {break}
   }
@@ -337,9 +339,10 @@ if (level == "none" || lineage_file == "none") {
   for (ntree in 1:length(treenames)) {
     #extract filename information:
     treename <- treenames[ntree]
-    bs_cutoff <- strsplit(strsplit(treename,"/")[[1]][2],"-")[[1]][2]
+    bs_cutoff <- strsplit(strsplit(treename,"cutoff-")[[1]][2],"/")[[1]][1]
     algorithm <- strsplit(treename,"/")[[1]][3]
-    alitrim <- strsplit(treename,"/")[[1]][4]
+    alitrim <- strsplit(strsplit(treename,".",fixed=T)[[1]][1], "/")[[1]][5]
+    hash <- strsplit(strsplit(treename,".",fixed=T)[[1]][2], "/")[[1]][1]
     prefix <- paste( algorithm, alitrim, bs_cutoff, sep="-")
     cat(paste0("Plot tree: ", prefix, "\n"))
     #reroot tree
@@ -469,7 +472,7 @@ if (level == "none" || lineage_file == "none") {
     if (single == "yes") {
         print(t2 + theme(legend.position="none") + plot_annotation(title = prefix, caption=paste0("Taxonomy: ", level,". Random seed:", seed,".\nOutgroup: ", outgroups), theme=theme(plot.title=element_text(hjust=0.5, size=16))))#+ plot_layout(guides = 'none')# & theme(legend.position='bottom')
     } else {
-        print(t2 + t1 + theme(legend.position="none") +  plot_annotation(tag_levels = 'A', title = prefix, caption=paste0("A) Topology without branch lengths. B) Topology with branch lengths and collapsed taxonomic groups. Taxonomy: ", level,". Random seed:", seed,".\nOutgroup: ", outgroups), theme=theme(plot.title=element_text(hjust=0.5, size=16))) + plot_layout(design = layout))#+ plot_layout(guides = 'none')# & theme(legend.position='bottom')
+        print(t2 + t1 + theme(legend.position="none") +  plot_annotation(tag_levels = 'A', title = prefix, caption=paste0("A) Topology without branch lengths. B) Topology with branch lengths and collapsed taxonomic groups. Taxonomy: ", level,". Random seed: ", seed,".\nOutgroup: ", outgroups, "\nHash: ", hash), theme=theme(plot.title=element_text(hjust=0.5, size=16))) + plot_layout(design = layout))#+ plot_layout(guides = 'none')# & theme(legend.position='bottom')
     } 
     garbage <- dev.off()
     if (single == TRUE) {break}

phylociraptor

@@ -889,6 +889,7 @@ elif args.command == "util":
 		qs_parser.add_argument("-l", "--lineagefile", action="store")
 		qs_parser.add_argument("-b", "--stopby", action="store", default=None)
 		qs_parser.add_argument("-a", "--selecttaxa", action="store", default="")
+		qs_parser.add_argument("--config-file", dest="configfile", action="store", default="data/config.yaml")
 		qs_parser.add_argument("--quiet", action="store_true", default=False) 
 		qs_args = qs_parser.parse_args(args.arguments) 
 		if qs_args.help: 
@@ -904,10 +905,34 @@ elif args.command == "util":
 			sys.exit(1)
 		elif qs_args.intrees == "all":
 			print(now(), "Will compare all trees, this expects phylociraptor mltree and/or speciestree to be finished.")
-			iqtree_trees = glob.glob("results/phylogeny-*/iqtree/*/concat.contree")
-			raxml_trees =  glob.glob("results/phylogeny-*/raxml/*/raxmlng.raxml.support")
-			astral_trees = glob.glob("results/phylogeny-*/astral/*/species_tree.tre")
-			all_trees = ",".join(iqtree_trees + raxml_trees + astral_trees)
+			print(now(), "Calculating analysis hashes for mltree, this may take a few seconds...")
+			hashes = {}
+			hashes["mltree"] = collect_hashes("mltree", parse_config_file(qs_args.configfile), qs_args.configfile, debug=qs_args.debug, wd=os.getcwd())
+			print(now(), "Calculating analysis hashes for speciestree, this may take a few seconds...")
+			hashes["speciestree"] = collect_hashes("speciestree", parse_config_file(qs_args.configfile), qs_args.configfile, debug=qs_args.debug, wd=os.getcwd())	
+			combinations = []
+			#first calculate all possible combinations:
+			for bs in hashes["mltree"]["mltree"]["per"]:
+				for treeprog in hashes["mltree"]["mltree"]["per"][bs]:
+					for trimmer in hashes["mltree"]["mltree"]["per"][bs][treeprog]["iqtree"]:
+						for aligner in hashes["mltree"]["mltree"]["per"][bs][treeprog]["iqtree"][trimmer]:
+							if treeprog == "iqtree":
+								filename = "concat.contree"
+							elif treeprog == "raxml":
+								filename = "raxmlng.raxml.support"
+							combinations.append("results/phylogeny/" + treeprog + "/bootstrap-cutoff-" + bs + "/" + aligner + "-" + trimmer + "." + hashes["mltree"]["mltree"]["per"][bs][treeprog]["iqtree"][trimmer][aligner] + "/" + filename)
+			for bs in hashes["speciestree"]["speciestree"]["per"]:
+				for treeprog in hashes["speciestree"]["speciestree"]["per"][bs]:
+					for trimmer in hashes["speciestree"]["speciestree"]["per"][bs][treeprog]["iqtree"]:
+						for aligner in hashes["speciestree"]["speciestree"]["per"][bs][treeprog]["iqtree"][trimmer]:
+							combinations.append("results/phylogeny/" + treeprog + "/bootstrap-cutoff-" + bs + "/" + aligner + "-" + trimmer + "." + hashes["speciestree"]["speciestree"]["per"][bs][treeprog]["iqtree"][trimmer][aligner] + "/species_tree.tre")
+			# now check which files are there in the file system:
+			iqtree_trees = glob.glob("results/phylogeny/iqtree/bootstrap-cutoff-*/*/concat.contree")
+			raxml_trees =  glob.glob("results/phylogeny/raxml/bootstrap-cutoff-*/*/raxmlng.raxml.support")
+			astral_trees = glob.glob("results/phylogeny/astral/bootstrap-cutoff-*/*/species_tree.tre")
+			all_trees = iqtree_trees + raxml_trees + astral_trees
+			all_trees = [tree for tree in all_trees if tree in combinations] #keep only files which are specified in config file.
+			all_trees = ",".join(all_trees)
 		else:
 			all_trees = qs_args.intrees	
 		cmd = ["singularity", "exec", "-B", os.getcwd(), "docker://reslp/phylopy:2", "python", "bin/estimate_conflict.py", "-i", all_trees,  "-o", qs_args.outfile, "-s", qs_args.seed, "-t", qs_args.threads]
@@ -937,6 +962,7 @@ elif args.command == "util":
 		qs_parser.add_argument("-l", "--lineagefile", action="store", default="none")
 		qs_parser.add_argument("-e", "--level", action="store", default="none")
 		qs_parser.add_argument("-s", "--seed", action="store", default="random") 
+		qs_parser.add_argument("--config-file", dest="configfile", action="store", default="data/config.yaml")
 		qs_parser.add_argument("--quiet", action="store_true", default=False) 
 		qs_parser.add_argument("--single", action="store_true", default=False)
 		qs_args = qs_parser.parse_args(args.arguments) 
@@ -949,10 +975,34 @@ elif args.command == "util":
 			sys.exit(1)
 		elif qs_args.intrees == "all":
 			print(now(), "Will compare all trees, this expects phylociraptor mltree and/or speciestree to be finished.")
-			iqtree_trees = glob.glob("results/phylogeny-*/iqtree/*/concat.contree")
-			raxml_trees =  glob.glob("results/phylogeny-*/raxml/*/raxmlng.raxml.support")
-			astral_trees = glob.glob("results/phylogeny-*/astral/*/species_tree.tre")
-			all_trees = ",".join(iqtree_trees + raxml_trees + astral_trees)
+			print(now(), "Calculating analysis hashes for mltree, this may take a few seconds...")
+			hashes = {}
+			hashes["mltree"] = collect_hashes("mltree", parse_config_file(qs_args.configfile), qs_args.configfile, debug=qs_args.debug, wd=os.getcwd())
+			print(now(), "Calculating analysis hashes for speciestree, this may take a few seconds...")
+			hashes["speciestree"] = collect_hashes("speciestree", parse_config_file(qs_args.configfile), qs_args.configfile, debug=qs_args.debug, wd=os.getcwd())	
+			combinations = []
+			#first calculate all possible combinations:
+			for bs in hashes["mltree"]["mltree"]["per"]:
+				for treeprog in hashes["mltree"]["mltree"]["per"][bs]:
+					for trimmer in hashes["mltree"]["mltree"]["per"][bs][treeprog]["iqtree"]:
+						for aligner in hashes["mltree"]["mltree"]["per"][bs][treeprog]["iqtree"][trimmer]:
+							if treeprog == "iqtree":
+								filename = "concat.contree"
+							elif treeprog == "raxml":
+								filename = "raxmlng.raxml.support"
+							combinations.append("results/phylogeny/" + treeprog + "/bootstrap-cutoff-" + bs + "/" + aligner + "-" + trimmer + "." + hashes["mltree"]["mltree"]["per"][bs][treeprog]["iqtree"][trimmer][aligner] + "/" + filename)
+			for bs in hashes["speciestree"]["speciestree"]["per"]:
+				for treeprog in hashes["speciestree"]["speciestree"]["per"][bs]:
+					for trimmer in hashes["speciestree"]["speciestree"]["per"][bs][treeprog]["iqtree"]:
+						for aligner in hashes["speciestree"]["speciestree"]["per"][bs][treeprog]["iqtree"][trimmer]:
+							combinations.append("results/phylogeny/" + treeprog + "/bootstrap-cutoff-" + bs + "/" + aligner + "-" + trimmer + "." + hashes["speciestree"]["speciestree"]["per"][bs][treeprog]["iqtree"][trimmer][aligner] + "/species_tree.tre")
+			# now check which files are there in the file system:
+			iqtree_trees = glob.glob("results/phylogeny/iqtree/bootstrap-cutoff-*/*/concat.contree")
+			raxml_trees =  glob.glob("results/phylogeny/raxml/bootstrap-cutoff-*/*/raxmlng.raxml.support")
+			astral_trees = glob.glob("results/phylogeny/astral/bootstrap-cutoff-*/*/species_tree.tre")
+			all_trees = iqtree_trees + raxml_trees + astral_trees
+			all_trees = [tree for tree in all_trees if tree in combinations] #keep only files which are specified in config file.	
+			all_trees = ",".join(all_trees)
 		else:
 			all_trees = qs_args.intrees	
 		if qs_args.single: