checkpoint_1Nov

Author: sinamajidian

FastOMA/_config.py

@@ -36,11 +36,11 @@
 # threshold_dubious_sd = float(os.getenv(VARIABLE_threshold_dubious_sd, 0.1))
 
 ## output writing files
-gene_trees_write = True
-msa_write = True
-gene_trees_write_all = True
-msa_write_all = True
-keep_subhog_each_pickle = True
+gene_trees_write = False
+msa_write = False
+gene_trees_write_all = False
+msa_write_all = False
+keep_subhog_each_pickle = False
 
 big_rhog_size = 60 * 1000
 omamer_family_threshold = 90

FastOMA/_hog_class.py

@@ -201,6 +201,7 @@ def _sorter_key(sh):
         element_list = []
         for sub_clade, sub_hogs in itertools.groupby(self._subhogs, key=_sorter_key):  # sub_clade is the taxrange
             list_of_subhogs_of_same_clade = list(sub_hogs)
+            # list_of_subhogs_of_same_clade is  [object of HOGclass hogID=HOG_D0634402_sub10094,len=12, tax_least=Amniota, tax_now= Amniota, object of HOGclass hogID=HOG_D0634402_sub10096,len=20, tax_least=Amniota, tax_now= Amniota, object of HOGclass hogID=HOG_D0634402_sub10081,len=1, tax_least=MONDO, tax_now= Amniota, object of HOGclass hogID=HOG_D0634402_sub10093,len=12, tax_least=Amniota, tax_now= Amniota, object of HOGclass hogID=HOG_D0634402_sub10095,len=7, tax_least=Amniota, tax_now= Amniota]
             # following only for debugging, can be deleted later
             for subhog in list_of_subhogs_of_same_clade:
                 if len(subhog._members) == 0:

FastOMA/_utils_frag_SO_detection.py

@@ -191,14 +191,20 @@ def find_prot_dubious_sd_remove(gene_tree, all_species_dubious_sd_dic):
                                 prot_dubious_list.append(prot_name)
                     try:
                         subhogs_list = [i.split("|_|")[1] for i in prot_dubious_list]  # subhog id at child level
-                        if len(set(subhogs_list)) > 1:
-                            # we are removing all sequences of this species on the the side of internal node (gene tree), with least leaves
-                            child_size_min_indx = child_size.index(min(child_size))
-                            prot_dubious_sd_remove_list.append(prot_dubious_list[child_size_min_indx])
-
-                        else:
-                            logger_hog.debug( "This species (protein from the same subhog) is safe to keep "+ str(node_name)+" "+str(species_dubious_sd))
-                            #all of them are from the same subhog, so it doesn't matter, a duplication event doesn't affect when all are from the same subhog at children level
+
+                        # todo check !, is it safe or not!
+                        # if len(set(subhogs_list)) > 1:
+                        #     # we are removing all sequences of this species on the the side of internal node (gene tree), with least leaves
+                        #     child_size_min_indx = child_size.index(min(child_size))
+                        #     prot_dubious_sd_remove_list.append(prot_dubious_list[child_size_min_indx])
+                        #
+                        # else:
+                        #     logger_hog.debug( "This species (protein from the same subhog) is safe to keep "+ str(node_name)+" "+str(species_dubious_sd))
+                        #     #all of them are from the same subhog, so it doesn't matter, a duplication event doesn't affect when all are from the same subhog at children level
+
+                        child_size_min_indx = child_size.index(min(child_size))
+                        prot_dubious_sd_remove_list.append(prot_dubious_list[child_size_min_indx])
+
                     except:
                         logger_hog.warning("issue 2495869: prot_dubious_list doesnt include the hog id . so we'll  keep it" + str(gene_tree) + " " + str(prot_dubious_list))
 

FastOMA/_utils_roothog.py

@@ -35,7 +35,7 @@ def parse_proteomes(folder="./"):  # list_oma_species
         prot_recs_lists[species_name]=prots_record
 
     logger_hog.info("The are "+str(len(species_names))+" species in the proteome folder.")
-    return species_names, prot_recs_lists
+    return species_names, prot_recs_lists, fasta_format_keep
 
 
 
@@ -78,7 +78,7 @@ def add_species_name_prot_id(species_names, prot_recs_lists):
     return prot_recs_all
 
 
-def parse_hogmap_omamer(species_names, folder="./"):
+def parse_hogmap_omamer(species_names, fasta_format_keep,  folder="./"):
     """
      function for parsing output of omamer (hogmap files) located in /hogmap/
     Each hogmap file correspond to one fasta file of species, with the same name.
@@ -90,7 +90,7 @@ def parse_hogmap_omamer(species_names, folder="./"):
     """
     hogmaps = {}
     for species_name in species_names:
-        hogmap_address = folder + "/hogmap/" + species_name + ".fa.hogmap"
+        hogmap_address = folder + "/hogmap/" + species_name + "."+fasta_format_keep+".hogmap"
         hogmap_file = open(hogmap_address, 'r')
 
         for line in hogmap_file:
@@ -256,20 +256,23 @@ def roothogs_postprocess(hogmaps, rhogs_prots):
         sp_prot_list_filt = []
         hogids = []
 
