update

Author: sinamajidian

FastOMA/_config.py

@@ -28,16 +28,19 @@
 threshold_dubious_sd = 1/10
 overlap_fragments = 0.15
 
+mergHOG_ratioMax_thresh = 0.8
+mergHOG_ratioMin_thresh = 0.9
+mergHOG_shared_thresh = 50
 
 # VARIABLE_threshold_dubious_sd
 # threshold_dubious_sd = float(os.getenv(VARIABLE_threshold_dubious_sd, 0.1))
 
 ## output writing files
-gene_trees_write = False
-msa_write = False
-gene_trees_write_all = False
-msa_write_all = False
-keep_subhog_each_pickle = False
+gene_trees_write = True
+msa_write = True
+gene_trees_write_all = True
+msa_write_all = True
+keep_subhog_each_pickle = True
 
 big_rhog_size = 60 * 1000
 omamer_family_threshold = 90

FastOMA/_utils_frag_SO_detection.py

@@ -37,7 +37,7 @@ def read_msa(input_msa):
         if len(ii):
             coords.append((np.min(ii), np.max(ii)))
         else:
-            logger_hog.warning("issue 1321230 all of the seq is gap"+str(rec))
+            logger_hog.warning("issue 1321230 all of the seq is gap "+str(rec))
             coords.append((0, 0))
 
     ids = np.array(ids)
