DSSP_visualizer.py

import seaborn as sns
import matplotlib.pyplot as plt
import numpy as np
import matplotlib.colors as mcolors
import os


def load_dssp_matrix(file_path):
    try:
        file_content = []
        with open(file_path, 'r') as file:
            # Skip header lines if any
            line_number = 0
            for line in file:
                line_number += 1
                if line_number % 10 != 0:  # Skip every line that is not a multiple of ten
                    continue
                if line.startswith('#'):
                    continue
                row = [ord(character) for character in line.rstrip()]
                file_content.append(row)
                if line_number >= 10030:  # Stop reading lines after reaching line 1000
                    break

            dssp_matrix = np.array(file_content)

        return dssp_matrix
    except Exception as e:
        print(f"Error loading DSSP matrix: {e}")
        raise ValueError("Error loading DSSP matrix. Check the file structure.")

def plot_dssp_heatmap(dssp_matrix):
    # Create a heatmap using seaborn
    sns.set()
    plt.figure(figsize=(10, 8))

    # Define colors for P, S, and ~
    colors = ['purple', 'green', 'white']
    # Create a ListedColormap
    cmap = mcolors.ListedColormap(colors)
    dssp_matrix=np.transpose(dssp_matrix)
    sns.heatmap(dssp_matrix, cmap=cmap, cbar=False, xticklabels=True, yticklabels=True)
    plt.ylabel('Residue Number')
    plt.xlabel('Frame Number')
    plt.title('DSSP Heatmap')

    # Add legend with secondary structure types
    legend_labels = {'P': 'κ-helix', 'S': 'bend', '~': 'loop'}
    legend_handles = [plt.Line2D([0], [0], marker='o', color='w', label=f"{label} - {legend_labels[label]}", markerfacecolor=color, markersize=10) for
                      label, color in zip(legend_labels.keys(), colors)]
    plt.legend(handles=legend_handles, title='Label - Secondary Structure Type')

    # # Add custom colorbar tick labels
    # cbar = plt.colorbar()
    # cbar.set_ticks([80, 83, 126])
    # cbar.set_ticklabels(['P', 'S', '~'])

    # Customize ticks and tick labels
    frame_step = 1000
    frame_numbers = range(0, len(dssp_matrix[0]) + 1, frame_step)
    plt.xticks(np.arange(0.5, len(dssp_matrix[0]) + 0.5, frame_step), frame_numbers)
    #plt.yticks(np.arange(0.5, len(dssp_matrix) + 0.5), range(1, len(dssp_matrix) + 1))
    residue_step = 50
    residue_numbers = range(0, len(dssp_matrix), residue_step)
    plt.yticks(np.arange(0.5, len(dssp_matrix) + 0.5, residue_step), residue_numbers)

    plt.show()


def plot_secondary_structure_percentages(dssp_matrix):
    total_counts = np.sum(dssp_matrix != ord('-'), axis=1) # Count residues that are not '-'

    # Grouping similar secondary structures and counting them
    counts_alpha_helix = np.sum((dssp_matrix == ord('H')) | (dssp_matrix == ord('G')) | (dssp_matrix == ord('I')) | (dssp_matrix == ord('P')), axis=1)
    counts_beta_structures = np.sum((dssp_matrix == ord('B')) | (dssp_matrix == ord('E')), axis=1)
    counts_other_structures = np.sum((dssp_matrix == ord('S')) | (dssp_matrix == ord('T')) | (dssp_matrix == ord('=')) | (dssp_matrix == ord('~')), axis=1)

    # Calculating percentages
    percentages_alpha_helix = (counts_alpha_helix / total_counts) * 100
    percentages_beta_structures = (counts_beta_structures / total_counts) * 100
    percentages_other_structures = (counts_other_structures / total_counts) * 100

    plt.figure(figsize=(10, 6))
    plt.stackplot(range(len(dssp_matrix)), percentages_alpha_helix, percentages_beta_structures, percentages_other_structures,
                  labels=['Alpha Helix', 'Beta Structures', 'Other Structures'], colors=['red', 'orange', 'green'])
    plt.xlabel('Time')
    plt.ylabel('Percentage')
    plt.title('Secondary Structure Percentages')
    plt.legend(loc='upper left')
    plt.show()


