2
\$\begingroup\$

A protein is composed of amino acids (also called residues). The amide nitrogen and hydrogens (N and H in the script), carbonyl carbon (C), alpha carbon (Ca), beta carbon (Cb), and alpha hydrogen (HA) form the backbone. Each amino acid has the N,H,C,CA,HA atoms, with most amino acids containing the Cb as well. Each amino acid will have a particular value (chemical shift) determined from experiments. However the experiments just give you raw values for each atom type, they do not tell you which value belongs to which amino acid (as a scientist, it is my job to determine which value fits to what amino acid) There is a program that can predict these values, for each amino acid (SPARTA). I have created a program to calculate the RMSD from the experimental value to the predicted value for each amino acid.

Both the experimental values, and predicted values, have a particular format (NMRSTAR and SPARTA). I have decided to convert the format of each file so that each amino acid has 6 atom types in both files (I use placeholders with values of 1000 if that atom type is not there, makes it easy to ignore when calculating RMSD), and filter both files to one another so they are the same size. This makes it much easier to calculate RMSDs between the 2.

To ensure ease of use, I have created a GUI. Initially I had the GUI script, and the functions that did all the convertions for both files, in the same script. This turned out to be a pain to troubleshoot, and difficult to read. Thus, I have split the GUI and convertions into separate files, import these separate files into the GUI script.

This is the first time that I've tried using functions and splitting into separate files and importing (I've always done everything in one long script, no functions or imports). As well as first time adding comments. Thus any feedback on structure and function use would also be greatly appreciated!

