<p>using your code :</p>
<pre><code>import numpy as np
import pandas as pd
from Bio import PDB
from Bio.PDB.Polypeptide import is_aa

def calc_distance_matrix(atoms):
    &quot;&quot;&quot;Calculate the distance matrix for a list of atoms.&quot;&quot;&quot;
    num_atoms = len(atoms)
    distance_matrix = np.zeros((num_atoms, num_atoms))

    for i, atom1 in enumerate(atoms):
        for j, atom2 in enumerate(atoms):
            distance_matrix[i, j] = atom1 - atom2  # This uses the overridden '-' operator to calculate distance

    return distance_matrix

def main(pdb_filename, output_csv_filename):
    # Create a parser and parse the structure
    parser = PDB.PDBParser(QUIET=True)
    structure = parser.get_structure('protein', pdb_filename)

    # Get the first model and first chain
    model = structure[0]
    chain = list(model.get_chains())[0]

    # Extract C-alpha atoms
    ca_atoms = [atom for atom in chain.get_atoms() if atom.get_id() == 'CA' and is_aa(atom.get_parent())]

    # Calculate the distance matrix
    distance_matrix = calc_distance_matrix(ca_atoms)
    
    #print dis\tance matrix
    
    print(distance_matrix)

    # Save the distance matrix to a CSV file
    df = pd.DataFrame(distance_matrix)
    df.to_csv(output_csv_filename, index=False)
    
    
    