def plot_multiple_dssp_files_across_takes(folder_path):
    file_names = [f for f in os.listdir(folder_path) if f.endswith('.dat')]
    num_files = len(file_names)

    num_cols = 1  # Number of columns for subplots
    num_rows = (num_files + num_cols - 1) // num_cols  # Calculate number of rows

    fig, axes = plt.subplots(num_rows, num_cols, figsize=(7, 4))
    axes = axes.flatten()  # Flatten the axes array for easier iteration

    for i, file_name in enumerate(file_names):
        file_path = os.path.join(folder_path, file_name)
        dssp_matrix = load_dssp_matrix(file_path)
        total_counts = np.sum(dssp_matrix != ord('-'), axis=1)  # Count residues that are not '-'
        counts_alpha_helix = np.sum(
            (dssp_matrix == ord('H')) | (dssp_matrix == ord('G')) | (dssp_matrix == ord('I')) | (
                        dssp_matrix == ord('P')), axis=1)
        counts_beta_structures = np.sum((dssp_matrix == ord('B')) | (dssp_matrix == ord('E')), axis=1)
        counts_other_structures = np.sum(
            (dssp_matrix == ord('S')) | (dssp_matrix == ord('T')) | (dssp_matrix == ord('=')) | (
                        dssp_matrix == ord('~')), axis=1)
        percentages_alpha_helix = (counts_alpha_helix / total_counts) * 100
        percentages_beta_structures = (counts_beta_structures / total_counts) * 100
        percentages_other_structures = (counts_other_structures / total_counts) * 100

        ax = axes[i]
        ax.stackplot(range(len(dssp_matrix)), percentages_alpha_helix, percentages_beta_structures,
                     percentages_other_structures,
                     labels=['Alpha Helix', 'Beta Structures', 'Other Structures'],
                     colors=['red', 'orange', 'green'])
        ax.set_ylabel(file_name[0:6])

        ax.set_yticks([0, 50, 100])  # Set y-axis ticks
        # Set y-axis tick labels with the % symbol
        y_tick_labels = [f'{tick} %' for tick in ax.get_yticks()]
        ax.set_yticklabels(y_tick_labels)

        # Remove x-axis labels except for the bottom subplot
        if i < len(axes) - 1:
            ax.set_xticklabels([])

    # Remove any unused subplots
    for j in range(i + 1, len(axes)):
        fig.delaxes(axes[j])

    # Set x-axis label for the bottom subplot
    axes[-1].set_xlabel('Time (ns)')

    # Create a single legend outside the loop
    handles, labels = ax.get_legend_handles_labels()
    fig.legend(handles, labels, loc='upper right')
    ax.set_xlabel('Time (ns)')

    plt.tight_layout()
    plt.show()


def plot_multiple_dssp_files_across_ForceFields(folder_path):
    # List of files to plot in the correct order
    list_of_files = [
        'AMBER_1_RbsB_dssp.dat', 'CHARMM_1_RbsB_dssp.dat', 'SIRAH_1_RbsB_dssp.dat',
        'AMBER_1_FkpA_dssp.dat', 'CHARMM_1_FkpA_dssp.dat', 'SIRAH_1_FkpA_dssp.dat',
        'AMBER_1_SpeA_dssp.dat', 'CHARMM_1_SpeA_dssp.dat', 'SIRAH_1_SpeA_dssp.dat'
    ]

    num_files = len(list_of_files)
    num_cols = 3  # Number of columns for subplots
    num_rows = (num_files + num_cols - 1) // num_cols  # Calculate number of rows

    fig, axes = plt.subplots(num_rows, num_cols, figsize=(6, 5))
    axes = axes.flatten()  # Flatten the axes array for easier iteration

    for i, file_name in enumerate(list_of_files):
        file_path = os.path.join(folder_path, file_name)
        dssp_matrix = load_dssp_matrix(file_path)
        total_counts = np.sum(dssp_matrix != ord('-'), axis=1)  # Count residues that are not '-'
        counts_alpha_helix = np.sum(
            (dssp_matrix == ord('H')) | (dssp_matrix == ord('G')) | (dssp_matrix == ord('I')) | (
                    dssp_matrix == ord('P')), axis=1)
        counts_beta_structures = np.sum((dssp_matrix == ord('B')) | (dssp_matrix == ord('E')), axis=1)
        counts_other_structures = np.sum(
            (dssp_matrix == ord('S')) | (dssp_matrix == ord('T')) | (dssp_matrix == ord('=')) | (
                    dssp_matrix == ord('~')), axis=1)
        percentages_alpha_helix = (counts_alpha_helix / total_counts) * 100
        percentages_beta_structures = (counts_beta_structures / total_counts) * 100
        percentages_other_structures = (counts_other_structures / total_counts) * 100

        ax = axes[i]
        ax.stackplot(range(len(dssp_matrix)), percentages_alpha_helix, percentages_beta_structures,
                     percentages_other_structures,
                     labels=['DSSP - Alpha Helix', 'DSSP - Beta Structures', 'DSSP - Other Structures'],
                     colors=['red', 'orange', 'green'])

        # Set y-axis label to the protein name only for the subplots in the first column
        if i % num_cols == 0:
            ax.yaxis.set_label_position('left')  # Set y-axis label position to the left for the first column
            ax.set_ylabel(file_name.split('_')[2].split('.')[0])

        # Set y-axis ticks and labels on the right-hand side
        ax.yaxis.tick_right()
        ax.yaxis.set_label_position('right')

        # Set y-axis ticks and labels
        ax.set_yticks([0, 50, 100])  # Set y-axis ticks
        ax.set_yticklabels([f'{tick}%' for tick in [0, 50, 100]])  # Set y-axis tick labels with the % symbol

        # # Set x-axis label for the bottom subplot
        # if i >= (num_rows - 1) * num_cols:
        #     ax.set_xlabel('Time (ns)')
        ax.set_xticks([0, 500, 1000])  # Set x-axis ticks
        # ax.set_xticklabels([0, 500, 1000])  # Set x-axis ticks
        # ax.set_xlabel('Time (ns)')


        # # Remove x-axis labels except for the bottom subplots
        # if i < num_files - num_cols:
        #     ax.set_xticklabels([])
        ax.set_xticklabels([])



    # Remove any unused subplots
    for j in range(num_files, len(axes)):
        fig.delaxes(axes[j])

    # Create a single legend outside the loop
    handles, labels = ax.get_legend_handles_labels()
    fig.legend(handles, labels, loc='upper right', bbox_to_anchor=(0.9, 0.9))

    plt.tight_layout()
    plt.show()