#The GUI Script (only the parts relevant to the code)
def nmrstarrun3():
    text_area.delete(1.0,END)
    #user inputs
    if sparta_file == ():
        text_area.insert(tk.INSERT,'please upload your sparta file (make sure to use browse)\n')
    if seq_file == ():
        text_area.insert(tk.INSERT,'please upload your seq file (make sure to use browse)\n')
    if save_file_sparta == ():
        text_area.insert(tk.INSERT,'please indicate sparta save file (make sure to use browse)\n')
    if save_file_peaklist == ():
        text_area.insert(tk.INSERT,'please indicate peaklist save file (make sure to use browse)\n')
    if set_threshold == ():
        text_area.insert(tk.INSERT,'please enter a threshold (make sure to hit enter)\n')
    if seq_start == ():
        text_area.insert(tk.INSERT,'please enter a seq number (make sure to hit enter)\n')
    if nmrstarfile == ():
        text_area.insert(tk.INSERT,'please upload your nmrstar file (make sure to use browse)\n')
    else:
        text_area.insert(tk.INSERT,'Starting Program\n')
        text_area.insert(tk.INSERT,'Creating Sparta File\n')
        text_area.update_idletasks()

        acid_map = {
                  'ASP':'D', 'THR':'T', 'SER':'S', 'GLU':'E',
                  'PRO':'P', 'GLY':'G', 'ALA':'A', 'CYS':'C',
                  'VAL':'V', 'MET':'M', 'ILE':'I', 'LEU':'L',
                  'TYR':'Y', 'PHE':'F', 'HIS':'H', 'LYS':'K',
                  'ARG':'R', 'TRP':'W', 'GLN':'Q', 'ASN':'N'
                }

        os.chdir(nmrstarfile_directory)
        #NMRSTAR files contain a variety of information, and side chain chemical shift values
        #We only want residues with backbone N,HA,C,CA,CB,H chemical shifts
        #Additionally, NMRSTAR file amino acids numbers are not always correct (they contain additional values). Thus the user defines what the starting value should be
        #NMRSTAR uses 3 letter amino acid abbreviations, we want single-letter, the acid map is used to convert
        exctracted_and_compiled_data=[]
        with open(nmrstarfile) as file:
            for lines in file:
                modifier=lines.strip()
                extract_data_only=re.search(r'\b\d+\s+[A-Z]{3}\s+\w+\s+\w+\s+\d+\s+\d+',modifier)
                if extract_data_only != None:
                    atom_search=extract_data_only.string
                    split_data=atom_search.split()
                    amino_acid_number=str(int(split_data[5])+int(seq_start)-1)
                    residue_type=split_data[6]
                    atom_type=split_data[7]
                    converted=acid_map[residue_type]
                    chemical_shift=split_data[10]
                    compile_data=[amino_acid_number]+[converted]+[atom_type]+[chemical_shift]
                    if atom_type == 'N' or atom_type == 'HA' or atom_type =='CA' or atom_type == 'CB' or atom_type=='H' or atom_type=='C':
                        joined=' '.join(compile_data)
                        exctracted_and_compiled_data.append(joined)
        from sparta_file_formatter import check_sparta_file_boundaries
        from nmrstar import dict_create
        from nmrstar import fill_missing_data
        dict_create(seq_file,seq_start,seq_directory)
        sparta_file_boundaries=check_sparta_file_boundaries(seq_file,seq_directory,mutation_list1,mutation_list2,sparta_file,sparta_directory,seq_start)
        data_files=fill_missing_data(final_list,seq_start)
        #The peaklist may have additional chemical shifts not present in the crystal structure, and thus sparta file
        #We filter out and create a new list containing only the residues found in the sparta file
        peaklist_filtered_to_match_sparta=[]
        count=0
        for lines in data_files:
            modify=lines.strip()
            splitting=modify.split()
            number_search=re.search('^-*\d+[A-Z]',splitting[0])
            r=re.compile(number_search.group(0))
            comparison_to_sparta=list(filter(r.match,sparta_file_boundaries))
            if comparison_to_sparta != []:
                peaklist_filtered_to_match_sparta.append(modify)
            else:
                count+=1
                if count==6:
                    #if any amino acid is the peaklist, but not SPARTA file, it will be excluded and printed out here
                    count=0
                    text_area.insert(tk.INSERT,f'{splitting[0]} was excluded\n')
        #RMSD values are calculated summing the deviations of the experimental with predicted values, and dividing it by the number of atoms used in the calculation
        amino_acid_square_deviation_values=[]
        number=0
        for experimental,predictions in zip(peaklist_filtered_to_match_sparta,sparta_file_boundaries):
            number+=1
            experimental_split=experimental.split()
            predictions_split=predictions.split()
            square_deviation=((float(predictions_split[1])-float(experimental_split[1]))**2)/((float(predictions_split[2]))**2)
            if square_deviation>100:
                square_deviation=0
            else:
                amino_acid_square_deviation_values.append(square_deviation)
            if number%6 ==0:
                if len(amino_acid_square_deviation_values)==0:
                    continue
                else:
                    rmsd=math.sqrt((1/int(len(amino_acid_square_deviation_values)))*sum(amino_acid_square_deviation_values))
                    amino_acid_square_deviation_values.clear()
                    if rmsd>float(set_threshold):
                        text_area.insert(tk.INSERT,f'{experimental_split[0]} had a rmsd of {rmsd}\n')
        #Both files are saved for use in other programs
        os.chdir(save_directory)
        with open(save_file_sparta,'w') as file, open(save_file_peaklist,'w') as file2:
            for stuff_to_write in sparta_file_boundaries:
                file.write(stuff_to_write+'\n')
            for stuff_to_write2 in peaklist_filtered_to_match_sparta:
                    file2.write(stuff_to_write2+'\n')
#sparta file formatter
import re
import os


#This creates a sequence list that will later be used to filter residues in the sparta file outside the range we want
def create_seq_list(seq_file,seq_directory,seq_start):
    os.chdir(seq_directory)
    amino_acid_count=(0+seq_start)-1
    sequence_list=[]
    with open(seq_file) as sequence_file:
        for amino_acid in sequence_file:
            stripped_amino_acid=amino_acid.strip().upper()
            for word in stripped_amino_acid:
                amino_acid_count+=1
                sequence_list.append(str(amino_acid_count)+word)
    return sequence_list
#SPARTA files contain a lot of miscellanious info, this removes that and only extracts the residue type, number, atom type, chemical shift, and error values
#Additioanlly, prolines only contain info for 4 atom types, placeholders are set in for the nitrogen and hydrogen
def format_sparta(sparta_file,sparta_directory):
    os.chdir(sparta_directory)
    sparta_file_list1=[]
    proline_counter=0
    with open(sparta_file) as sparta_predictions:
        for line in sparta_predictions:
            modifier=line.strip().upper()
            if re.findall('^\d+',modifier):
                A=modifier.split()
                del A[5:8]
                del A[3]
                A[0:3]=["".join(A[0:3])]
                joined=" ".join(A)
                proline_searcher=re.search('\BP',joined)
                if proline_searcher != None:
                    proline_counter+=1
                    proline_count=re.search('^\d+',joined)
                    if proline_counter<2:
                        sparta_file_list1.append(f'{proline_count.group(0)}PN'+' 1000'+' 1000')
                    else:
                        if proline_counter == 4:
                            sparta_file_list1.append(joined)
                            sparta_file_list1.append(f'{proline_count.group(0)}PHN'+' 1000'+' 1000')
                            proline_counter=0
                            continue
                sparta_file_list1.append(joined)
    return sparta_file_list1

