Fix error caused by exessive trimming by trimal (#66)

src/fastdtlmapper/scripts/FastDTLmapper.py

@@ -146,13 +146,19 @@ def trimal_run(outpath: OutPath, cmd: Cmd) -> None:
         trimal_cmd_list, outpath.tmp_parallel_cmds_file, outpath.trimal_log_file
     )
 
-    # If trimal removed sequence, use raw mafft alignment in next step
+    # If trimal removed sequence or exessive trimming
+    # use raw mafft alignment in next step
     for aln_file in sorted(outpath.dtl_rec_dir.glob("**/*_aln.fa")):
         group_id = aln_file.parent.name
         aln_trim_file = aln_file.parent / (group_id + "_aln_trim.fa")
         if UtilFasta(aln_file).seq_num > UtilFasta(aln_trim_file).seq_num:
             shutil.copy(aln_file, aln_trim_file)
 
+        aln_length = UtilFasta.get_aln_length(aln_file)
+        aln_trim_length = UtilFasta.get_aln_length(aln_trim_file)
+        if aln_trim_length < 5 or aln_length > aln_trim_length * 10:
+            shutil.copy(aln_file, aln_trim_file)
+
 
 @print_runtime
 def iqtree_run(outpath: OutPath, cmd: Cmd) -> None:

src/fastdtlmapper/util/fasta.py

@@ -3,7 +3,7 @@
 from pathlib import Path
 from typing import Dict, List, Union
 
-from Bio import SeqIO
+from Bio import AlignIO, SeqIO
 
 
 @dataclass
@@ -72,3 +72,15 @@ def add_serial_id(
             record.id = f"{id_prefix}_GENE{cnt:06d}_{record_id}"
             fix_records.append(record)
         SeqIO.write(fix_records, fasta_outfile, "fasta-2line")
+
+    @classmethod
+    def get_aln_length(cls, aln_file: Union[str, Path]) -> int:
+        """Get alignment length
+
+        Args:
+            aln_file (Union[str, Path]): Alignment file
+
+        Returns:
+            int: alignment length
+        """
+        return AlignIO.read(aln_file, "fasta").get_alignment_length()