[FT] use full seqs at every level for alignments

Author: alpae

FastOMA.nf

@@ -372,6 +372,7 @@ process hog_big{
                                --gap-ratio-row ${params.filter_gap_ratio_row} \
                                --gap-ratio-col ${params.filter_gap_ratio_col} \
                                --number-of-samples-per-hog ${params.nr_repr_per_hog} \
+                               --align-subhogs-seqs full-seqs \
                                ${ params.write_msas ? "--msa-write" : ""} \
                                ${ params.write_genetrees ? "--gene-trees-write" : ""}
     """
@@ -404,6 +405,7 @@ process hog_rest{
                               --gap-ratio-row ${params.filter_gap_ratio_row} \
                               --gap-ratio-col ${params.filter_gap_ratio_col} \
                               --number-of-samples-per-hog ${params.nr_repr_per_hog} \
+                              --align-subhogs-seqs full-seqs \
                               ${ params.write_msas ? "--msa-write" : ""} \
                               ${ params.write_genetrees ? "--gene-trees-write" : ""}
     """

FastOMA/_infer_subhog.py

@@ -387,8 +387,21 @@ def write_msa_or_tree_if_necessary(self, elem, fn_suffix="", features=[]):
                         for member in rep_val.representative.get_subelements():
                             writer.writerow([n.name, member])
 
+    def get_sequences_from_subhogs(self):
+        if self.conf.align_subhogs_seqs == "msa":
+            return [hog.get_msa() for hog in self.subhogs.values() if len(hog.get_msa()) > 0]
+        elif self.conf.align_subhogs_seqs == "full-seqs":
+            seqs = []
+            for hog in self.subhogs.values():
+                for rep in hog.get_representatives():
+                    seqs.append(MultipleSeqAlignment([rep.get_record()]))
+            assert len(seqs) == sum(len(subhog.get_msa()) for subhog in self.subhogs.values())
+            return seqs
+        else:
+            raise ValueError("Unknown align_subhogs_seqs value")
+
     def align_subhogs(self):
-        sub_msas = [hog.get_msa() for hog in self.subhogs.values() if len(hog.get_msa()) > 0]
+        sub_msas = self.get_sequences_from_subhogs()
         logger.debug(f"Merging {len(sub_msas)} MSAs for rhog: {self.rhogid}, level: {self.node_species_tree.name}")
         if len(sub_msas) == 0:
             logger.info(

FastOMA/_wrappers.py

@@ -41,7 +41,7 @@ def merge_msa(list_msas):
     #     wrapper_mafft_merge.options['--merge'].active = False
 
     #mafft --auto Automatically selects an appropriate strategy from L-INS-i, FFT-NS-i and FFT-NS-2, according to data size. Default: off (always FFT-NS-2)
-    wrapper_mafft_merge.options['--auto'].set_value(False) # todo we can make it as an argument in fastoma-infer-subhogs. 
+    wrapper_mafft_merge.options['--auto'].set_value(False)    # todo we can make it as an argument in fastoma-infer-subhogs.
     # wrapper_mafft_merge.options['--anysymbol'].active = True
     wrapper_mafft_merge.options['--anysymbol'].set_value(True)
     wrapper_mafft_merge.options['--thread'].set_value(-1) # -1 uses a largely appropriate number of threads in each step, after automatically counting the number of physical cores the computer has.

FastOMA/infer_subhogs.py

@@ -49,6 +49,10 @@ def fastoma_infer_subhogs():
                         help="For trimming the MSA, the threshold of ratio of gaps for each column.")
     parser.add_argument("--min-col-trim", required=False, type=int, default=50,  # todo min rows trim
                         help="min no. columns in msa to consider for filtering")
+    parser.add_argument("--align-subhogs-seqs", default="msa", choices=("msa", "full-seqs"),
+                        help="Sequences used in ancestral nodes to build MSA from. 'msa' uses the trimmed, "
+                             "aligned sequences from the more recent taxonomic level, where as 'full-seqs' "
+                             "uses always the initial full protein sequences.")
     parser.add_argument('-v', action="count", default=0, help="Increase verbosity to info/debug")
     conf_infer_subhhogs = parser.parse_args()