#The user may have a protein that has a mutation, causing the sequence of the sparta file to differ from theirs
#The sparta predicted value for that mutant is useless, thus it is replaced with a placeholder
def add_mutation(mutation_list1,mutation_list2,sparta_file,sparta_directory):
    sparta_file_list2=[]
    if mutation_list1==() or mutation_list2==():
        for amino_acids in format_sparta(sparta_file,sparta_directory):
            sparta_file_list2.append(amino_acids)
    else:
        for mutations,mutations2 in zip(mutation_list1,mutation_list2):
            for amino_acids in format_sparta(sparta_file,sparta_directory):
                if re.findall(mutations,amino_acids):
                    splitting=amino_acids.split()
                    mutation=re.sub(mutations,mutations2,splitting[0])
                    mutation_value=re.sub('\d+.\d+',' 1000',splitting[1])
                    mutation_value2=re.sub('\d+.\d+',' 1000',splitting[2])
                    mutation_replacement=mutation+mutation_value+mutation_value2
                    sparta_file_list2.append(mutation_replacement)
                else:
                    sparta_file_list2.append(amino_acids)
    return sparta_file_list2
#The SPARTA file may have residues beyond the scope of the users protein, those residues are filtered out
def filter_sparta_using_seq(seq_file,seq_directory,mutation_list1,mutation_list2,sparta_file,sparta_directory,seq_start):
    sparta_file_list3=[]
    sparta_comparison=create_seq_list(seq_file,seq_directory,seq_start)
    for aa in add_mutation(mutation_list1,mutation_list2,sparta_file,sparta_directory):
        modifiers=aa.strip()
        splitter=modifiers.split()
        searcher=re.search('^\d+[A-Z]',splitter[0])
        compiler=re.compile(searcher.group(0))
        sparta_sequence_comparison=list(filter(compiler.match,sparta_comparison))
        if sparta_sequence_comparison != []:
            sparta_file_list3.append(aa)

    return sparta_file_list3

#The first amino acid and last amino acid will only have 4 and 5 atom respectively, breaking the rule of 6
#If the user picks somewhere in the middle of the protein, than this is not the case, thus a check is done, and if the entire protein is not divisible by 6
#The sides are removed
def check_sparta_file_boundaries(seq_file,seq_directory,mutation_list1,mutation_list2,sparta_file,sparta_directory,seq_start):
    residue_number=[]
    number_of_residues_looped_through=0
    sparta_filtered_list=filter_sparta_using_seq(seq_file,seq_directory,mutation_list1,mutation_list2,sparta_file,sparta_directory,seq_start)
    for checker in sparta_filtered_list:
        remove_whitespace=checker.strip()
        split_values=remove_whitespace.split()
        exctract_residue_number=re.search('^\d+',split_values[0])
        residue_number.append(exctract_residue_number.group(0))
        number_of_residues_looped_through+=1
        if number_of_residues_looped_through==5:
            if int(exctract_residue_number.group(0))==int(residue_number[0]):
                break
            else:
                del sparta_filtered_list[0:4]
                break
    if len(sparta_filtered_list)%6 != 0:
        del sparta_filtered_list[-5:-1]

    return sparta_filtered_list
#nmrstar
import re
import os


