update

Author: sinamajidian

FastOMA.nf

@@ -154,6 +154,7 @@ workflow {
     (pickle_rest_rhog,  msas_out_rest, genetrees_out_test, ready_hog_rest) = hog_rest(rhogsrest_tree_ready)
 
     (orthoxml_file, OrthologousGroupsFasta, OrthologousGroups_tsv, rootHOGs_tsv)  = collect_subhogs(ready_hog_rest.collect(), ready_hog_big.collect(), pickles_temp, gene_id_dic_xml, omamer_rhogs)
+    omamer_rhogs.view{" output omamer_rhogs ${it}"}
     orthoxml_file.view{" output orthoxml file ${it}"}
 
 }

FastOMA/_config.py

@@ -40,7 +40,7 @@
 keep_subhog_each_pickle = False
 
 big_rhog_size = 60 * 1000
-omamer_family_threshold = 110
+omamer_family_threshold = 90
 #
 # omamer_fscore_treshold_big_rhog = 0.04 # 0.5 # means no thresold #0.2 #0.5  #  to have more proteins in the ortho groups 0.05  considering for big rhogs
 # omamer_treshold_big_rhog_szie = 100 #9000 #100
@@ -49,7 +49,7 @@
 # omamer_treshold_big_rhog_szie2 = 50*1000
 # omamer_fscore_treshold_big_rhog2 = 0.6 #0.9
 
-hogclass_max_num_seq = 20  # subsampling in msa # ver very 2
+hogclass_max_num_seq = 40  # subsampling in msa # ver very 2
 hogclass_min_cols_msa_to_filter = hogclass_max_num_seq * 50
 hogclass_tresh_ratio_gap_col = 0.6  # 0.8 for very very big
 # old code after samplign if there are 2 seq sampled, then at least one nongap
@@ -70,7 +70,7 @@
 rooting_mad_executable_path = "mad"  # /work/FAC/FBM/DBC/cdessim2/default/smajidi1/software/installers/mad/
 
 ##inferhog
-inferhog_tresh_ratio_gap_row =0.4 # 0.6   # to have more proteins in the ortho groups 0.1
+inferhog_tresh_ratio_gap_row =0.1 # 0.6   # to have more proteins in the ortho groups 0.1
 inferhog_tresh_ratio_gap_col =0.5  # 0.6   # ver very 0.8
 inferhog_min_cols_msa_to_filter = 50 #300 #50  # used for msa before gene tree inference and  saving msa in hog class
 

README.md

@@ -109,7 +109,7 @@ Then, cd to the `testdata` folder and download the omamer database and change it
 ```
 cd FastOMA/testdata
 wget https://omabrowser.org/All/Primates-v2.0.0.h5     # 105MB
-mv Primates.h5  in_folder/omamerdb.h5 
+mv Primates-v2.0.0.h5    in_folder/omamerdb.h5 
 ```
 (This is for the test however, I would suggest downloading the `LUCA-v2.0.0.h5` instead of `Primates-v2.0.0.h5` for your real analysis.). Check the item 2 in the [input section](https://github.com/sinamajidian/FastOMA#input) for details.
 

archive/test_curn.py

@@ -8,8 +8,82 @@
 
 # --input-rhog-folder ./bb/ --parrallel True  --species-tree species_tree.nwk
 
-a=2
-infer_subhogs()
+#a=2
+#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()
+
+
+#
+#
+# # orthoxml_handle= open(orthoxml_file,"r")
+# # orthoxml =""
+# # for line in orthoxml_handle:
+# #     orthoxml+=line
+#
+#
+# from xml.etree.ElementTree import XMLParser
+#
+# parser = XMLParser()
+# with open(orthoxml_file, 'rb') as xml:
+#     for chunk in xml:
+#         parser.feed(chunk)
+# parser.close()
+#
+#
+# lxml.etree.parse(oxml)
+#
+# orthoxm= lxml.etree.parse(orthoxml)
+#
+# # expected = [("1", "2"), ("1", "3"), ("1", "4"), ("1", "5"), ("1", "6"),
+# #             ("2", "5"), ("2", "6"), ("3", "4"), ("3", "5"), ("3", "6"),
+# #             ("4", "5"), ("4", "6"), ("5", "6")]
+# #    self.assertEqual(expected, pw_orthologs)
+#
+# from xml.etree import ElementTree
+# tree = ElementTree.parse(orthoxml_file)
+# root = tree.getroot()

setup.py

@@ -10,7 +10,7 @@
             exec(line.rstrip())
 
 # TODO
-requirements = ['biopython', 'ete3', 'omamer>=2.0.0.dev8', 'nextflow', 'pyparsing' , 'DendroPy', 'future', 'lxml','pyham']
+requirements = ['biopython', 'ete3', 'omamer>=2.0.0', 'nextflow', 'pyparsing' , 'DendroPy', 'future', 'lxml','pyham']
 
 desc = 'FastOM - a package to infer orthology information '