def calculate_deviation_per_protein(folder_path):
    list_of_files = [
        'AMBER_1_RbsB_dssp.dat', 'CHARMM_1_RbsB_dssp.dat', 'SIRAH_1_RbsB_dssp.dat',
        'AMBER_1_FkpA_dssp.dat', 'CHARMM_1_FkpA_dssp.dat', 'SIRAH_1_FkpA_dssp.dat',
        'AMBER_1_SpeA_dssp.dat', 'CHARMM_1_SpeA_dssp.dat', 'SIRAH_1_SpeA_dssp.dat'
    ]

    for file_name in list_of_files:
        file_path = os.path.join(folder_path, file_name)
        dssp_matrix = load_dssp_matrix(file_path)

        total_counts = np.sum(dssp_matrix != ord('-'), axis=1)
        counts_alpha_helix = np.sum(
            (dssp_matrix == ord('H')) | (dssp_matrix == ord('G')) | (dssp_matrix == ord('I')) | (
                        dssp_matrix == ord('P')), axis=1)
        counts_beta_structures = np.sum((dssp_matrix == ord('B')) | (dssp_matrix == ord('E')), axis=1)
        counts_other_structures = np.sum(
            (dssp_matrix == ord('S')) | (dssp_matrix == ord('T')) | (dssp_matrix == ord('=')) | (
                        dssp_matrix == ord('~')), axis=1)

        percentages_alpha_helix = (counts_alpha_helix / total_counts) * 100
        percentages_beta_structures = (counts_beta_structures / total_counts) * 100
        percentages_other_structures = (counts_other_structures / total_counts) * 100

        std_alpha_helix = np.std(percentages_alpha_helix)
        std_beta_structures = np.std(percentages_beta_structures)
        std_other_structures = np.std(percentages_other_structures)

        protein_name = file_name.split('_')[2].split('.')[0]
        force_field = file_name.split('_')[0]

        print(f"Deviation for {protein_name} with {force_field}:")
        print(f"  Alpha Helix: {std_alpha_helix:.2f}%")
        print(f"  Beta Structures: {std_beta_structures:.2f}%")
        print(f"  Other Structures: {std_other_structures:.2f}%")
        print("\n")

if __name__ == "__main__":
    file_path = "/ceph/users/akv16273/34_Elcock_13abundant_inEColi/21_Analysis_Datastructure/26_RbsB_only/1ba2_A_AMBER/Take_1/RbsB_dssp.dat"

    # dssp_matrix = load_dssp_matrix(file_path)
    # plot_dssp_heatmap(dssp_matrix)
    # plot_secondary_structure_percentages(dssp_matrix)

    folder_path = "/ceph/users/akv16273/34_Elcock_13abundant_inEColi/21_Analysis_Datastructure/plots/dssp_results/"
    # plot_multiple_dssp_files_across_takes(folder_path)
    #plot_multiple_dssp_files_across_ForceFields(folder_path)

    # get numbers for text
    calculate_deviation_per_protein(folder_path)