Protein Docking: load protein data

load_protein_data.jl

using BioStructures

proteins_dict = Dict()
protein_complex_pdb_id = "3WD5"
figures_path = mkpath(string(protein_complex_pdb_id, "_figures"))

for designation in ("l_u", "r_u", "l_b", "r_b")
    protein = read(string("./protein_data/", protein_complex_pdb_id, '_', designation, ".pdb"), PDBFormat)
    length(protein.models) == 1 || throw(DomainError(protein.models), "pdb files should only contain one model")
    proteins_dict[designation] = protein[1]
end

# spatial midpoints of unbound proteins
ligand_midpt, receptor_midpt = calc_midpoint(ligand_coords), calc_midpoint(receptor_coords)
ligand_coords, receptor_coords = [c - ligand_midpt for c in ligand_coords], [c - receptor_midpt for c in receptor_coords]

# gather all amino acids into a dictionary where the keys are (chain ID, amino acid ID) so the same amino acids can be compared
# between the bound and unbound versions of the proteins
all_residues_dict = Dict()

for k in keys(proteins_dict)
    protein = proteins_dict[k]
    residues_dict = Dict()
    for chain_id in keys(protein.chains)
        for res in collectresidues(protein[chain_id])
            # ignore non-protein molecules
            res.het_res && continue
            coords = [atom.coords - receptor_midpt for atom in values(res.atoms)]
            residues_dict[(chain_id, res.number)] = Dict("coords" => coords, "res_num" => res.number, "is_interface" => false, "c_a_coords" => res.atoms[" CA "].coords - receptor_midpt )
        end
    end
    all_residues_dict[k] = residues_dict
end