main('test.pdb'  , 'test_out.csv')
</code></pre>
<p>with input `test.pdb' :</p>
<pre><code>ATOM      1  CA  ALA A   1      00.000  00.000  00.000  1.00  0. 0           C
ATOM      2  CA  ALA A   2      01.000  00.000  00.000  1.00  0. 0           C
ATOM      3  CA  ALA A   3      01.000  01.000  00.000  1.00  0. 0           C
ATOM      3  CA  ALA A   4      01.000  01.000  01.000  1.00  0. 0           C

</code></pre>
<p>got :</p>
<pre><code>[[0.         1.         1.41421354 1.73205078]
 [1.         0.         1.         1.41421354]
 [1.41421354 1.         0.         1.        ]
 [1.73205078 1.41421354 1.         0.        ]]
</code></pre>
<p>using 'test_2.pdb' :</p>
<pre><code>ATOM      1  CA  ALA A   1      00.000  00.000  00.000  1.00  0. 0           C
ATOM      2  CA  ALA A   2      00.000  00.000  00.000  1.00  0. 0           C
ATOM      3  CA  ALA A   3      00.000  00.000  00.000  1.00  0. 0           C
ATOM      3  CA  ALA A   4      00.000  00.000  00.000  1.00  0. 0           C
</code></pre>
<p>got :</p>
<pre><code>[[0. 0. 0. 0.]
 [0. 0. 0. 0.]
 [0. 0. 0. 0.]
 [0. 0. 0. 0.]]
</code></pre>
<p>seems OK to me</p>
<p>as long as <code>distance_matrix</code> is , like WIKIPEDIA says,</p>
<blockquote>
<p>In mathematics, computer science and especially graph theory, a
distance matrix is a square matrix (two-dimensional array) containing
the distances, taken pairwise, between the elements of a set.</p>
</blockquote>
<p><strong>EDITED</strong></p>
<p>To prove that your results are correct, you could check that the output produce an Euclidean Distance Matrix (<a href="https://en.wikipedia.org/wiki/Euclidean_distance_matrix" rel="nofollow noreferrer">Euclidean distance matrix Wikipedia</a> , checking that the result fulfils its properties:</p>
<p>By the fact that Euclidean distance is a metric, the matrix A has the following properties.</p>
<p>All elements on the diagonal of A are zero (i.e. it is a hollow atrix); hence the trace of A is zero</p>
<p>A is symmetric</p>
<p>........
.....</p>
<p>see this updated code:</p>
<pre><code>import numpy as np
import pandas as pd
from Bio import PDB
from Bio.PDB.Polypeptide import is_aa

def calc_distance_matrix(atoms):
    &quot;&quot;&quot;Calculate the distance matrix for a list of atoms.&quot;&quot;&quot;
    num_atoms = len(atoms)
    distance_matrix = np.zeros((num_atoms, num_atoms))

    for i, atom1 in enumerate(atoms):
        for j, atom2 in enumerate(atoms):
            distance_matrix[i, j] = atom1 - atom2  # This uses the overridden '-' operator to calculate distance

    return distance_matrix

def check_result(matrix):
    
    ## all elements are &gt;= 0
    if not (matrix &gt;= 0).all() :
    
        return False
    
    
    # matrix is symmetric
    elif not (matrix == matrix.T).all():
        
        return False
    
    # matrix diagonal is [0 , 0, ,0 ...]
    elif not (matrix.diagonal() == 0).all():
        
        return False
    
    
    else:
        
        return True

def main(pdb_filename, output_csv_filename):
    # Create a parser and parse the structure
    parser = PDB.PDBParser(QUIET=True)
    structure = parser.get_structure('protein', pdb_filename)

    # Get the first model and first chain
    model = structure[0]
    chain = list(model.get_chains())[0]

    # Extract C-alpha atoms
    ca_atoms = [atom for atom in chain.get_atoms() if atom.get_id() == 'CA' and is_aa(atom.get_parent())]

    # Calculate the distance matrix
    distance_matrix = calc_distance_matrix(ca_atoms)
    
    #print dis\tance matrix
    
    print(distance_matrix)

    # Save the distance matrix to a CSV file
    df = pd.DataFrame(distance_matrix)
    df.to_csv(output_csv_filename, index=False)
    
    print('is correct ? :' , check_result(distance_matrix))
    
main('test.pdb'  , 'test_out.csv')
</code></pre>

CC BY-SA 4.0
2024-07-05T11:49:50.300Z
what logic have u used?
user366312
CC BY-SA 4.0
2024-07-05T12:38:24.990Z
How do you verify that an alghorythm is correct ??
pippo1980
CC BY-SA 4.0
2024-07-05T13:20:44.377Z
Depends on the algorithm.
user366312
CC BY-SA 4.0
2024-07-05T15:17:47.467Z
On your aghorithm ?
pippo1980
CC BY-SA 4.0
2024-07-05T16:46:22.857Z
Updated my answer
pippo1980
CC BY-SA 4.0
2024-07-05T11:15:39.907Z
2024-07-05T16:29:14.877Z
2024-07-05T16:29:14.877Z
pippo1980
<p>I want to compute the distance matrix from PDB files.</p>
<p>So, I wrote the following Python code using <code>BioPython</code>.</p>
<pre><code>import numpy as np
import pandas as pd
from Bio import PDB
from Bio.PDB.Polypeptide import is_aa

def calc_distance_matrix(atoms):
    &quot;&quot;&quot;Calculate the distance matrix for a list of atoms.&quot;&quot;&quot;
    num_atoms = len(atoms)
    distance_matrix = np.zeros((num_atoms, num_atoms))

    for i, atom1 in enumerate(atoms):
        for j, atom2 in enumerate(atoms):
            distance_matrix[i, j] = atom1 - atom2  # This uses the overridden '-' operator to calculate distance

    return distance_matrix

def main(pdb_filename, output_csv_filename):
    # Create a parser and parse the structure
    parser = PDB.PDBParser(QUIET=True)
    structure = parser.get_structure('protein', pdb_filename)

    # Get the first model and first chain
    model = structure[0]
    chain = list(model.get_chains())[0]

    # Extract C-alpha atoms
    ca_atoms = [atom for atom in chain.get_atoms() if atom.get_id() == 'CA' and is_aa(atom.get_parent())]

    # Calculate the distance matrix
    distance_matrix = calc_distance_matrix(ca_atoms)

    # Save the distance matrix to a CSV file
    df = pd.DataFrame(distance_matrix)
    df.to_csv(output_csv_filename, index=False)
</code></pre>
<p>How can I verify that this script is computing distance matrices correctly?</p>

<proteins><distance-matrix>
How can I compute distance matrices from a PDB file?
# How can I compute distance matrices from a PDB file?

(Asked by: username_0 on 2024-07-04)

I want to compute the distance matrix from PDB files.
So, I wrote the following Python code using BioPython.
import numpy as np
import pandas as pd
from Bio import PDB
from Bio.PDB.Polypeptide import is_aa

def calc_distance_matrix(atoms):
    """Calculate the distance matrix for a list of atoms."""
    num_atoms = len(atoms)
    distance_matrix = np.zeros((num_atoms, num_atoms))

    for i, atom1 in enumerate(atoms):
        for j, atom2 in enumerate(atoms):
            distance_matrix[i, j] = atom1 - atom2  # This uses the overridden '-' operator to calculate distance

    return distance_matrix

def main(pdb_filename, output_csv_filename):
    # Create a parser and parse the structure
    parser = PDB.PDBParser(QUIET=True)
    structure = parser.get_structure('protein', pdb_filename)

    # Get the first model and first chain
    model = structure[0]
    chain = list(model.get_chains())[0]

    # Extract C-alpha atoms
    ca_atoms = [atom for atom in chain.get_atoms() if atom.get_id() == 'CA' and is_aa(atom.get_parent())]

    # Calculate the distance matrix
    distance_matrix = calc_distance_matrix(ca_atoms)

    # Save the distance matrix to a CSV file
    df = pd.DataFrame(distance_matrix)
    df.to_csv(output_csv_filename, index=False)

How can I verify that this script is computing distance matrices correctly?

**Tags:** proteins, distance-matrix

---

## 1 answers

### Answer by username_1 on 2024-07-05. Score: 4 | Accepted answer

using your code :
import numpy as np
import pandas as pd
from Bio import PDB
from Bio.PDB.Polypeptide import is_aa

def calc_distance_matrix(atoms):
    """Calculate the distance matrix for a list of atoms."""
    num_atoms = len(atoms)
    distance_matrix = np.zeros((num_atoms, num_atoms))

    for i, atom1 in enumerate(atoms):
        for j, atom2 in enumerate(atoms):
            distance_matrix[i, j] = atom1 - atom2  # This uses the overridden '-' operator to calculate distance

    return distance_matrix

def main(pdb_filename, output_csv_filename):
    # Create a parser and parse the structure
    parser = PDB.PDBParser(QUIET=True)
    structure = parser.get_structure('protein', pdb_filename)

    # Get the first model and first chain
    model = structure[0]
    chain = list(model.get_chains())[0]

    # Extract C-alpha atoms
    ca_atoms = [atom for atom in chain.get_atoms() if atom.get_id() == 'CA' and is_aa(atom.get_parent())]

    # Calculate the distance matrix
    distance_matrix = calc_distance_matrix(ca_atoms)
    
    #print dis\tance matrix
    
    print(distance_matrix)

    # Save the distance matrix to a CSV file
    df = pd.DataFrame(distance_matrix)
    df.to_csv(output_csv_filename, index=False)
    
    
    
main('test.pdb'  , 'test_out.csv')

with input `test.pdb' :
ATOM      1  CA  ALA A   1      00.000  00.000  00.000  1.00  0. 0           C
ATOM      2  CA  ALA A   2      01.000  00.000  00.000  1.00  0. 0           C
ATOM      3  CA  ALA A   3      01.000  01.000  00.000  1.00  0. 0           C
ATOM      3  CA  ALA A   4      01.000  01.000  01.000  1.00  0. 0           C