#The NMRSTAR file is sorted HA,C,CA,CB,H,N, we want to format it N,HA,C<CA,CB,H
#The below function stores the residue number of each amino acid, then stores the appropriate atom in the appropriate list
#Using the residue_number_list we will know when we have moved on to the next amino acids
#When you move onto the next amino acid, the previous amino acids atoms are sorted into the appropriate order
def atom_ordering(exctracted_and_compiled_data):
    sorted_atom_types=[]
    residue_number_list=[]
    hydrogen_value=[]
    nitrogen_value=[]
    side_chain_cabonyl_values=[]
    x=0
    for amino_acids in exctracted_and_compiled_data:
        splitter2=amino_acids.split()
        x+=1
        if x >= 2:
            if splitter2[0] != residue_number_list[0]:
                list_compiler=nitrogen_value+side_chain_cabonyl_values+hydrogen_value
                sorted_atom_types.append(list_compiler)
                residue_number_list.clear()
                hydrogen_value.clear()
                nitrogen_value.clear()
                side_chain_cabonyl_values.clear()
                residue_number_list.append(splitter2[0])
                if splitter2[2] == 'H':
                    hydrogen_value.append(amino_acids)
                elif splitter2[2] == 'N':
                    nitrogen_value.append(amino_acids)
                else:
                    side_chain_cabonyl_values.append(amino_acids)
            else:
                if splitter2[2] == 'H':
                    hydrogen_value.append(amino_acids)
                elif splitter2[2] == 'N':
                    nitrogen_value.append(amino_acids)
                else:
                    side_chain_cabonyl_values.append(amino_acids)
        else:
            residue_number_list.append(splitter2[0])
            if splitter2[2] == 'H':
                hydrogen_value.append(amino_acids)
            elif splitter2[2] == 'N':
                nitrogen_value.append(amino_acids)
            else:
                side_chain_cabonyl_values.append(amino_acids)
    return sorted_atom_types

#Due to the above concatenation of lists, we form a list of lists that needs to be flattened_list
#Additionally, we wish to add a hyphen between the residue number and atom type that will be used for regex later
def flatten_list(exctracted_and_compiled_data):
    flattened_list=[]
    for lists in atom_ordering(exctracted_and_compiled_data):
        for elements in lists:
            splitting=elements.split()
            joined=''.join(splitting[0:2])
            flattened_list.append(joined+'-'+splitting[2]+ ' ' + splitting[3])
    return flattened_list

#Not every residue will have a chemical shift value for every atom types
#We want to fill in placeholders for all the missing data, but maintain that N,HA,C,CA,CB,H format
#At this point, every atom will only have the 6 desired atom types, in the appropriate atom order
#Therefore, we go through every atom for each amino acid, and check to see if we have data for that atom types in the N,HA,C order
def fill_empty_data(exctracted_and_compiled_data):
    missing_values_added=[]
    atom_value_holder=[]
    count=0
    for values in flatten_list(exctracted_and_compiled_data):
        atom_find=re.search('^-*\d+[A-Z]',values)
        count+=1
        atom_value_holder.append(atom_find.group(0))
        if count == 1:
            if re.findall('-N',values) != []:
                missing_values_added.append(values+'\n')
            else:
                missing_values_added.append(atom_value_holder[0]+'-N'+' 1000'+'\n')
                count+=1
        if count == 2:
            if re.findall('-HA',values) != []:
                missing_values_added.append(values+'\n')
            else:
                missing_values_added.append(atom_value_holder[0]+'-HA'+' 1000'+'\n')
                count+=1
        if count == 3:
            if re.findall('-C\s',values) != []:
                missing_values_added.append(values+'\n')
            else:
                missing_values_added.append(atom_value_holder[0]+'-C'+' 1000'+'\n')
                count+=1
        if count == 4:
            if re.findall('-CA',values) != []:
                missing_values_added.append(values+'\n')
            else:
                missing_values_added.append(atom_value_holder[0]+'-CA'+' 1000'+'\n')
                count+=1
        if count == 5:
            if re.findall('-CB',values) != []:
                missing_values_added.append(values+'\n')
            else:
                missing_values_added.append(atom_value_holder[0]+'-CB'+' 1000'+'\n')
                count+=1
        if count == 6:
            if re.findall('-H\s',values) != []:
                missing_values_added.append(values+'\n')
                count=0
                atom_value_holder.clear()
            else:
                missing_values_added.append(atom_value_holder[0]+'-H'+' 1000'+'\n')
                atom_value_holder.clear()
                if re.findall('-N',values) != []:
                    missing_values_added.append(values+'\n')
                    count=1
                if re.findall('-HA',values) != []:
                    missing_values_added.append(atom_find.group(0)+'-N'+' 1000'+'\n')
                    missing_values_added.append(values+'\n')
                    count=2
                if re.findall('-C',values) != []:
                    missing_values_added.append(atom_find.group(0)+'-N'+' 1000'+'\n')
                    missing_values_added.append(atom_find.group(0)+'-HA'+' 1000'+'\n')
                    missing_values_added.append(values+'\n')
                    count=3
                if re.findall('-CA',values) != []:
                    missing_values_added.append(atom_find.group(0)+'-N'+' 1000'+'\n')
                    missing_values_added.append(atom_find.group(0)+'-HA'+' 1000'+'\n')
                    missing_values_added.append(atom_find.group(0)+'-C'+' 1000'+'\n')
                    missing_values_added.append(values+'\n')
                    count=4
                if re.findall('-CB',values) != []:
                    missing_values_added.append(atom_find.group(0)+'-N'+' 1000'+'\n')
                    missing_values_added.append(atom_find.group(0)+'-HA'+' 1000'+'\n')
                    missing_values_added.append(atom_find.group(0)+'-C'+' 1000'+'\n')
                    missing_values_added.append(atom_find.group(0)+'-CA'+' 1000'+'\n')
                    missing_values_added.append(values+'\n')
                    count=5
    return missing_values_added

