Solved issue #201

src/pdbcleanresiduestandardizationutils.py

@@ -3,6 +3,8 @@
 from Bio.PDB.MMCIFParser import FastMMCIFParser
 from PDBClean.alignmentutils import *
 from PDBClean.listutils import *
+import string
+
 #
 
 ####################
@@ -104,9 +106,36 @@ def perform_multiple_alignment(Structure_Sequences, ChID_ResiNum_Vector, structi
     check : str
         Updated string representing the state of the main menu, set to '1' to indicate a state change.
     """
+
+    def get_gap_letter(n):
+        """
+            Convert a 0-based index into a spreadsheet-style alphanumeric label.
+
+            Cycles through all single uppercase letters (A-Z), lowercase letters (a-z), and digits (0-9)
+            before moving into multi-letter combinations (AA, AB, ..., 99, AAA...).
+
+            Parameters:
+                n (int): The 0-based index to convert.
+
+            Returns:
+                str: The corresponding label.
+            """
+        alphabet = string.ascii_uppercase + string.ascii_lowercase + string.digits
+        base = len(alphabet)
+
+        label = ""
+        n += 1
+        while n > 0:
+            n, remainder = divmod(n - 1, base)
+            label = alphabet[remainder] + label
+
+        return label
+
     Structure_Sequences_Aligned = {}
     Structure_ConversionTemplate = {}
     Structure_Sequences_GAPS = {}
+    big_structures_warning = {}
+
     input_submenu = ""
     while(input_submenu != "QUIT"):
         print("    Perform multiple alignments to identify residues",
@@ -145,8 +174,10 @@ def perform_multiple_alignment(Structure_Sequences, ChID_ResiNum_Vector, structi
                 i = 0
                 for I in range(len(structid_list)):
                     key = str(structid_list[I]) + "_" + chid
+
                     if key in Structure_Sequences:
                         #Structure_Sequences_Aligned[key] = this_chainsseq_aligned_list[i]
+
                         Structure_Sequences_Aligned[key] = this_chainsseq_aligned_list_map[str(structid_list[I])]
                         Structure_Sequences_GAPS[key] = this_chainseq_gap_percentages #FAPA
                         i += 1
@@ -189,12 +220,13 @@ def perform_multiple_alignment(Structure_Sequences, ChID_ResiNum_Vector, structi
                         new_res_num=[]
                         freq_tracker=1
                         gap_tracker=0
-                        gap_letter=['A','B','C','D','E','F','G','H','I','J','K','L','M','N','O','P','Q','R','S','T','U','V','W','X','Y','Z','a','b','c','d','e','f','g','h','i','j','k','l','m','n','o','p','q','r','s','t','u','v','w','x','y','z']
-                        for freq in gaps:
+                        for freq in gaps: # accepted gap percentage based on user defined occupancy threshold
+
+                        # Changing the functionality of gap_letter
                             #freq_tracker = 1
                             #gap_tracker = 0
                             #print(freq)
-                            if freq < 100-user_gap: # accepted gap percentage based on user defined occupancy threshold
+                            if freq < 100-user_gap:
                                 new_res_num.append(counter)
                                 counter += 1
                                 freq_tracker=freq
@@ -203,7 +235,11 @@ def perform_multiple_alignment(Structure_Sequences, ChID_ResiNum_Vector, structi
                             else:
                                 #print(gap_tracker)
                                 #new_res_num.append(str(counter)+"_"+str(gap_tracker))
-                                new_res_num.append(str(counter-1) +" "+str(gap_letter[gap_tracker]))
+                                new_res_num.append(str(counter-1) +" "+str(get_gap_letter(gap_tracker)))
+
+                                if gap_tracker == 62:
+                                    structure = (structid_list[I].split("/")[-1])
+                                    big_structures_warning.setdefault(structure, []).append(chain)
                                 gap_tracker+=1
                                 freq_tracker=freq
 
@@ -266,6 +302,19 @@ def perform_multiple_alignment(Structure_Sequences, ChID_ResiNum_Vector, structi
 
             check = "1"
             input_submenu = "QUIT"
+
+    if len(big_structures_warning) > 0:
+        print(f"\nWARNING: Some structures were affected by large gaps in most of the MSA (>62 positions), "
+              f"resulting in residue renumbering that requires two-character alternate identifiers. This may lead "
+              f"to problems when opening the structure with molecular visualization software. These structures "
+              f"will be saved in residue_warning_structures.csv in the format 'structure:chain'\n")
+
+        with open("residue_warning_structures.csv", "w") as f:
+            f.write("structure:chain\n")
+            for structure in big_structures_warning:
+                for chain in big_structures_warning[structure]:
+                    f.write(f"{structure}:{chain}\n")
+
     return Structure_Sequences_Aligned, Structure_ConversionTemplate, chid_list, check
 
 def show_conversiontemplate(Structure_ConversionTemplate):
@@ -489,4 +538,4 @@ def conversiontemplate_to_pdb_FAPA(filelist, Structure_ConversionTemplate, targe
                     else:
                         newciffile.write(line)
 
-# FAPA MAY 2024 TEST ENDS
+# FAPA MAY 2024 TEST ENDS