-        # removing protines with low omamer score from big rootHOG
+        # removing proteins with low omamer score from big rootHOG
         for species_name, prot_id in sp_prot_list:
             prot_maps = hogmaps[species_name][prot_id]
-            if len(prot_maps) > 1:  # probably for big rootHOG, there won't be multi-hits
-                scores = [float(i[1]) for i in prot_maps]  # (hogid,score,seqlen,subfamily_medianseqlen)
-                hogids = [i[0] for i in prot_maps]
-                max_score = max(scores)
-                max_index = scores.index(max_score)
-                hogid = hogids[max_index]
-            else:
-                hogid = prot_maps[0][0]
-                score = float(prot_maps[0][1])
+            # if len(prot_maps) > 1:  # probably for big rootHOG, there won't be multi-hits
+            #     scores = [float(i[1]) for i in prot_maps]  # (hogid,score,seqlen,subfamily_medianseqlen)
+            #     hogids = [i[0] for i in prot_maps]
+            #     max_score = max(scores)
+            #     max_index = scores.index(max_score)
+            #     hogid = hogids[max_index]
+            # else:
+            #     hogid = prot_maps[0][0]
+            #     max_score = float(prot_maps[0][1])
+            # todo: highest normcount or pavalue , default of omamer ?
+            hogid = prot_maps[0][0]
+            max_score = float(prot_maps[0][1])
 
-            if score > _config.omamer_family_threshold:
+            if max_score > _config.omamer_family_threshold:
                 sp_prot_list_filt.append((species_name, prot_id))
                 hogids.append(hogid)
 
@@ -279,21 +282,23 @@ def roothogs_postprocess(hogmaps, rhogs_prots):
         if len(sp_prot_list_filt) < _config.big_rhog_size:
             sp_prot_list_filt2 = sp_prot_list_filt
 
-
-        else:  # removing protines that are mapped to rootHOG (=HOGC123 not the levelsHOGC123.1a) from big rootHOG
-            hogids2 = []
+        else:  # removing proteins that are mapped to rootHOG (= HOGC123 not the levelsHOGC123.1a) from big rootHOG
+            #hogids2 = []
             sp_prot_list_filt2 = []
-            for prot_idx, sp_prot in enumerate(sp_prot_list_filt):
-                hogid = hogids[prot_idx]
-
-                if len(hogid.split(".")) > 1:
+            for prot_idx, sp_prot in enumerate(sp_prot_list_filt): #[('UP000192223_224129', 'tr|A0A1W4WAU6|A0A1W4WAU6_AGRPL'), ('UP000192223_224129', 'tr|A0A1W4WU99|A0A1W4WU99_AGRPL'),
+                species , protein = sp_prot
+                prot_maps =hogmaps[species][protein]
+                hogid = prot_maps[0][0]
+                if len(hogid.split(".")) > 1: # if mapped to subHOG (not to the rootHOG)
                     sp_prot_list_filt2.append(sp_prot)
-                    hogids2.append(hogid)
+                    #hogids2.append(hogid)
 
         logger_hog.info("For big rootHOG " + rhogid + ", " + str(
             len(sp_prot_list_filt2)) + " proteins left after removing non-child subhogs")
-
-        rhogs_prots[rhogid] = sp_prot_list_filt2
+        if len(sp_prot_list_filt2):
+            rhogs_prots[rhogid] = sp_prot_list_filt2
+        else:
+            del rhogs_prots[rhogid]
 
     return rhogs_prots
 

FastOMA/collect_subhogs.py

@@ -60,7 +60,7 @@ def collect_subhogs():
     pickle_files_adress = [i for i in pickle_files_adress_raw if i.endswith(".pickle") and i.startswith("file_")]
 
     logger_hog.info("number of pickle files is "+str(len(pickle_files_adress))+".")
-    logger_hog.debug("pickle files are " + str(pickle_files_adress) + ".")
+    logger_hog.debug("pickle files are " + str(len(pickle_files_adress)) + ".")
     hogs_a_rhog_xml_all = []
     for pickle_file_adress in pickle_files_adress:
         with open(pickle_folder + pickle_file_adress, 'rb') as handle:

FastOMA/infer_roothogs.py

@@ -33,10 +33,10 @@ def infer_roothogs():
     # sys.exit(0)
 
     #folder = "/scratch/smajidi1/euk_omamer200.dev8_2/test/hogmap/FastOMA-main/testdata/in_folder/"
-    species_names, prot_recs_lists = _utils_roothog.parse_proteomes()# folder
+    species_names, prot_recs_lists,fasta_format_keep = _utils_roothog.parse_proteomes() # optional input folder
     prot_recs_all = _utils_roothog.add_species_name_prot_id(species_names, prot_recs_lists)
 
-    hogmaps = _utils_roothog.parse_hogmap_omamer(species_names)#folder
+    hogmaps = _utils_roothog.parse_hogmap_omamer(species_names,fasta_format_keep) # optional input folder
 
     splice_files =  os.path.exists("./splice/")
     if splice_files:
@@ -48,7 +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,8 +9,8 @@
 # --input-rhog-folder ./bb/ --parrallel True  --species-tree species_tree.nwk
 
 #a=2
-infer_subhogs()
-#infer_roothogs()
+#infer_subhogs()
+infer_roothogs()
 
 #
 # from FastOMA.zoo.hog import transform