#Glycines do not have CBs, and they have additional HA. The above script will add an CB, this creates a new list without it
def add_glycine_HA(exctracted_and_compiled_data):
    glycine_search_list=[]
    for stuff in fill_empty_data(exctracted_and_compiled_data):
        if re.findall('\BG-HA',stuff) != []:
            splitting=stuff.split()
            glycine_search_list.append(stuff)
            glycine_search_list.append(splitting[0]+'2'+' 1000'+'\n')
        elif re.findall('\BG-CB',stuff) != []:
            pass
        else:
            glycine_search_list.append(stuff)
    return glycine_search_list


#This function creates a dictionary of residue numbers to residue type, that will be used below
dict={}
def dict_create(seq_file,seq_start,seq_directory):
    os.chdir(seq_directory)
    x=(0+seq_start)-1
    global dict
    dict={}
    with open(seq_file) as sequence_file:
        for line in sequence_file:
            white_spaces_removed=line.strip().upper()
            for word in white_spaces_removed:
                x+=1
                dict[x]=word

#The above function filled in missing data only for amino acids that had some data, but were missing data for other atom types
#This fills in placeholders for amino acids that have no data for any atom type
def fill_missing_data(exctracted_and_compiled_data,seq_start):
    outskirts_added=[]
    current_amino_acid=[]
    x=0
    y=0
    for atoms in add_glycine_HA(exctracted_and_compiled_data):
        A=re.search('^-*\d+',atoms)
        outskirts_added.append(atoms)
        x+=1
        y+=1
        if x == 6:
            if len(current_amino_acid)>0:
                if int(current_aa_residue_number) == (int(current_amino_acid[0])+1):
                    x=0
                    current_amino_acid.clear()
                    current_amino_acid.append(current_aa_residue_number)
                    pass
                else:
                    number_of_missing_amino_acid=int(current_amino_acid[0])+1
                    offset=0
                    while number_of_missing_amino_acid != int(current_aa_residue_number):
                        outskirts_added.insert((y+offset-6),f'{number_of_missing_amino_acid}{dict[number_of_missing_amino_acid]}N-H' + ' 1000' +'\n')
                        outskirts_added.insert((y+offset-6),f'{number_of_missing_amino_acid}{dict[number_of_missing_amino_acid]}N-CB' + ' 1000' +'\n')
                        outskirts_added.insert((y+offset-6),f'{number_of_missing_amino_acid}{dict[number_of_missing_amino_acid]}N-CA' + ' 1000' +'\n')
                        outskirts_added.insert((y+offset-6),f'{number_of_missing_amino_acid}{dict[number_of_missing_amino_acid]}N-C' + ' 1000' +'\n')
                        outskirts_added.insert((y+offset-6),f'{number_of_missing_amino_acid}{dict[number_of_missing_amino_acid]}N-HA' + ' 1000' +'\n')
                        outskirts_added.insert((y+offset-6),f'{number_of_missing_amino_acid}{dict[number_of_missing_amino_acid]}N-HN' + ' 1000' + '\n')
                        number_of_missing_amino_acid+=1
                        offset+=6
                    x=0
                    y+=offset
                    current_amino_acid.clear()
                    current_amino_acid.append(current_aa_residue_number)
            else:
                current_amino_acid.append(current_aa_residue_number)
                x=0
    return outskirts_added
