bfactor_plot_residue.py

#bfactor_plot_residue

import os
import matplotlib.pyplot as plt

ensemble_paths = ["model_5_multimer_v2__115",  "model_5_multimer_v2__156", "model_5_multimer_v2_no_templates_194", "model_5_multimer_v2_no_templates_r9_76", "model_5_multimer_v2__129", "model_5_multimer_v2__188", "model_5_multimer_v2_no_templates_86", "model_5_multimer_v2_no_templates_r9_85", "model_5_multimer_v2__130", "model_5_multimer_v2_no_templates_170", "model_5_multimer_v2_no_templates_r9_133", "model_5_multimer_v2_no_templates_r9_92", "model_5_multimer_v2__150", "model_5_multimer_v2_no_templates_184", "model_5_multimer_v2_no_templates_r9_65"]


def get_average_bfactor(pdb_path, residue_name):
    found = False
    for p in ensemble_paths:
        if pdb_path.find(p) != -1:
            found = True
            break
    if not(found):
        x = 0
        #return None
    """
    Extracts the average B-factor for a specific residue in a PDB file.
    
    Parameters:
    - pdb_path: Path to the PDB file.
    - residue_name: Name of the residue (e.g., "SER A  65").
    
    Returns:
    - The average B-factor as a float, or None if the residue is not found.
    """
    b_factors = []
    with open(pdb_path, 'r') as pdb_file:
        for line in pdb_file:
            if line.startswith("ATOM") and residue_name in line:
                b_factor = float(line[60:66].strip())
                b_factors.append(b_factor)
    
    if b_factors:
        return sum(b_factors) / len(b_factors)
    else:
        return None




def plot_bfactor_distribution(directory, residue_name):
    """
    Plots the distribution of average B-factor values for a specific residue across PDB files in a directory.
    
    Parameters:
    - directory: Directory containing PDB files.
    - residue_name: Name of the residue (e.g., "SER A  65").
    """
    average_bfactors = []
    for filename in os.listdir(directory):
        if filename.endswith(".pdb"):
            pdb_path = os.path.join(directory, filename)
            avg_bfactor = get_average_bfactor(pdb_path, residue_name)
            if avg_bfactor is not None:
                average_bfactors.append(avg_bfactor)
    
    # Plotting
    residue_name = 'SER A  65'  # Specify the residue of interest
    if len(average_bfactors) > 1000:
        plt.hist(average_bfactors, bins=100, alpha=0.75)
        plt.title(f'Distribution of plddt for {residue_name} in 6000 AFS structures')
    else:
        plt.hist(average_bfactors, bins=10, alpha=0.75)
        plt.title(f'Distribution of plddt for {residue_name} in final AFS ensemble')
    plt.xlim((50, 100))
    plt.xlabel('Average B-factor')
    plt.ylabel('Frequency')
    plt.show()


# Example usage
directory = './7jq8/'  # Update with the actual directory path
residue_name = 'SER A  40'  # Specify the residue of interest
plot_bfactor_distribution(directory, residue_name)