@@ -219,7 +219,10 @@ def handle_fragment_sd(node_species_tree, gene_tree, genetree_msa_file_addr, all
         if prot_dubious_sd_remove_list:
             rest_leaves = set([i.name for i in gene_tree.get_leaves()]) - set(prot_dubious_sd_remove_list)
             gene_tree.prune(rest_leaves, preserve_branch_length=True)
-            (gene_tree, all_species_dubious_sd_dic_updated) = _utils_subhog.genetree_sd(node_species_tree, gene_tree, genetree_msa_file_addr + "_dubious_sd"+str(itr_so))
+            if _config.gene_trees_write_all or _config.rooting_method=="mad":
+                gene_tree.write(format=1, outfile=genetree_msa_file_addr + "_dubious_sd"+str(itr_so)+".nwk" )
+
+            (gene_tree, all_species_dubious_sd_dic_updated) = _utils_subhog.genetree_sd(node_species_tree, gene_tree, genetree_msa_file_addr + "_dubious_sd"+str(itr_so)+".nwk")
 
             hogs_children_level_list_raw = hogs_children_level_list
             for prot_dubious_sd_remove in prot_dubious_sd_remove_list:
@@ -409,8 +412,11 @@ def handle_fragment_msa(prot_dubious_msa_list, seq_dubious_msa_list, gene_tree,
             else:
                 gene_tree.prune(rest_leaves, preserve_branch_length=True)
 
+        if _config.gene_trees_write_all or _config.rooting_method == "mad":
+            gene_tree.write(outfile=genetree_msa_file_addr+"_dubiousMSA.nwk",format=1)
+
         if len(gene_tree) > 1:
-            (gene_tree, all_species_dubious_sd_dic2) = _utils_subhog.genetree_sd(node_species_tree, gene_tree, genetree_msa_file_addr+"_dubiousMSA")
+            (gene_tree, all_species_dubious_sd_dic2) = _utils_subhog.genetree_sd(node_species_tree, gene_tree, genetree_msa_file_addr+"_dubiousMSA.nwk")
 
 
 

FastOMA/_utils_roothog.py

@@ -136,6 +136,7 @@ def group_prots_roothogs(hogmaps):
         for prot_id, prot_map in prots_map.items():
             # omamer output is sorted based on normcount. but that's ok
             # this helps me in other functions like handle_singleton in this
+            #  this should be commented
             # if len(prot_map)>1:
             #     scores = [float(i[1]) for i in prot_map] # (hogid,score,seqlen,subfamily_medianseqlen)
             #     rhogids =[i[0] for i in prot_map]
@@ -328,7 +329,7 @@ def write_rhog(rhogs_prot_records, prot_recs_all, address_rhogs_folder, min_rhog
 def find_rhog_candidate_pairs(hogmaps, rhogs_prots):
     threshod_f_score_merging = 70
     pair_rhogs_count = {}
-    for species, prt_prot_maps in hogmaps.items():
+    for rhog, prt_prot_maps in hogmaps.items():
         for prot, prot_maps in prt_prot_maps.items():
             # [('HOG:D0017631.5a', '1771.7328874713658', '253', '234'), ('HOG:D0863448', '163.60700392437903', '253', '244'),
             rhogids = []
@@ -359,8 +360,8 @@ def find_rhog_candidate_pairs(hogmaps, rhogs_prots):
             ratioMax = count_shared / max(rhogs_size[hogi], rhogs_size[hogj])
             ratioMin = count_shared / min(rhogs_size[hogi], rhogs_size[hogj])
 
-            if (ratioMax > 0.8 or ratioMin > 0.9) and count_shared > 20 and rhogs_size[hogi] < _config.big_rhog_size / 2 and \
-                    rhogs_size[hogj] < _config.big_rhog_size / 2:
+            if (ratioMax > _config.mergHOG_ratioMax_thresh or ratioMin >  _config.mergHOG_ratioMin_thresh ) and _config.mergHOG_shared_thresh > 20 and\
+                    rhogs_size[hogi] < _config.big_rhog_size / 2 and  rhogs_size[hogj] < _config.big_rhog_size / 2:
                 if rhogs_size[hogi] >= rhogs_size[hogj]:
                     candidates_pair.append((hogi, hogj))  # bigger first
                 else:
@@ -428,19 +429,32 @@ def merge_rhogs(hogmaps, rhogs_prots):
     logger_hog.debug("There are " + str(len(rhogs_prots)) + " rhogs before merging.")
     print(len(rhogs_prots))
 
+    file_out_merge =  open("merging_rhogs.txt","w")
+    file_out_merge.write("#first column is the host hog, the rest will be merged here.\n")
     for cluster in cluster_rhogs_list:
-
-        prots = [rhogs_prots[hog] for hog in cluster]
-
+        for cluster_i in cluster:
+            file_out_merge.write(cluster_i+"\t")
+        file_out_merge.write("\n")
+    file_out_merge.close()
+    counter_merged_prots=0
+    for cluster in cluster_rhogs_list:
+        #prots = [rhogs_prots[hog] for hog in cluster]
         host_hog = cluster[0]
+        rhogs_host_size= len(rhogs_prots[host_hog])
         all_prots = []
         for hog in cluster:
             all_prots += rhogs_prots[hog]
             del rhogs_prots[hog]
-        rhogs_prots[host_hog] = all_prots
-
-    print(len(rhogs_prots))
-    logger_hog.debug("There are " + str(len(rhogs_prots)) + " rhogs after merging.")
+        all_prots_uniq = list(set(all_prots))  # there might be repeated prots
+        rhogs_prots[host_hog] = all_prots_uniq
+        # merging D0562038 and D0559070
+        # tr | C3ZG56 | C3ZG56_BRAFL HOG: D0562038, tr | H2Y1V7 | H2Y1V7_CIOIN HOG: D0562038, tr | C3ZG56 | C3ZG56_BRAFL HOG: D0559070, tr | H2Y1V7 | H2Y1V7_CIOIN HOG: D0559070
+        if len(all_prots_uniq) > rhogs_host_size:
+            counter_merged_prots += len(all_prots_uniq)
+        # otherwise, merging didn't help
+
+    print(len(rhogs_prots),counter_merged_prots)
+    logger_hog.debug("There are " + str(len(rhogs_prots)) + " rhogs by merging "+str(counter_merged_prots)+" proteins in total.")
 
     return rhogs_prots
 

FastOMA/_wrappers.py

@@ -137,7 +137,7 @@ def mad_rooting(input_tree_file_path: str, mad_executable_path: str = "./mad"):
     #os.system(mad_command)
 
     import subprocess
-
+    logger_hog.debug("MAD rooting started")
     #bashCommand = f"mad {gene_tree_address}"
     process = subprocess.Popen(mad_command.split(), stdout=subprocess.PIPE)
     output, error = process.communicate()
@@ -146,8 +146,21 @@ def mad_rooting(input_tree_file_path: str, mad_executable_path: str = "./mad"):
     #    print("error:\n", error)
 
     if "Error analyzing file" in str(output) or error:
-        rooted_tree = Tree(input_tree_file_path)
-    else:
+        try:
+            rooted_tree = Tree(input_tree_file_path)
+        except:
+            rooted_tree = Tree(input_tree_file_path,format=1)
+
+
+    elif "rooted trees written" in str(output): # 3 rooted trees written to tree_664187_Eukaryota.nwk.rooted Warning: Trees with repeating branch lengths are suspicious (3 repeating values).
+        file_multiple_tree_handle= open(input_tree_file_path + ".rooted",'r')
+        trees=[]
+        for line_tree1 in file_multiple_tree_handle:
+            if line_tree1.strip():
+                trees.append(line_tree1.strip())
+        # todo which one to choose ?
+        rooted_tree = Tree(trees[0])
+    else:#Rooted tree written to 'tree_672375_Theria.nwk.rooted'
         rooted_tree = Tree(input_tree_file_path + ".rooted")
-
+    logger_hog.debug("MAD rooting finished")
     return rooted_tree

FastOMA/infer_roothogs.py

@@ -48,6 +48,7 @@ def infer_roothogs():
 
 
     rhogs_prots = _utils_roothog.group_prots_roothogs(hogmaps)
+    # todo
     rhogs_prots = _utils_roothog.handle_singleton(rhogs_prots,hogmaps)
     rhogs_prots = _utils_roothog.merge_rhogs(hogmaps, rhogs_prots)
     rhogs_prots = _utils_roothog.roothogs_postprocess(hogmaps, rhogs_prots)

archive/test_curn.py

@@ -9,53 +9,53 @@
 # --input-rhog-folder ./bb/ --parrallel True  --species-tree species_tree.nwk
 
 #a=2
-#infer_subhogs()
+infer_subhogs()
 #infer_roothogs()
 
-
-from FastOMA.zoo.hog import transform
-
-#from zoo.tree_utils import collapse, gene_species, transform, HOG_coverages
-
-import io
-import lxml.etree
-orthoxml_file = "/work/FAC/FBM/DBC/cdessim2/default/smajidi1/gethog3_qfo/benchmark-webservice3/orthoxml/euk_omamer200.dev8_13oct.orthoxml"
-
-
-orthxml_str = []
-with open(orthoxml_file, "r") as f:
-    for i in f:
-        orthxml_str.append(i)
-print(len(orthxml_str))
-dic_gene_integer={}
-for line in orthxml_str:
-    if "gene id" in line:
-        found=False
-        gene_int= line.split("\"")[1]
-        gene_name = line.split("\"")[3]
-        dic_gene_integer[gene_int] = gene_name
-
-
-
-orthoxml_etree=lxml.etree.parse(orthoxml_file)
-
-pw_orthologs_integer = sorted(list(transform.iter_pairwise_relations(orthoxml_etree)))
-# iter_pairwise_relations(obj, rel_type=None    (def:'ortholog' , but possible to use 'paralog')
-print(len(pw_orthologs_integer))
-print(pw_orthologs_integer[:2])
-pw_orthologs_gene =[]
-for pair in pw_orthologs_integer:
-    pw_orthologs_gene.append((dic_gene_integer[pair[0]],dic_gene_integer[pair[1]]))
-
-
-
-print(len(pw_orthologs_gene))
-
-output_file = open(orthoxml_file+"_pairs.tsv","w")
-for  pair in pw_orthologs_gene:
-    output_file.write(pair[0]+"\t"+pair[1]+"\n")
-
-output_file.close()
+#
+# from FastOMA.zoo.hog import transform
+#
+# #from zoo.tree_utils import collapse, gene_species, transform, HOG_coverages
+#
+# import io
+# import lxml.etree
+# orthoxml_file = "/work/FAC/FBM/DBC/cdessim2/default/smajidi1/gethog3_qfo/benchmark-webservice3/orthoxml/euk_omamer200.dev8_13oct.orthoxml"
+#
+#
+# orthxml_str = []
+# with open(orthoxml_file, "r") as f:
+#     for i in f:
+#         orthxml_str.append(i)
+# print(len(orthxml_str))
+# dic_gene_integer={}
+# for line in orthxml_str:
+#     if "gene id" in line:
+#         found=False
+#         gene_int= line.split("\"")[1]
+#         gene_name = line.split("\"")[3]
+#         dic_gene_integer[gene_int] = gene_name
+#
+#
+#
+# orthoxml_etree=lxml.etree.parse(orthoxml_file)
+#
+# pw_orthologs_integer = sorted(list(transform.iter_pairwise_relations(orthoxml_etree)))
+# # iter_pairwise_relations(obj, rel_type=None    (def:'ortholog' , but possible to use 'paralog')
+# print(len(pw_orthologs_integer))
+# print(pw_orthologs_integer[:2])
+# pw_orthologs_gene =[]
+# for pair in pw_orthologs_integer:
+#     pw_orthologs_gene.append((dic_gene_integer[pair[0]],dic_gene_integer[pair[1]]))
+#
+#
+#
+# print(len(pw_orthologs_gene))
+#
+# output_file = open(orthoxml_file+"_pairs.tsv","w")
+# for  pair in pw_orthologs_gene:
+#     output_file.write(pair[0]+"\t"+pair[1]+"\n")
+#
+# output_file.close()
 
 
 #