Owens_Expression_data_graph.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "provenance": [],
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/malcolmfisher103/Bioinformatic-Scripts/blob/main/Owens_Expression_data_graph.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "Z02AmrrWqKSF"
      },
      "outputs": [],
      "source": [
        "# @title Specify Genes\n",
        "genes_to_plot = [\"nadsyn1\", \"wnt7a\",\"shh\"] # @param {type:\"raw\"}"
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "# @title Imports\n",
        "import argparse\n",
        "import pandas as pd\n",
        "import seaborn as sns\n",
        "import matplotlib.pyplot as plt\n",
        "import numpy as np  # Add this line\n",
        "from matplotlib.ticker import LogLocator\n",
        "from scipy.interpolate import interp1d"
      ],
      "metadata": {
        "id": "VUrSChfgemcw"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# @title Load data\n",
        "def load_tpm_matrix():\n",
        "    return pd.read_csv('https://bigfrog.xenbase.org/xenbase/genomics/GEO/GSE65785/XENTR_10.0/RNA-Seq/ExpressionFiles/Genes_TPM_Matrix.txt', sep='\\t', index_col=0)\n",
        "\n",
        "def load_gene_info():\n",
        "    return pd.read_csv('https://xenbase-bio1.ucalgary.ca/cgi-bin/reports/models_gene_info.cgi', sep='\\t', index_col=0)\n",
        "\n",
        "def load_gsm_sample_mapping():\n",
        "    return pd.read_csv('https://bigfrog.xenbase.org/xenbase/genomics/GEO/GSE65785/XENTR_10.0/RNA-Seq/gsm_to_track.txt', sep='\\t')\n",
        "\n",
        "def load_gsm_srr_mapping():\n",
        "    return pd.read_csv('https://download.xenbase.org/xenbase/GenePageReports/geo_srr_metadata_chd.txt', sep='\\t')\n",
        "\n",
        "def create_gsm_to_srr_mapping(gsm_srr_mapping):\n",
        "    return gsm_srr_mapping[['SRR', 'GSM']]"
      ],
      "metadata": {
        "id": "mJQZoLF5jIb9"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# @title substitute_gene_symbols\n",
        "def substitute_gene_symbols(tpm_matrix, gene_info):\n",
        "\n",
        "# Drop rows with null values in column 1 and column 2 in file 2\n",
        "  gene_info = gene_info.dropna(subset=['GENE_SYMBOL', 'MODEL_NAME'])\n",
        "\n",
        "# Drop duplicates in file 2 based on column 2\n",
        "  gene_info = gene_info.drop_duplicates(subset=['MODEL_NAME'])\n",
        "\n",
        "# Create a mapping dictionary from column 2 to column 1 in file 2\n",
        "  mapping = dict(zip(gene_info['MODEL_NAME'], gene_info['GENE_SYMBOL']))\n",
        "\n",
        "# Replace values in column 1 of file 1 with corresponding values from file 2\n",
        "  tpm_matrix['Gene'] = tpm_matrix['Gene'].map(mapping).fillna(tpm_matrix['Gene'])\n",
        "\n",
        "# Save the modified DataFrame back to file1.csv or use it as needed\n",
        "#tpm_matrix.to_csv('modified_file1.csv', index=False)\n",
        "  return tpm_matrix"
      ],
      "metadata": {
        "id": "1KrEBM-6wPhC"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# @title process_data\n",
        "def process_data(tpm_matrix, gene_info, gsm_sample_mapping, gsm_srr_mapping):\n",
        "    # Substituting gene symbols\n",
        "    tpm_matrix = substitute_gene_symbols(tpm_matrix, gene_info) #This works\n",
        "\n",
        "    # Mapping GSMs to SRRs\n",
        "    gsm_to_srr = create_gsm_to_srr_mapping(gsm_srr_mapping)\n",
        "\n",
        "    # Select columns 2 and 6 to create a DataFrame with Track Name and GSMs\n",
        "    gsm_mapping = gsm_sample_mapping[['Track Name', 'GSMs']]\n",
        "\n",
        "    # Convert GSMs column to list if it contains comma-separated values\n",
        "    gsm_mapping.loc[:, 'GSMs'] = gsm_mapping['GSMs'].str.split(',')\n",
        "\n",
        "    # Explode the list of GSMs to create multiple rows for each track name\n",
        "    gsm_mapping = gsm_mapping.explode('GSMs')\n",
        "    # Merge gsm_mapping with gsm_srr_mapping on the 'GSM' column to get corresponding SRRs\n",
        "    track_srr_mapping = pd.merge(left=gsm_mapping, right=gsm_srr_mapping, left_on='GSMs', right_on='GSM', validate=\"1:m\")\n",
        "    # Group by 'Track Name' and aggregate the corresponding SRRs into lists\n",
        "    track_srr_mapping = track_srr_mapping.groupby('Track Name')['SRR'].apply(list).reset_index()\n",
        "    track_srr_mapping_expanded = track_srr_mapping.explode('SRR')\n",
        "    #print(track_srr_mapping_expanded.head(10))\n",
        "\n",
        "    srr_columns = tpm_matrix.columns.intersection(track_srr_mapping_expanded['SRR'])\n",
        "    #print(srr_columns)\n",
        "    srr_to_track = dict(zip(track_srr_mapping_expanded['SRR'], track_srr_mapping_expanded['Track Name']))\n",
        "\n",
        "    # Grouping TPM matrix by SRRs\n",
        "    tpm_matrix.rename(columns={col: srr_to_track.get(col, col) for col in srr_columns}, inplace=True)\n",
        "    #print(tpm_matrix.head(10))\n",
        "\n",
        "    # Melt the DataFrame to have a single column for the gene and the rest for values\n",
        "    tpm_matrix = pd.melt(tpm_matrix, id_vars=['Gene'], var_name='Column')\n",
        "\n",
        "    # Extract the unique column names excluding the first column (Gene)\n",
        "    columns_to_merge = tpm_matrix['Column'].unique()[1:]\n",
        "\n",
        "    # Group by track names and take the mean\n",
        "    tpm_matrix_grouped = tpm_matrix.groupby(['Gene', 'Column']).mean().unstack()\n",
        "\n",
        "    # Display the first 10 rows and first 5 columns\n",
        "   # print(tpm_matrix_grouped.iloc[:10, :5])\n",
        "\n",
        "    # Remove the hierarchical index and reset index\n",
        "    tpm_matrix_grouped.columns = tpm_matrix_grouped.columns.droplevel()\n",
        "    tpm_matrix_grouped.reset_index(inplace=True)\n",
        "\n",
        "    return tpm_matrix_grouped"
      ],
      "metadata": {
        "id": "xEL5TjcHwVtL"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# @title Main\n",
        "def main():\n",
        "    # Load data\n",
        "    tpm_matrix = load_tpm_matrix()\n",
        "    tpm_matrix = tpm_matrix.reset_index()\n",
        "    # print(tpm_matrix.columns)\n",
        "    gene_info = load_gene_info()\n",
        "    gene_info = gene_info.reset_index()\n",
        "    # print(gene_info.head(10))\n",
        "    gsm_sample_mapping = load_gsm_sample_mapping()\n",
        "    gsm_srr_mapping = load_gsm_srr_mapping()\n",
        "\n",
        "    # Process data\n",
        "    new_tpm_matrix = process_data(tpm_matrix, gene_info, gsm_sample_mapping, gsm_srr_mapping)\n",
        "\n",
        "    # Output new_tpm_matrix to file\n",
        "    new_tpm_matrix.to_csv('output.txt', sep='\\t')\n",
        "\n",
        "if __name__ == \"__main__\":\n",
        "    main()"
      ],
      "metadata": {
        "id": "9cp1TNvulrKl"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# @title Plot Graph\n",
        "# Load DataFrame\n",
        "new_tpm_matrix = pd.read_csv('output.txt', sep='\\t', index_col=0)\n",
        "\n",
        "df = pd.DataFrame(new_tpm_matrix)\n",
        "df = df[df.columns.drop(list(df.filter(regex='ClutchA rd|ClutchB')))]\n",
        "#print(df.head(1))\n",
        "df = df[['Gene','WE - NF1 (0hpf) - ClutchA','WE - NF1 (0.5hpf) - ClutchA','WE - NF2- (1hpf) - ClutchA','WE - NF3 (1.5hpf) - ClutchA','WE - NF5 (2hpf) - ClutchA','WE - NF6 (2.5hpf) - ClutchA','WE - NF7 (3hpf) - ClutchA','WE - NF7 (3.5hpf) - ClutchA','WE - NF8 (4hpf) - ClutchA','WE - NF9 (4.5hpf) - ClutchA','WE - NF9 (5hpf) - ClutchA','WE - NF9 (5.5hpf) - ClutchA','WE - NF9 (6hpf) - ClutchA','WE - NF10 (6.5hpf) - ClutchA','WE - NF10 (7hpf) - ClutchA','WE - NF10 (7.5hpf) - ClutchA','WE - NF11 (8hpf) - ClutchA','WE - NF12 (8.5hpf) - ClutchA','WE - NF12 (9hpf) - ClutchA','WE - NF12 (9.5hpf) - ClutchA','WE - NF12 (10hpf) - ClutchA','WE - NF12.5 (10.5hpf) - ClutchA','WE - NF12.5 (11hpf) - ClutchA','WE - NF12.5 (11.5hpf) - ClutchA','WE - NF13 (12hpf) - ClutchA','WE - NF13 (12.5hpf) - ClutchA','WE - NF13 (13hpf) - ClutchA','WE - NF14 (13.5hpf) - ClutchA','WE - NF14 (14hpf) - ClutchA','WE - NF15 (14.5hpf) - ClutchA','WE - NF16-17 (15hpf) - ClutchA','WE - NF18-19 (15.5hpf) - ClutchA','WE - NF20 (16hpf) - ClutchA','WE - NF21 (16.5hpf) - ClutchA','WE - NF21 (17hpf) - ClutchA','WE - NF22 (17.5hpf) - ClutchA','WE - NF22 (18hpf) - ClutchA','WE - NF23 (18.5hpf) - ClutchA','WE - NF24 (19hpf) - ClutchA','WE - NF24 (19.5hpf) - ClutchA','WE - NF24 (20hpf) - ClutchA','WE - NF24 (20.5hpf) - ClutchA','WE - NF25 (21hpf) - ClutchA','WE - NF25 (21.5hpf) - ClutchA','WE - NF25 (22hpf) - ClutchA','WE - NF26 (22.5hpf) - ClutchA','WE - NF26 (23hpf) - ClutchA','WE - NF26 (23.5hpf) - ClutchA','WE - NF27 (24hpf) - ClutchA','WE - NF27 (25hpf) - ClutchA','WE - NF27 (26hpf) - ClutchA','WE - NF28 (27hpf) - ClutchA','WE - NF28 (28hpf) - ClutchA','WE - NF29/30 (29hpf) - ClutchA','WE - NF31 (30hpf) - ClutchA','WE - NF32 (31hpf) - ClutchA','WE - NF32 (32hpf) - ClutchA','WE - NF33/34 (33hpf) - ClutchA','WE - NF33/34 (34hpf) - ClutchA','WE - NF33/34 (35hpf) - ClutchA','WE - NF35/36 (36hpf) - ClutchA','WE - NF35/36 (37hpf) - ClutchA','WE - NF37/38 (38hpf) - ClutchA','WE - NF37/38 (39hpf) - ClutchA','WE - NF37/38 (40hpf) - ClutchA','WE - NF39 (41hpf) - ClutchA','WE - NF39 (42hpf) - ClutchA','WE - NF39 (44hpf) - ClutchA','WE - NF39 (45hpf) - ClutchA','WE - NF40 (46hpf) - ClutchA','WE - NF40 (47hpf) - ClutchA','WE - NF40 (48hpf) - ClutchA','WE - NF40 (49hpf) - ClutchA','WE - NF40 (50hpf) - ClutchA','WE - NF41 (51hpf) - ClutchA','WE - NF41 (52hpf) - ClutchA','WE - NF41 (53hpf) - ClutchA','WE - NF41 (54hpf) - ClutchA','WE - NF41 (55hpf) - ClutchA','WE - NF42 (56hpf) - ClutchA','WE - NF42 (57hpf) - ClutchA','WE - NF42 (58hpf) - ClutchA','WE - NF42 (59hpf) - ClutchA','WE - NF42 (60hpf) - ClutchA','WE - NF42 (61hpf) - ClutchA','WE - NF42 (62hpf) - ClutchA','WE - NF43 (63hpf) - ClutchA','WE - NF43 (64hpf) - ClutchA','WE - NF43 (65hpf) - ClutchA','WE - NF43 (66hpf) - ClutchA']]\n",
        "#print(df.head(10))\n",
        "# List of genes to include in the plot\n",
        "#genes_to_plot = genes_to_plot[\"pax2.L\", \"nadsyn1.L\", \"slc7a5.L\"]\n",
        "\n",
        "# Filter the DataFrame to include only the specified genes\n",
        "filtered_df = df[df['Gene'].isin(genes_to_plot)]\n",
        "#print(filtered_df.head(10))\n",
        "\n",
        "# Plotting\n",
        "plt.figure(figsize=(10, 6))\n",
        "\n",
        "for index, row in filtered_df.iterrows():\n",
        "    gene = row['Gene']\n",
        "    expression_values = row.drop('Gene').values\n",
        "#    time_points = row.drop('Gene').index\n",
        "    #time_points =['NF stage 1','NF stage 2 (2-cell)','NF stage 4 (8-cell)','NF stage 5 (16-cell)','NF stage 6 (32-cell)','NF stage 7','NF stage 8','NF stage 9','NF stage 10','NF stage 11','NF stage 12','NF stage 12.5','NF stage 13','NF stage 14','NF stage 15','NF stage 17','NF stage 19','NF stage 20','NF stage 21','NF stage 22','NF stage 23','NF stage 24','NF stage 25','NF stage 26','NF stage 27','NF stage 28','NF stage 29 and 30','NF stage 31','NF stage 32','NF stage 33 and 34','NF stage 35 and 36','NF stage 37 and 38','NF stage 39','NF stage 40','NF stage 41','NF stage 42']\n",
        "    time_points =['WE - NF1 (0hpf) - ClutchA','WE - NF1 (0.5hpf) - ClutchA','WE - NF2- (1hpf) - ClutchA','WE - NF3 (1.5hpf) - ClutchA','WE - NF5 (2hpf) - ClutchA','WE - NF6 (2.5hpf) - ClutchA','WE - NF7 (3hpf) - ClutchA','WE - NF7 (3.5hpf) - ClutchA','WE - NF8 (4hpf) - ClutchA','WE - NF9 (4.5hpf) - ClutchA','WE - NF9 (5hpf) - ClutchA','WE - NF9 (5.5hpf) - ClutchA','WE - NF9 (6hpf) - ClutchA','WE - NF10 (6.5hpf) - ClutchA','WE - NF10 (7hpf) - ClutchA','WE - NF10 (7.5hpf) - ClutchA','WE - NF11 (8hpf) - ClutchA','WE - NF12 (8.5hpf) - ClutchA','WE - NF12 (9hpf) - ClutchA','WE - NF12 (9.5hpf) - ClutchA','WE - NF12 (10hpf) - ClutchA','WE - NF12.5 (10.5hpf) - ClutchA','WE - NF12.5 (11hpf) - ClutchA','WE - NF12.5 (11.5hpf) - ClutchA','WE - NF13 (12hpf) - ClutchA','WE - NF13 (12.5hpf) - ClutchA','WE - NF13 (13hpf) - ClutchA','WE - NF14 (13.5hpf) - ClutchA','WE - NF14 (14hpf) - ClutchA','WE - NF15 (14.5hpf) - ClutchA','WE - NF16-17 (15hpf) - ClutchA','WE - NF18-19 (15.5hpf) - ClutchA','WE - NF20 (16hpf) - ClutchA','WE - NF21 (16.5hpf) - ClutchA','WE - NF21 (17hpf) - ClutchA','WE - NF22 (17.5hpf) - ClutchA','WE - NF22 (18hpf) - ClutchA','WE - NF23 (18.5hpf) - ClutchA','WE - NF24 (19hpf) - ClutchA','WE - NF24 (19.5hpf) - ClutchA','WE - NF24 (20hpf) - ClutchA','WE - NF24 (20.5hpf) - ClutchA','WE - NF25 (21hpf) - ClutchA','WE - NF25 (21.5hpf) - ClutchA','WE - NF25 (22hpf) - ClutchA','WE - NF26 (22.5hpf) - ClutchA','WE - NF26 (23hpf) - ClutchA','WE - NF26 (23.5hpf) - ClutchA','WE - NF27 (24hpf) - ClutchA','WE - NF27 (25hpf) - ClutchA','WE - NF27 (26hpf) - ClutchA','WE - NF28 (27hpf) - ClutchA','WE - NF28 (28hpf) - ClutchA','WE - NF29/30 (29hpf) - ClutchA','WE - NF31 (30hpf) - ClutchA','WE - NF32 (31hpf) - ClutchA','WE - NF32 (32hpf) - ClutchA','WE - NF33/34 (33hpf) - ClutchA','WE - NF33/34 (34hpf) - ClutchA','WE - NF33/34 (35hpf) - ClutchA','WE - NF35/36 (36hpf) - ClutchA','WE - NF35/36 (37hpf) - ClutchA','WE - NF37/38 (38hpf) - ClutchA','WE - NF37/38 (39hpf) - ClutchA','WE - NF37/38 (40hpf) - ClutchA','WE - NF39 (41hpf) - ClutchA','WE - NF39 (42hpf) - ClutchA','WE - NF39 (44hpf) - ClutchA','WE - NF39 (45hpf) - ClutchA','WE - NF40 (46hpf) - ClutchA','WE - NF40 (47hpf) - ClutchA','WE - NF40 (48hpf) - ClutchA','WE - NF40 (49hpf) - ClutchA','WE - NF40 (50hpf) - ClutchA','WE - NF41 (51hpf) - ClutchA','WE - NF41 (52hpf) - ClutchA','WE - NF41 (53hpf) - ClutchA','WE - NF41 (54hpf) - ClutchA','WE - NF41 (55hpf) - ClutchA','WE - NF42 (56hpf) - ClutchA','WE - NF42 (57hpf) - ClutchA','WE - NF42 (58hpf) - ClutchA','WE - NF42 (59hpf) - ClutchA','WE - NF42 (60hpf) - ClutchA','WE - NF42 (61hpf) - ClutchA','WE - NF42 (62hpf) - ClutchA','WE - NF43 (63hpf) - ClutchA','WE - NF43 (64hpf) - ClutchA','WE - NF43 (65hpf) - ClutchA','WE - NF43 (66hpf) - ClutchA']\n",
        "# Perform cubic spline interpolation\n",
        "    f = interp1d(range(len(time_points)), expression_values, kind='cubic')\n",
        "    x_new = np.linspace(0, len(time_points) - 1, 100)  # Adjust 100 for desired smoothness\n",
        "    y_new = np.maximum(f(x_new), 0)\n",
        "    plt.plot(time_points, expression_values, 'o', label=gene)  # Plot original data points\n",
        "    plt.plot(x_new, y_new, label=f'{gene} (smoothed)')\n",
        "\n",
        "# Customize x-axis ticks\n",
        "tick_positions = range(len(time_points))\n",
        "tick_labels = [time_points[i] if i % 4 == 0 else '' for i in range(len(time_points))]  # Skip every other label\n",
        "plt.xticks(tick_positions, tick_labels, rotation=45, ha='right')\n",
        "\n",
        "plt.xlabel('Time Points')\n",
        "plt.ylabel('Gene Expression Level (TPM)')\n",
        "plt.title('Gene Expression Data')\n",
        "#plt.xticks(rotation=45, ha='right')\n",
        "plt.legend()\n",
        "plt.grid(True)\n",
        "#plt.yscale('log')  # Set y-axis to logarithmic scale\n",
        "#plt.gca().set_yscale('function', functions=(lambda x: np.log2(x), lambda x: np.exp2(x))) # Set y-axis to logarithmic scale with base 2\n",
        "#plt.gca().yaxis.set_major_locator(LogLocator(base=2.0)) # Use LogLocator for y-axis ticks\n",
        "plt.tight_layout()\n",
        "plt.show()"
      ],
      "metadata": {
        "id": "CmTCEDytecp2"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# @title Graph with proper timing\n",
        "# Load DataFrame\n",
        "new_tpm_matrix = pd.read_csv('output.txt', sep='\\t', index_col=0)\n",
        "# Given timings\n",
        "timings = {\n",
        "\t'WE - NF1 (0hpf) - ClutchA': 0,\n",
        "\t'WE - NF1 (0.5hpf) - ClutchA': 0.5,\n",
        "\t'WE - NF2- (1hpf) - ClutchA': 1,\n",
        "\t'WE - NF3 (1.5hpf) - ClutchA': 1.5,\n",
        "\t'WE - NF5 (2hpf) - ClutchA': 2,\n",
        "\t'WE - NF6 (2.5hpf) - ClutchA': 2.5,\n",
        "\t'WE - NF7 (3hpf) - ClutchA': 3,\n",
        "\t'WE - NF7 (3.5hpf) - ClutchA': 3.5 ,\n",
        "\t'WE - NF8 (4hpf) - ClutchA': 4,\n",
        "\t'WE - NF9 (4.5hpf) - ClutchA': 4.5,\n",
        "\t'WE - NF9 (5hpf) - ClutchA': 5,\n",
        "\t'WE - NF9 (5.5hpf) - ClutchA': 5.5,\n",
        "\t'WE - NF9 (6hpf) - ClutchA': 6,\n",
        "\t'WE - NF10 (6.5hpf) - ClutchA': 6.5,\n",
        "\t'WE - NF10 (7hpf) - ClutchA': 7,\n",
        "\t'WE - NF10 (7.5hpf) - ClutchA': 7.5,\n",
        "\t'WE - NF11 (8hpf) - ClutchA': 8,\n",
        "\t'WE - NF12 (8.5hpf) - ClutchA': 8.5,\n",
        "\t'WE - NF12 (9hpf) - ClutchA': 9,\n",
        "\t'WE - NF12 (9.5hpf) - ClutchA': 9.5,\n",
        "\t'WE - NF12 (10hpf) - ClutchA': 10,\n",
        "\t'WE - NF12.5 (10.5hpf) - ClutchA': 10.5,\n",
        "\t'WE - NF12.5 (11hpf) - ClutchA': 11,\n",
        "\t'WE - NF12.5 (11.5hpf) - ClutchA': 11.5,\n",
        "\t'WE - NF13 (12hpf) - ClutchA': 12,\n",
        "\t'WE - NF13 (12.5hpf) - ClutchA': 12.5,\n",
        "\t'WE - NF13 (13hpf) - ClutchA': 13,\n",
        "\t'WE - NF14 (13.5hpf) - ClutchA': 13.5,\n",
        "\t'WE - NF14 (14hpf) - ClutchA': 14,\n",
        "\t'WE - NF15 (14.5hpf) - ClutchA': 14.5,\n",
        "\t'WE - NF16-17 (15hpf) - ClutchA': 15,\n",
        "\t'WE - NF18-19 (15.5hpf) - ClutchA': 15.5,\n",
        "\t'WE - NF20 (16hpf) - ClutchA': 16,\n",
        "\t'WE - NF21 (16.5hpf) - ClutchA': 16.5,\n",
        "\t'WE - NF21 (17hpf) - ClutchA': 17,\n",
        "\t'WE - NF22 (17.5hpf) - ClutchA': 17.5,\n",
        "\t'WE - NF22 (18hpf) - ClutchA': 18,\n",
        "\t'WE - NF23 (18.5hpf) - ClutchA': 18.5,\n",
        "\t'WE - NF24 (19hpf) - ClutchA': 19,\n",
        "\t'WE - NF24 (19.5hpf) - ClutchA': 19.5,\n",
        "\t'WE - NF24 (20hpf) - ClutchA': 20,\n",
        "\t'WE - NF24 (20.5hpf) - ClutchA': 20.5,\n",
        "\t'WE - NF25 (21hpf) - ClutchA': 21,\n",
        "\t'WE - NF25 (21.5hpf) - ClutchA': 21.5,\n",
        "\t'WE - NF25 (22hpf) - ClutchA': 22,\n",
        "\t'WE - NF26 (22.5hpf) - ClutchA': 22.5,\n",
        "\t'WE - NF26 (23hpf) - ClutchA': 23,\n",
        "\t'WE - NF26 (23.5hpf) - ClutchA': 23.5,\n",
        "\t'WE - NF27 (24hpf) - ClutchA': 24,\n",
        "\t'WE - NF27 (25hpf) - ClutchA': 25,\n",
        "\t'WE - NF27 (26hpf) - ClutchA': 26,\n",
        "\t'WE - NF28 (27hpf) - ClutchA': 27,\n",
        "\t'WE - NF28 (28hpf) - ClutchA': 28,\n",
        "\t'WE - NF29/30 (29hpf) - ClutchA': 29,\n",
        "\t'WE - NF31 (30hpf) - ClutchA': 30,\n",
        "\t'WE - NF32 (31hpf) - ClutchA': 31,\n",
        "\t'WE - NF32 (32hpf) - ClutchA': 32,\n",
        "\t'WE - NF33/34 (33hpf) - ClutchA': 33,\n",
        "\t'WE - NF33/34 (34hpf) - ClutchA': 34,\n",
        "\t'WE - NF33/34 (35hpf) - ClutchA': 35,\n",
        "\t'WE - NF35/36 (36hpf) - ClutchA': 36,\n",
        "\t'WE - NF35/36 (37hpf) - ClutchA': 37,\n",
        "\t'WE - NF37/38 (38hpf) - ClutchA': 38,\n",
        "\t'WE - NF37/38 (39hpf) - ClutchA': 39,\n",
        "\t'WE - NF37/38 (40hpf) - ClutchA': 40,\n",
        "\t'WE - NF39 (41hpf) - ClutchA': 41,\n",
        "\t'WE - NF39 (42hpf) - ClutchA': 42,\n",
        "\t'WE - NF39 (44hpf) - ClutchA': 44,\n",
        "\t'WE - NF39 (45hpf) - ClutchA': 45,\n",
        "\t'WE - NF40 (46hpf) - ClutchA': 46,\n",
        "\t'WE - NF40 (47hpf) - ClutchA': 47,\n",
        "\t'WE - NF40 (48hpf) - ClutchA': 48,\n",
        "\t'WE - NF40 (49hpf) - ClutchA': 49,\n",
        "\t'WE - NF40 (50hpf) - ClutchA': 50,\n",
        "\t'WE - NF41 (51hpf) - ClutchA': 51,\n",
        "\t'WE - NF41 (52hpf) - ClutchA': 52,\n",
        "\t'WE - NF41 (53hpf) - ClutchA': 53,\n",
        "\t'WE - NF41 (54hpf) - ClutchA': 54,\n",
        "\t'WE - NF41 (55hpf) - ClutchA': 55,\n",
        "\t'WE - NF42 (56hpf) - ClutchA': 56,\n",
        "\t'WE - NF42 (57hpf) - ClutchA': 57,\n",
        "\t'WE - NF42 (58hpf) - ClutchA': 58,\n",
        "\t'WE - NF42 (59hpf) - ClutchA': 59,\n",
        "\t'WE - NF42 (60hpf) - ClutchA': 60,\n",
        "\t'WE - NF42 (61hpf) - ClutchA': 61,\n",
        "\t'WE - NF42 (62hpf) - ClutchA': 62,\n",
        "\t'WE - NF43 (63hpf) - ClutchA': 63,\n",
        "\t'WE - NF43 (64hpf) - ClutchA': 64,\n",
        "\t'WE - NF43 (65hpf) - ClutchA': 65,\n",
        "\t'WE - NF43 (66hpf) - ClutchA': 66\n",
        "}\n",
        "\n",
        "# Create DataFrame\n",
        "df = pd.DataFrame(new_tpm_matrix)\n",
        "df = df[['Gene', 'WE - NF1 (0hpf) - ClutchA','WE - NF1 (0.5hpf) - ClutchA','WE - NF2- (1hpf) - ClutchA','WE - NF3 (1.5hpf) - ClutchA','WE - NF5 (2hpf) - ClutchA','WE - NF6 (2.5hpf) - ClutchA','WE - NF7 (3hpf) - ClutchA','WE - NF7 (3.5hpf) - ClutchA','WE - NF8 (4hpf) - ClutchA','WE - NF9 (4.5hpf) - ClutchA','WE - NF9 (5hpf) - ClutchA','WE - NF9 (5.5hpf) - ClutchA','WE - NF9 (6hpf) - ClutchA','WE - NF10 (6.5hpf) - ClutchA','WE - NF10 (7hpf) - ClutchA','WE - NF10 (7.5hpf) - ClutchA','WE - NF11 (8hpf) - ClutchA','WE - NF12 (8.5hpf) - ClutchA','WE - NF12 (9hpf) - ClutchA','WE - NF12 (9.5hpf) - ClutchA','WE - NF12 (10hpf) - ClutchA','WE - NF12.5 (10.5hpf) - ClutchA','WE - NF12.5 (11hpf) - ClutchA','WE - NF12.5 (11.5hpf) - ClutchA','WE - NF13 (12hpf) - ClutchA','WE - NF13 (12.5hpf) - ClutchA','WE - NF13 (13hpf) - ClutchA','WE - NF14 (13.5hpf) - ClutchA','WE - NF14 (14hpf) - ClutchA','WE - NF15 (14.5hpf) - ClutchA','WE - NF16-17 (15hpf) - ClutchA','WE - NF18-19 (15.5hpf) - ClutchA','WE - NF20 (16hpf) - ClutchA','WE - NF21 (16.5hpf) - ClutchA','WE - NF21 (17hpf) - ClutchA','WE - NF22 (17.5hpf) - ClutchA','WE - NF22 (18hpf) - ClutchA','WE - NF23 (18.5hpf) - ClutchA','WE - NF24 (19hpf) - ClutchA','WE - NF24 (19.5hpf) - ClutchA','WE - NF24 (20hpf) - ClutchA','WE - NF24 (20.5hpf) - ClutchA','WE - NF25 (21hpf) - ClutchA','WE - NF25 (21.5hpf) - ClutchA','WE - NF25 (22hpf) - ClutchA','WE - NF26 (22.5hpf) - ClutchA','WE - NF26 (23hpf) - ClutchA','WE - NF26 (23.5hpf) - ClutchA','WE - NF27 (24hpf) - ClutchA','WE - NF27 (25hpf) - ClutchA','WE - NF27 (26hpf) - ClutchA','WE - NF28 (27hpf) - ClutchA','WE - NF28 (28hpf) - ClutchA','WE - NF29/30 (29hpf) - ClutchA','WE - NF31 (30hpf) - ClutchA','WE - NF32 (31hpf) - ClutchA','WE - NF32 (32hpf) - ClutchA','WE - NF33/34 (33hpf) - ClutchA','WE - NF33/34 (34hpf) - ClutchA','WE - NF33/34 (35hpf) - ClutchA','WE - NF35/36 (36hpf) - ClutchA','WE - NF35/36 (37hpf) - ClutchA','WE - NF37/38 (38hpf) - ClutchA','WE - NF37/38 (39hpf) - ClutchA','WE - NF37/38 (40hpf) - ClutchA','WE - NF39 (41hpf) - ClutchA','WE - NF39 (42hpf) - ClutchA','WE - NF39 (44hpf) - ClutchA','WE - NF39 (45hpf) - ClutchA','WE - NF40 (46hpf) - ClutchA','WE - NF40 (47hpf) - ClutchA','WE - NF40 (48hpf) - ClutchA','WE - NF40 (49hpf) - ClutchA','WE - NF40 (50hpf) - ClutchA','WE - NF41 (51hpf) - ClutchA','WE - NF41 (52hpf) - ClutchA','WE - NF41 (53hpf) - ClutchA','WE - NF41 (54hpf) - ClutchA','WE - NF41 (55hpf) - ClutchA','WE - NF42 (56hpf) - ClutchA','WE - NF42 (57hpf) - ClutchA','WE - NF42 (58hpf) - ClutchA','WE - NF42 (59hpf) - ClutchA','WE - NF42 (60hpf) - ClutchA','WE - NF42 (61hpf) - ClutchA','WE - NF42 (62hpf) - ClutchA','WE - NF43 (63hpf) - ClutchA','WE - NF43 (64hpf) - ClutchA','WE - NF43 (65hpf) - ClutchA','WE - NF43 (66hpf) - ClutchA']]\n",
        "\n",
        "# Filter the DataFrame to include only the specified genes\n",
        "filtered_df = df[df['Gene'].isin(genes_to_plot)]\n",
        "\n",
        "# Plotting\n",
        "plt.figure(figsize=(10, 6))\n",
        "\n",
        "for index, row in filtered_df.iterrows():\n",
        "    gene = row['Gene']\n",
        "    expression_values = row.drop('Gene').values\n",
        "    time_points = row.drop('Gene').index\n",
        "\n",
        "    # Interpolating time points to align with the given timings\n",
        "    stage_positions = [timings[stage] for stage in time_points]\n",
        "\n",
        "    # Perform cubic spline interpolation\n",
        "    f = interp1d(stage_positions, expression_values, kind='cubic')\n",
        "    x_new = np.linspace(min(stage_positions), max(stage_positions), 100)  # Adjust 100 for desired smoothness\n",
        "    y_new = np.maximum(f(x_new), 0)\n",
        "\n",
        "    stage_labels =  ['WE - NF1 (0hpf) - ClutchA','WE - NF1 (0.5hpf) - ClutchA','WE - NF2- (1hpf) - ClutchA','WE - NF3 (1.5hpf) - ClutchA','WE - NF5 (2hpf) - ClutchA','WE - NF6 (2.5hpf) - ClutchA','WE - NF7 (3hpf) - ClutchA','WE - NF7 (3.5hpf) - ClutchA','WE - NF8 (4hpf) - ClutchA','WE - NF9 (4.5hpf) - ClutchA','WE - NF9 (5hpf) - ClutchA','WE - NF9 (5.5hpf) - ClutchA','WE - NF9 (6hpf) - ClutchA','WE - NF10 (6.5hpf) - ClutchA','WE - NF10 (7hpf) - ClutchA','WE - NF10 (7.5hpf) - ClutchA','WE - NF11 (8hpf) - ClutchA','WE - NF12 (8.5hpf) - ClutchA','WE - NF12 (9hpf) - ClutchA','WE - NF12 (9.5hpf) - ClutchA','WE - NF12 (10hpf) - ClutchA','WE - NF12.5 (10.5hpf) - ClutchA','WE - NF12.5 (11hpf) - ClutchA','WE - NF12.5 (11.5hpf) - ClutchA','WE - NF13 (12hpf) - ClutchA','WE - NF13 (12.5hpf) - ClutchA','WE - NF13 (13hpf) - ClutchA','WE - NF14 (13.5hpf) - ClutchA','WE - NF14 (14hpf) - ClutchA','WE - NF15 (14.5hpf) - ClutchA','WE - NF16-17 (15hpf) - ClutchA','WE - NF18-19 (15.5hpf) - ClutchA','WE - NF20 (16hpf) - ClutchA','WE - NF21 (16.5hpf) - ClutchA','WE - NF21 (17hpf) - ClutchA','WE - NF22 (17.5hpf) - ClutchA','WE - NF22 (18hpf) - ClutchA','WE - NF23 (18.5hpf) - ClutchA','WE - NF24 (19hpf) - ClutchA','WE - NF24 (19.5hpf) - ClutchA','WE - NF24 (20hpf) - ClutchA','WE - NF24 (20.5hpf) - ClutchA','WE - NF25 (21hpf) - ClutchA','WE - NF25 (21.5hpf) - ClutchA','WE - NF25 (22hpf) - ClutchA','WE - NF26 (22.5hpf) - ClutchA','WE - NF26 (23hpf) - ClutchA','WE - NF26 (23.5hpf) - ClutchA','WE - NF27 (24hpf) - ClutchA','WE - NF27 (25hpf) - ClutchA','WE - NF27 (26hpf) - ClutchA','WE - NF28 (27hpf) - ClutchA','WE - NF28 (28hpf) - ClutchA','WE - NF29/30 (29hpf) - ClutchA','WE - NF31 (30hpf) - ClutchA','WE - NF32 (31hpf) - ClutchA','WE - NF32 (32hpf) - ClutchA','WE - NF33/34 (33hpf) - ClutchA','WE - NF33/34 (34hpf) - ClutchA','WE - NF33/34 (35hpf) - ClutchA','WE - NF35/36 (36hpf) - ClutchA','WE - NF35/36 (37hpf) - ClutchA','WE - NF37/38 (38hpf) - ClutchA','WE - NF37/38 (39hpf) - ClutchA','WE - NF37/38 (40hpf) - ClutchA','WE - NF39 (41hpf) - ClutchA','WE - NF39 (42hpf) - ClutchA','WE - NF39 (44hpf) - ClutchA','WE - NF39 (45hpf) - ClutchA','WE - NF40 (46hpf) - ClutchA','WE - NF40 (47hpf) - ClutchA','WE - NF40 (48hpf) - ClutchA','WE - NF40 (49hpf) - ClutchA','WE - NF40 (50hpf) - ClutchA','WE - NF41 (51hpf) - ClutchA','WE - NF41 (52hpf) - ClutchA','WE - NF41 (53hpf) - ClutchA','WE - NF41 (54hpf) - ClutchA','WE - NF41 (55hpf) - ClutchA','WE - NF42 (56hpf) - ClutchA','WE - NF42 (57hpf) - ClutchA','WE - NF42 (58hpf) - ClutchA','WE - NF42 (59hpf) - ClutchA','WE - NF42 (60hpf) - ClutchA','WE - NF42 (61hpf) - ClutchA','WE - NF42 (62hpf) - ClutchA','WE - NF43 (63hpf) - ClutchA','WE - NF43 (64hpf) - ClutchA','WE - NF43 (65hpf) - ClutchA','WE - NF43 (66hpf) - ClutchA']\n",
        "\n",
        "    plt.plot(stage_positions, expression_values, 'o', label=gene)  # Plot original data points\n",
        "    plt.plot(x_new, y_new, label=f'{gene} (smoothed)')  # Plot smoothed line\n",
        "\n",
        "# Customize x-axis ticks\n",
        "tick_positions = range(len(stage_labels))\n",
        "tick_labels = [stage_labels[i] if i % 4 == 0 else '' for i in range(len(stage_labels))]  # Skip every other label\n",
        "\n",
        "plt.xlabel('NF stage')  # X-axis represents time in hours\n",
        "plt.ylabel('Gene Expression Level (TPM)')\n",
        "plt.title('Gene Expression Data')\n",
        "plt.xticks(list(timings.values()), tick_labels, rotation=45, ha='right')  # Set custom ticks and labels\n",
        "plt.legend()\n",
        "plt.grid(True)\n",
        "plt.tight_layout()\n",
        "plt.show()\n"
      ],
      "metadata": {
        "id": "82OrIp-nOkeK"
      },
      "execution_count": null,
      "outputs": []
    }
  ]
}