Protein Docking: get accessible surface area

get_asa.jl

function get_sphere_points(n::Int)
    dl = pi * (3.0 - sqrt(5.0))
    dz = 2.0 / n

    longitude = 0
    z = 1 - dz / 2

    coords = zeros(n, 3)
    for k in 1:n
        r = sqrt((1 - z * z))
        coords[k, 1] = cos(longitude) * r
        coords[k, 2] = sin(longitude) * r
        coords[k, 3] = z
        z -= dz
        longitude += dl
    end

    coords
end

# based on https://biopython.org/docs/dev/api/Bio.PDB.SASA.html#Bio.PDB.SASA.ShrakeRupley
function get_asa(atoms::Vector, radius_list::Vector, dist_matrix::Symmetric{Float64, Matrix{Float64}}, sphere_points = nothing)
    if sphere_points == nothing
        sphere_points = get_sphere_points(100)
    end

    probe_radius = 1.4

    max_radius = 2.0 * ( maximum(radius_list) + probe_radius )
    neighbors = [ findall(<=(max_radius + 0.01), dist_matrix[i,:]) for i in 1:(size(dist_matrix)[1])]

    atom_asa = zeros(length(atoms))

    for atom_idx in 1:length(atoms)
        atom_radius = radius_list[atom_idx]
        neighbor_indexes = neighbors[atom_idx]

        atom_sphere_points = [atoms[atom_idx] + ( (atom_radius + probe_radius) * sphere_points[i,:]) for i in 1:(size(sphere_points)[1])]
        
        points_inaccessible = [false for n in 1:n_points]

        for neighbor_idx in neighbor_indexes
            (atom_idx == neighbor_idx) && continue

            neighbor_atom = atoms[neighbor_idx]
            neighbor_radius = radius_list[neighbor_idx]

            neighbor_contains_points = [ (norm(neighbor_atom - pt) < (neighbor_radius + probe_radius)) for pt in atom_sphere_points]
            points_inaccessible = [ (neighbor_contains_points[pt_idx] || points_inaccessible[pt_idx]) for pt_idx in 1:n_points]
        end
        
        percent_accessible = 1.0 - ( sum(points_inaccessible) / n_points )

        surface_area = 4 * pi * ( (atom_radius + probe_radius) ^2 )
        atom_asa[atom_idx] = percent_accessible * surface_area
    end
    atom_asa
end

# Van der Waal radius of each atom type
atom_radius_dict = Dict("C"=> 1.700, "N"=> 1.550, "O"=> 1.520, "S"=> 1.800, "NA"=> 2.270, "MG"=> 1.730)

ligand_radius_list = [atom_radius_dict[strip(atom.name[1:2])] for atom in ligand_atoms]
receptor_radius_list = [atom_radius_dict[strip(atom.name[1:2])] for atom in receptor_atoms]

# calculate accessible surface area
n_points = 1000
sphere_points = get_sphere_points(n_points)

# calculate surface area
ligand_atoms_asa = get_asa(ligand_coords, ligand_radius_list, ligand_dist_matrix, sphere_points)
receptor_atoms_asa = get_asa(receptor_coords, receptor_radius_list, receptor_dist_matrix, sphere_points)