Merge branch 'dev' into 'main'

bug correction: considers aligned ref_seq

See merge request !2

influenza_resistance_annotation.py

@@ -211,7 +211,7 @@ def align_pairwise(ref_seq, target_seq,path_to_alignment_tool):
     os.remove(output_fasta)
     return alignment      
 
-def check_presence_mutations(protein_seq, mutations):
+def check_presence_mutations(reference_seq, protein_seq, mutations):
     mutations_present = []
     for mut in mutations:
         if "+" in mut:
@@ -221,7 +221,7 @@ def check_presence_mutations(protein_seq, mutations):
                 ref_aa = mutis[0]
                 pos = mutis[1:-1]
                 alt_aa = mutis[-1]
-                if protein_seq[int(pos)-1] == alt_aa:
+                if (protein_seq[int(pos)-1] == alt_aa) & (reference_seq[int(pos)-1] == ref_aa):
                     mutation_present+=1
 
             if mutation_present == len(muts): # meaning that all individual mutations are present
@@ -239,9 +239,14 @@ def check_presence_mutations(protein_seq, mutations):
             pos = mut[1:-1]
             alt_aa = mut[-1]
 
-            if protein_seq[int(pos)-1] == alt_aa:
+            if (protein_seq[int(pos)-1] == alt_aa) & (reference_seq[int(pos)-1] == ref_aa):
                 mutations_present.append(mut)
 
+    # if a mutation is part of a combination, remove the individual call, e.g. E119V, E119V+T148I should be only E119V+T148I
+    for mut in mutations_present:
+        if any(mut in item and mut != item for item in mutations_present):
+            mutations_present.remove(mut)
+
     return mutations_present
 
 def main():
@@ -307,17 +312,18 @@ def main():
     for protein_seq in protein_seqs:
         # align sequences
         alignment = align_pairwise(reference_seq,protein_seq,path_to_alignment_tool)
+        reference_seq = alignment[0]
         protein_seq_aligned = alignment[1]
         
         print("\n>Aligned ref sequence:")
-        print(alignment[0])
+        print(reference_seq)
 
         print("\n>Aligned sequence:")
         print(protein_seq_aligned)
 
         # Test presence of resistance mutations in the translated sequence
         print(f"\n>Identified mutations in {args.fasta_file} {args.protein}:")
-        mutations_present = check_presence_mutations(protein_seq_aligned,mutations)
+        mutations_present = check_presence_mutations(reference_seq,protein_seq_aligned,mutations)
         all_mutations_present.append(mutations_present)
 
     all_mutations_present = [item for sublist in all_mutations_present for item in sublist] # flatten list