got :
[[0.         1.         1.41421354 1.73205078]
 [1.         0.         1.         1.41421354]
 [1.41421354 1.         0.         1.        ]
 [1.73205078 1.41421354 1.         0.        ]]

using 'test_2.pdb' :
ATOM      1  CA  ALA A   1      00.000  00.000  00.000  1.00  0. 0           C
ATOM      2  CA  ALA A   2      00.000  00.000  00.000  1.00  0. 0           C
ATOM      3  CA  ALA A   3      00.000  00.000  00.000  1.00  0. 0           C
ATOM      3  CA  ALA A   4      00.000  00.000  00.000  1.00  0. 0           C

got :
[[0. 0. 0. 0.]
 [0. 0. 0. 0.]
 [0. 0. 0. 0.]
 [0. 0. 0. 0.]]

seems OK to me
as long as distance_matrix is , like WIKIPEDIA says,

In mathematics, computer science and especially graph theory, a
distance matrix is a square matrix (two-dimensional array) containing
the distances, taken pairwise, between the elements of a set.

EDITED
To prove that your results are correct, you could check that the output produce an Euclidean Distance Matrix (Euclidean distance matrix Wikipedia , checking that the result fulfils its properties:
By the fact that Euclidean distance is a metric, the matrix A has the following properties.
All elements on the diagonal of A are zero (i.e. it is a hollow atrix); hence the trace of A is zero
A is symmetric
........
.....
see this updated code:
import numpy as np
import pandas as pd
from Bio import PDB
from Bio.PDB.Polypeptide import is_aa

def calc_distance_matrix(atoms):
    """Calculate the distance matrix for a list of atoms."""
    num_atoms = len(atoms)
    distance_matrix = np.zeros((num_atoms, num_atoms))

    for i, atom1 in enumerate(atoms):
        for j, atom2 in enumerate(atoms):
            distance_matrix[i, j] = atom1 - atom2  # This uses the overridden '-' operator to calculate distance

    return distance_matrix

def check_result(matrix):
    
    ## all elements are >= 0
    if not (matrix >= 0).all() :
    
        return False
    
    
    # matrix is symmetric
    elif not (matrix == matrix.T).all():
        
        return False
    
    # matrix diagonal is [0 , 0, ,0 ...]
    elif not (matrix.diagonal() == 0).all():
        
        return False
    
    
    else:
        
        return True

def main(pdb_filename, output_csv_filename):
    # Create a parser and parse the structure
    parser = PDB.PDBParser(QUIET=True)
    structure = parser.get_structure('protein', pdb_filename)

    # Get the first model and first chain
    model = structure[0]
    chain = list(model.get_chains())[0]

    # Extract C-alpha atoms
    ca_atoms = [atom for atom in chain.get_atoms() if atom.get_id() == 'CA' and is_aa(atom.get_parent())]

    # Calculate the distance matrix
    distance_matrix = calc_distance_matrix(ca_atoms)
    
    #print dis\tance matrix
    
    print(distance_matrix)

    # Save the distance matrix to a CSV file
    df = pd.DataFrame(distance_matrix)
    df.to_csv(output_csv_filename, index=False)
    
    print('is correct ? :' , check_result(distance_matrix))
    
main('test.pdb'  , 'test_out.csv')

> **Comments:**
> **username_0**: what logic have u used?  
> **username_1**: How do you verify that an alghorythm is correct ??  
> **username_0**: Depends on the algorithm.  
> **username_1**: On your aghorithm ?  
> **username_1**: Updated my answer  

---