#NMRSTAR file input (this is only a portion to get an idea on the format
Content for NMR-STAR saveframe, "assigned_chem_shift_list_1"
    save_assigned_chem_shift_list_1
   _Assigned_chem_shift_list.Sf_category                   assigned_chemical_shifts
   _Assigned_chem_shift_list.Sf_framecode                  assigned_chem_shift_list_1
   _Assigned_chem_shift_list.Entry_ID                      26909
   _Assigned_chem_shift_list.ID                            1
   _Assigned_chem_shift_list.Sample_condition_list_ID      1
   _Assigned_chem_shift_list.Sample_condition_list_label   $sample_conditions_1
   _Assigned_chem_shift_list.Chem_shift_reference_ID       1
   _Assigned_chem_shift_list.Chem_shift_reference_label    $chemical_shift_reference_1
   _Assigned_chem_shift_list.Chem_shift_1H_err             .
   _Assigned_chem_shift_list.Chem_shift_13C_err            .
   _Assigned_chem_shift_list.Chem_shift_15N_err   
...
#part we are interested in
      1      .   1   1   2     2     SER   HA     H   1    4.477     0.003   .   1   .   .   .   .   .   -1    Ser   HA     .   26909   1
      2      .   1   1   2     2     SER   HB2    H   1    3.765     0.001   .   1   .   .   .   .   .   -1    Ser   HB2    .   26909   1
      3      .   1   1   2     2     SER   HB3    H   1    3.765     0.001   .   1   .   .   .   .   .   -1    Ser   HB3    .   26909   1
      4      .   1   1   2     2     SER   C      C   13   173.726   0.2     .   1   .   .   .   .   .   -1    Ser   C      .   26909   1
      5      .   1   1   2     2     SER   CA     C   13   58.16     0.047   .   1   .   .   .   .   .   -1    Ser   CA     .   26909   1
      6      .   1   1   2     2     SER   CB     C   13   64.056    0.046   .   1   .   .   .   .   .   -1    Ser   CB     .   26909   1
      7      .   1   1   3     3     HIS   H      H   1    8.357     0.004   .   1   .   .   .   .   .   0     His   H      .   26909   1
      8      .   1   1   3     3     HIS   HA     H   1    4.725     0.003   .   1   .   .   .   .   .   0     His   HA     .   26909   1
      9      .   1   1   3     3     HIS   HB2    H   1    3.203     0.003   .   2   .   .   .   .   .   0     His   HB2    .   26909   1
      10     .   1   1   3     3     HIS   HB3    H   1    2.996     0.005   .   2   .   .   .   .   .   0     His   HB3    .   26909   1
      11     .   1   1   3     3     HIS   C      C   13   174.33    0.2     .   1   .   .   .   .   .   0     His   C      .   26909   1
      12     .   1   1   3     3     HIS   CA     C   13   55.353    0.044   .   1   .   .   .   .   .   0     His   CA     .   26909   1
      13     .   1   1   3     3     HIS   CB     C   13   31.166    0.043   .   1   .   .   .   .   .   0     His   CB     .   26909   1
      14     .   1   1   3     3     HIS   N      N   15   120.402   0.041   .   1   .   .   .   .   .   0     His   N      .   26909   1
#SPARTA file format (again, only an excerpt)
REMARK SPARTA+ Protein Chemical Shift Prediction Table
REMARK  All chemical shifts are reported in ppm:
...
#part we are interested in
   1    M   HA     0.000     4.384     4.480    -0.161     0.000     0.227
   1    M    C     0.000   176.242   176.300    -0.096     0.000     1.140
   1    M   CA     0.000    55.217    55.300    -0.139     0.000     0.988
   1    M   CB     0.000    32.488    32.600    -0.187     0.000     1.302
   2    I    N     1.287   121.802   120.570    -0.092     0.000     2.680
   2    I   HA    -0.123     4.012     4.170    -0.058     0.000     0.286
   2    I    C    -0.818   175.259   176.117    -0.066     0.000     1.144
...
\$\endgroup\$

1 Answer 1

2
\$\begingroup\$

Working directory

It's not necessary to do this:

    os.chdir(nmrstarfile_directory)

and having other code rely on the working directory makes that code more fragile and debugging trickier. pathlib has excellent facilities for building full paths off of a base path.

Regular expressions

This regex:

           extract_data_only=re.search(r'\b\d+\s+[A-Z]{3}\s+\w+\s+\w+\s+\d+\s+\d+',modifier)

would benefit from being re.compile'd outside of your loops - maybe as a global constant, or at the least near the top of the function. That way you don't have to re-compile it on every loop iteration.

Unpacking

                amino_acid_number=str(int(split_data[5])+int(seq_start)-1)
                residue_type=split_data[6]
                atom_type=split_data[7]
                converted=acid_map[residue_type]
                chemical_shift=split_data[10]

if you only need items 5-10, then

amino_acid, residue_type, atom_type, _, _, chemical_shift = split_data[5:11]

Generally, you should avoid repeated references to difficult-to-understand index expressions like splitter2[0]. Attempt to give them their own meaningfully-named variable.

Set membership

if atom_type == 'N' or atom_type == 'HA' or atom_type =='CA' or atom_type == 'CB' or atom_type=='H' or atom_type=='C':
                

can be

if atom_type in {'N', 'HA', 'CA', 'CB', 'H', 'C'}:

That set should likely be stored outside of the function as a constant.

Imports

Don't do these:

    from sparta_file_formatter import check_sparta_file_boundaries
    from nmrstar import dict_create
    from nmrstar import fill_missing_data

in the middle of your function. Do them at the top of the file.

String interpolation

atom_value_holder[0]+'-C'+' 1000'+'\n'

can be

f'{atom_value_holder[0]}-C 1000\n'

Even if you didn't use an f-string, there is no need to separate those last three string literals into concatenations.

Extend

                missing_values_added.append(atom_find.group(0)+'-N'+' 1000'+'\n')
                missing_values_added.append(atom_find.group(0)+'-HA'+' 1000'+'\n')
                missing_values_added.append(atom_find.group(0)+'-C'+' 1000'+'\n')
                missing_values_added.append(atom_find.group(0)+'-CA'+' 1000'+'\n')

should be

atom = atom_find.group(0)
missing_values_added.extend((
    f'{atom}-N 1000\n',
    f'{atom}-HA 1000\n',
    f'{atom}-C 1000\n',
    f'{atom}-CA 1000\n',
))

Checking for any match

Do not use findall here:

re.findall('\BG-CB',stuff) != []

Use search. If it returns None, there are no hits; otherwise there is at least one hit; pair this with is not None.

Shadowing

This:

dict={}

is nasty, and setting you up for failure. dict is a (very commonly used) built-in name, so don't shadow it with your own variable - particularly at the global level.

\$\endgroup\$
7
  • \$\begingroup\$ A couple of questions. 1) Why would I use pathlib versus os? I use it to determine the directory, since the user may have files in different directories. 2) I'm a bit confused by what you mean by "recompile it every loop". 3) For unpacking, does putting an underscore _, mean that value will not be taken? Because I don't want everything from 5:11, only 6,7 and 10. 4) For membership, why would I leave these outside of the function as constants? What difference would it make inside vs outside? 5) I'm not too faimilar with extend, however it appears you can't index it (which is crucial). \$\endgroup\$
    – samman
    Commented Jun 29, 2020 at 22:25
  • \$\begingroup\$ Why would I use pathlib versus os - pathlib.Path is a nice, self-contained, object-oriented representation of a path, whereas the os methods use an older, more procedural style. It's a matter of convenience and aesthetics. \$\endgroup\$
    – Reinderien
    Commented Jun 29, 2020 at 22:26
  • \$\begingroup\$ recompile it every loop - if you call re.search / re.match, it has to parse your regex all over again, and in the middle of a loop (such as iterating over file lines) this may prove costly. \$\endgroup\$
    – Reinderien
    Commented Jun 29, 2020 at 22:28
  • \$\begingroup\$ does putting an underscore _, mean that value will not be taken? By convention, yes; it indicates you won't use that variable but it has to exist for the syntax to be correct. \$\endgroup\$
    – Reinderien
    Commented Jun 29, 2020 at 22:29
  • \$\begingroup\$ For membership, why would I leave these outside of the function as constants? If that set of strings is fundamental to the program and reused anywhere else; I'd be surprised if that's the only place that could use it. \$\endgroup\$
    – Reinderien
    Commented Jun 29, 2020 at 22:30

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.