analysis.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "initial_id",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "\n",
    "file_path = 'data/Supplemental_Data_Raw_Genecounts.csv'\n",
    "\n",
    "# Load the CSV file into a DataFrame.\n",
    "df = pd.read_csv(file_path)\n",
    "df = df.drop('Unnamed: 17', axis=1)\n",
    "df = df.rename(columns={\"Unnamed: 0\": \"mnra\",})\n",
    "df.set_index('mnra', inplace=True)\n",
    "# Convert every non-index column to int\n",
    "df = df.apply(pd.to_numeric, errors='coerce', downcast='integer')\n",
    "\n",
    "# Pruning rows where all the values are 0.\n",
    "# Assuming that you have columns 'x' and 'y', replace them with the actual names of your columns.\n",
    "print(df.shape)\n",
    "df = df[(df.loc[:, df.columns != 'index'] != 0).any(axis=1)]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6b5dd579777c0d6d",
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "# Pruning rows where all the values in columns_to_check are 0.\n",
    "df = df[(df != 0).any(axis=1)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "777c7d51c5dd2033",
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "df.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1fd64f5081067189",
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "df.columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "13c096b4-b098-4ffe-beb7-51e623842e43",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Step 1: Calculate pairwise correlation\n",
    "correlation_matrix = df.corr(method='pearson')\n",
    "\n",
    "import seaborn as sns\n",
    "import matplotlib.pyplot as plt\n",
    "sns.heatmap(correlation_matrix, annot=False, cmap='viridis')\n",
    "plt.title('sample-sample Correlation Heatmap')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3d4f80fb0bb91c4c",
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Calculating the mean and standard deviation for each row\n",
    "df['mean'] = df.mean(axis=1)\n",
    "df['std'] = df.std(axis=1)\n",
    "df['cv'] = (df.std(axis=1) / df.mean(axis=1)).abs() * 100\n",
    "\n",
    "import numpy as np\n",
    "\n",
    "# Generate a color map\n",
    "colors = plt.cm.viridis(np.linspace(0, 1, len(df)))\n",
    "\n",
    "plt.figure(figsize=(10, 6))\n",
    "for idx, color in zip(df.index, colors):\n",
    "    plt.errorbar(idx, df.loc[idx, 'mean'], yerr=df.loc[idx, 'std'], fmt='o', color=color)\n",
    "\n",
    "plt.title('Mean and Variability of Each Sample')\n",
    "plt.xlabel('Sample Index')\n",
    "plt.ylabel('Mean')\n",
    "plt.grid(True)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "52db739acf122215",
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "# Plotting histogram for df['cv']\n",
    "plt.figure(figsize=(10,6))\n",
    "plt.hist(df['cv'], bins=200, color='blue', edgecolor='black')\n",
    "plt.title('Histogram of Coefficient of Variation (CV)')\n",
    "plt.xlabel('CV')\n",
    "plt.ylabel('Frequency')\n",
    "plt.grid(True)\n",
    "plt.show()\n",
    "\n",
    "# Plotting histogram for df['mean']\n",
    "plt.figure(figsize=(10,6))\n",
    "plt.hist(df['mean'], bins=200, color='green', edgecolor='black')\n",
    "plt.title('Histogram of Mean')\n",
    "plt.xlabel('Mean')\n",
    "plt.ylabel('Frequency')\n",
    "plt.grid(True)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bb7a2b5bfba33a85",
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "33d4164c559ba84",
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "source": [
    "## Differential expression analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bbfdc95b-0115-498a-8995-e1dbf01eb55d",
   "metadata": {
    "jp-MarkdownHeadingCollapsed": true
   },
   "source": [
    "### Pairwise"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1fcd8ff4e2b9ede1",
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "from statsmodels.stats import multitest\n",
    "from scipy import stats\n",
    "import itertools\n",
    "from joblib import Parallel, delayed\n",
    "from tqdm import tqdm\n",
    "\n",
    "# Define the function to perform the t-test for a given pair of conditions and a given gene\n",
    "def perform_ttest(condition_1, condition_2, gene_row):\n",
    "    condition_1_data = gene_row[condition_1]\n",
    "    condition_2_data = gene_row[condition_2]\n",
    "    \n",
    "    # If the variance is zero for both conditions, skip the t-test\n",
    "    if np.var([condition_1_data]) == 0 and np.var([condition_2_data]) == 0:\n",
    "        return {\n",
    "            'condition_1': condition_1,\n",
    "            'condition_2': condition_2,\n",
    "            'gene': gene_row.name,\n",
    "            'mean_condition_1': condition_1_data,\n",
    "            'mean_condition_2': condition_2_data,\n",
    "            'p_value': np.nan  # or return some default p_value, like 1.0\n",
    "        }\n",
    "    \n",
    "    # Compute the mean expression level for each condition\n",
    "    mean_condition_1 = np.mean(condition_1_data)\n",
    "    mean_condition_2 = np.mean(condition_2_data)\n",
    "    \n",
    "    # Perform a t-test between the two conditions\n",
    "    t_stat, p_value = stats.ttest_ind([condition_1_data], [condition_2_data], equal_var=False, nan_policy='omit')\n",
    "    \n",
    "    return {\n",
    "        'condition_1': condition_1,\n",
    "        'condition_2': condition_2,\n",
    "        'gene': gene_row.name,\n",
    "        'mean_condition_1': mean_condition_1,\n",
    "        'mean_condition_2': mean_condition_2,\n",
    "        'p_value': p_value\n",
    "    }\n",
    "\n",
    "# Get all condition columns\n",
    "condition_columns = df.columns\n",
    "\n",
    "\n",
    "# Create a list to store the tasks\n",
    "tasks = []\n",
    "\n",
    "# Initialize the progress bar\n",
    "total_combinations = len(list(itertools.combinations(condition_columns, 2))) * len(df.index)\n",
    "progress_bar = tqdm(total=total_combinations, desc='Generating tasks', unit='task')\n",
    "\n",
    "# Generate all tasks (pairs of conditions and genes) with progress update\n",
    "for condition_1, condition_2 in itertools.combinations(condition_columns, 2):\n",
    "    for gene in df.index:\n",
    "        tasks.append((condition_1, condition_2, df.loc[gene]))\n",
    "        progress_bar.update(1)\n",
    "\n",
    "# Close the progress bar\n",
    "progress_bar.close()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b1e7c25e31952f5b",
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "\n",
    "# Perform the tasks in parallel with progress bar\n",
    "results_list = Parallel(n_jobs=-1)(delayed(perform_ttest)(condition_1, condition_2, gene_row) for condition_1, condition_2, gene_row in tqdm(tasks, desc='Processing', unit='task'))\n",
    "\n",
    "# Convert the results list to a DataFrame\n",
    "results_df = pd.DataFrame(results_list)\n",
    "\n",
    "# Correct for multiple testing (optional)\n",
    "results_df['adjusted_p_value'] = multitest.multipletests(results_df['p_value'], method='fdr_bh')[1]\n",
    "\n",
    "# You can filter the results based on a significance threshold, e.g. adjusted p-value < 0.05\n",
    "significant_results = results_df[results_df['adjusted_p_value'] < 0.05]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6dbf6d237d5cf343",
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "results_df.columns"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e7ac0ebf5eb36a03",
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "source": [
    "### A and B"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "39b7b9388a237add",
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "from scipy import stats\n",
    "import numpy as np\n",
    "from statsmodels.stats import multitest\n",
    "\n",
    "\n",
    "# Define the groups you want to compare\n",
    "group_a_cols = [col for col in df.columns if 'A' in col]\n",
    "group_b_cols = [col for col in df.columns if 'B' in col]\n",
    "\n",
    "# Create an empty DataFrame to store the results\n",
    "results_df = pd.DataFrame(columns=['gene', 'mean_group_a', 'mean_group_b', 'p_value'])\n",
    "\n",
    "for index, row in df.iterrows():\n",
    "    group_a_data = row[group_a_cols]\n",
    "    group_b_data = row[group_b_cols]\n",
    "    \n",
    "    # Compute the mean expression level for each group\n",
    "    mean_group_a = np.mean(group_a_data)\n",
    "    mean_group_b = np.mean(group_b_data)\n",
    "    \n",
    "    # Perform a t-test between the two groups\n",
    "    t_stat, p_value = stats.ttest_ind(group_a_data, group_b_data, equal_var=False, nan_policy='omit')\n",
    "    \n",
    "    # Append the results to the results DataFrame\n",
    "    new_row = pd.DataFrame({\n",
    "        'gene': [index],\n",
    "        'mean_group_a': [mean_group_a],\n",
    "        'mean_group_b': [mean_group_b],\n",
    "        'p_value': [p_value]\n",
    "    })\n",
    "    \n",
    "    results_df = pd.concat([results_df, new_row], ignore_index=True)\n",
    "\n",
    "\n",
    "# Correct for multiple testing (optional)\n",
    "results_df['adjusted_p_value'] = multitest.multipletests(results_df['p_value'], method='fdr_bh')[1]\n",
    "results_df.set_index('gene', inplace=True)\n",
    "\n",
    "# Define significance threshold\n",
    "sig_threshold = 0.05\n",
    "\n",
    "# Filter the results based on a significance threshold, for example, adjusted p-value < 0.05\n",
    "significant_genes = results_df[results_df['adjusted_p_value'] < sig_threshold]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "43aaca3d0c5afb53",
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "significant_genes.columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "97d69469e9647d90",
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "# Compute log2 fold change\n",
    "results_df['log2_fold_change'] = np.log2(results_df['mean_group_b'] / results_df['mean_group_a'])\n",
    "\n",
    "# Compute -log10 of the adjusted p-value\n",
    "results_df['neg_log10_p_value'] = -np.log10(results_df['adjusted_p_value'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3ee7f28cdd33e2d",
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "\n",
    "# Create a figure and a set of subplots.\n",
    "fig, ax = plt.subplots(figsize=(10, 6))\n",
    "\n",
    "# Color significant points differently\n",
    "colors = np.where(results_df['adjusted_p_value'] < sig_threshold, 'r', 'b')\n",
    "\n",
    "# Scatter plot\n",
    "ax.scatter(results_df['log2_fold_change'], results_df['neg_log10_p_value'], c=colors, alpha=0.5)\n",
    "\n",
    "# Title and labels\n",
    "ax.set_title('Volcano Plot')\n",
    "ax.set_xlabel('Log2 Fold Change')\n",
    "ax.set_ylabel('-Log10(Adjusted P-value)')\n",
    "\n",
    "# Draw horizontal line to show the significance threshold\n",
    "ax.axhline(y=-np.log10(sig_threshold), linestyle='--', color='grey')\n",
    "\n",
    "# Show the plot\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7fa1e84f-9787-4cd1-9c27-f6cafea66843",
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "source": [
    "## comparing with CV genes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6e6a7d18-43c1-47d7-bb95-ae7002e3c68a",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# Step 1: Get the same number of rows from df with the highest 'cv' values\n",
    "top_cv_genes_df = df.nlargest(significant_genes.shape[0], 'cv')\n",
    "\n",
    "# Step 2: Compare the genes (indexes) from the two DataFrames and find the common ones\n",
    "common_genes = set(significant_genes.index).intersection(set(top_cv_genes_df.index))\n",
    "\n",
    "# Print the common genes\n",
    "print(f\"Top CV genes: {len(top_cv_genes_df)}\")\n",
    "\n",
    "print(f\"Differential expression analysis genes: {significant_genes.shape[0]}\")\n",
    "\n",
    "print(f\"Common genes between significant genes and top CV genes: {len(common_genes)}\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0d84088b-3c76-4f6d-9bac-1c7378d58af9",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "# Create a figure and a set of subplots.\n",
    "fig, ax = plt.subplots(figsize=(10, 6))\n",
    "\n",
    "# Initialize an empty list to store the colors\n",
    "colors = []\n",
    "\n",
    "# Iterate over the rows of the DataFrame and assign colors\n",
    "for idx, row in results_df.iterrows():\n",
    "    if idx in common_genes:\n",
    "        colors.append('green')  # Green for genes in common_genes\n",
    "    elif row['adjusted_p_value'] < sig_threshold:\n",
    "        colors.append('r')  # Red for significant genes\n",
    "    else:\n",
    "        colors.append('b')  # Blue for other genes\n",
    "\n",
    "# Scatter plot\n",
    "ax.scatter(results_df['log2_fold_change'], results_df['neg_log10_p_value'], c=colors, alpha=0.5)\n",
    "\n",
    "# Title and labels\n",
    "ax.set_title('Volcano Plot')\n",
    "ax.set_xlabel('Log2 Fold Change')\n",
    "ax.set_ylabel('-Log10(Adjusted P-value)')\n",
    "\n",
    "# Draw horizontal line to show the significance threshold\n",
    "ax.axhline(y=-np.log10(sig_threshold), linestyle='--', color='grey')\n",
    "\n",
    "# Show the plot\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1a112f01-da02-435c-af3e-d4b3caacc0b1",
   "metadata": {},
   "source": [
    "# "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ffa28a86-f5b7-4635-9255-384afe239b82",
   "metadata": {},
   "source": [
    "# Clustering"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "68de7c0b-3cd1-4c54-8121-0644690e3d8f",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.preprocessing import StandardScaler\n",
    "from sklearn.decomposition import PCA\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "# Step 1: Extract rows corresponding to common_genes from df\n",
    "common_genes_df = df.loc[list(common_genes)]\n",
    "\n",
    "# Step 2: Standardize the data\n",
    "scaler = StandardScaler()\n",
    "standardized_data = scaler.fit_transform(common_genes_df.T)  # Transposing as we want to standardize the features (genes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dcd59254-ab34-4ea5-a987-d1439a9de5c6",
   "metadata": {},
   "outputs": [],
   "source": [
    "common_genes_df.columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7689bf25-1b23-4feb-a633-b64ff6de056e",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "from sklearn.cluster import HDBSCAN\n",
    "\n",
    "# Create a DBSCAN object with specific epsilon and minimum samples\n",
    "hdbscan = HDBSCAN(min_samples=3)\n",
    "\n",
    "# Fit the model to the scaled data\n",
    "clustering_labels = hdbscan.fit_predict(standardized_data.T)\n",
    "\n",
    "# Attach the DBSCAN labels to the filtered DataFrame\n",
    "common_genes_df.insert(0, \"clustering_labels\", pd.Series(clustering_labels, index=common_genes_df.index))\n",
    "common_genes_df.columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a9d5f346-c938-4218-acb1-87b7dbcdf831",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.unique(clustering_labels)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0734c850-14e2-477d-a263-569997e08cb6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Perform PCA on the standardized data\n",
    "pca = PCA(n_components=0.9, svd_solver='full')\n",
    "pca_model = pca.fit(standardized_data.T)\n",
    "principalComponents = pca_model.transform(standardized_data.T)\n",
    "\n",
    "pc_df = pd.DataFrame(principalComponents, columns=[f'PC{i}' for i in range(1, pca.n_components_ + 1)])\n",
    "pc_df.insert(0, \"clustering_labels\", pd.Series(clustering_labels))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f822f0ff-b1b1-4057-89fe-0dc7c8205acc",
   "metadata": {},
   "outputs": [],
   "source": [
    "pca.n_components_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b44eab07-0703-4557-9939-1bef8cf76f12",
   "metadata": {},
   "outputs": [],
   "source": [
    "pc_df['clustering_labels'].unique()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f136a3ab-d54b-4506-8636-76ca1a30bb25",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# Set the color palette as needed\n",
    "palette = sns.color_palette(\"colorblind\", len(set(clustering_labels)))\n",
    "\n",
    "# Create the pairplot colored by clustering_labels\n",
    "sns.pairplot(pc_df, hue='clustering_labels', palette=palette, diag_kind='kde', corner=True)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b504dbda-c599-40db-81c5-d667ee0bfed8",
   "metadata": {},
   "outputs": [],
   "source": [
    "common_genes_df[common_genes_df[\"clustering_labels\"]== 4]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d4ab1633-b36f-4115-8f34-56808a33162f",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.unique(clustering_labels)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4cb6b0c0-8841-4fdc-8ddc-cb63277b4b0d",
   "metadata": {},
   "source": [
    "# Graphs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cf8e7316-aedf-4a67-bec7-8fa3f0953af5",
   "metadata": {},
   "outputs": [],
   "source": [
    "expression_data = common_genes_df.drop(columns=['mean', 'std', 'cv'])\n",
    "correlation_matrix = expression_data.T.corr()\n",
    "\n",
    "# Make all values in correlation_matrix positive\n",
    "correlation_matrix = np.abs(correlation_matrix)\n",
    "\n",
    "# Set values below threshold to 0\n",
    "threshold = 0.5\n",
    "correlation_matrix.where(correlation_matrix >= threshold, 0, inplace=True)\n",
    "sns.heatmap(correlation_matrix)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "40b58c4f-8976-40c6-b588-3e972eb1d5d6",
   "metadata": {},
   "outputs": [],
   "source": [
    "from graph_tool.all import *\n",
    "import numpy as np\n",
    "\n",
    "# Initialize a graph\n",
    "g = Graph(directed=False)\n",
    "\n",
    "# Create a vertex property for gene names\n",
    "v_label = g.new_vertex_property(\"string\")\n",
    "\n",
    "# Create a dictionary to map gene names to vertex index\n",
    "gene_to_vertex = {}\n",
    "genes_to_index = {}\n",
    "\n",
    "# Create vertices for each gene\n",
    "for i, gene in enumerate(expression_data.index):\n",
    "    v = g.add_vertex()\n",
    "    v_label[v] = gene\n",
    "    gene_to_vertex[gene] = v\n",
    "    genes_to_index[gene]=i\n",
    "\n",
    "# Create an edge property for edge weights (correlation)\n",
    "edge_weight = g.new_edge_property(\"double\")\n",
    "\n",
    "# Iterate over the upper triangle of the correlation matrix and add edges\n",
    "for i, gene1 in enumerate(expression_data.index):\n",
    "    for j, gene2 in enumerate(expression_data.index):\n",
    "        if i >= j:  # Avoid adding duplicate edges\n",
    "            continue\n",
    "        correlation = correlation_matrix.loc[gene1, gene2]\n",
    "        \n",
    "        # Optional: You might want to set a threshold to add an edge, e.g. only add edges for correlation > 0.5\n",
    "        if correlation > 0.6:\n",
    "            v1 = gene_to_vertex[gene1]\n",
    "            v2 = gene_to_vertex[gene2]\n",
    "            e = g.add_edge(v1, v2)\n",
    "            edge_weight[e] = np.abs(correlation) * 3  # Scale the correlation for visualization purposes\n",
    "\n",
    "# Save vertex and edge properties to the graph\n",
    "g.vertex_properties[\"label\"] = v_label\n",
    "g.edge_properties[\"weight\"] = edge_weight"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8b0981f3-6545-41b9-a3c8-a143c30f715a",
   "metadata": {},
   "outputs": [],
   "source": [
    "from graph_tool.all import minimize_nested_blockmodel_dl, minimize_blockmodel_dl\n",
    "\n",
    "state = minimize_nested_blockmodel_dl(g)\n",
    "\n",
    "state.draw(output=\"power_nested_mdl.svg\", vertex_text=g.vp[\"label\"], vertex_text_position='centered')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "675393fe-6b2e-4b2e-b23d-5d423ef8fe4e",
   "metadata": {},
   "source": [
    "# Comparaison with litterature"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d144ff0d-d932-4a2a-832f-46809e0ae892",
   "metadata": {},
   "outputs": [],
   "source": [
    "state"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "266c2c3a-36ea-4443-bc7c-c927613c806e",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# given array\n",
    "clusters_array = state.levels[0].get_blocks().get_array() \n",
    "\n",
    "# create an empty dictionary to store the mapping of gene names to clusters\n",
    "gene_to_cluster = {}\n",
    "\n",
    "# map each gene name to its corresponding cluster\n",
    "for gene, index in genes_to_index.items():\n",
    "    if index < len(clusters_array):  # to ensure the index exists in clusters_array\n",
    "        gene_to_cluster[gene] = clusters_array[index]\n",
    "\n",
    "gene_to_cluster_df = pd.DataFrame(list(gene_to_cluster.items()), columns=['Gene', 'Cluster'])\n",
    "\n",
    "# available clusters\n",
    "gene_to_cluster_df['Cluster'].unique()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "29aa1d8f-312f-46f9-8ac4-cd5bccc12c79",
   "metadata": {},
   "outputs": [],
   "source": [
    "selected_cluster = gene_to_cluster_df[gene_to_cluster_df['Cluster'] == 873]\n",
    "list(selected_cluster['Gene'])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a9ea0bc1-b917-4899-b1c0-12934f120424",
   "metadata": {},
   "source": [
    "## Compare with literrature"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "788e24a5-4b03-4a38-9b8f-75ec5efdf0b0",
   "metadata": {},
   "source": [
    "### Prepare request"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b50f619b-b0c7-4473-b19a-9d64e8d220fd",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "import requests ## python -m pip install requests\n",
    "\n",
    "string_api_url = \"https://version-11-5.string-db.org/api\"\n",
    "output_format = \"tsv-no-header\"\n",
    "method = \"get_string_ids\"\n",
    "\n",
    "## Set parameters\n",
    "params = {\n",
    "    \"identifiers\" : \"\\r\".join(list(selected_cluster['Gene'])), # your protein list\n",
    "    \"species\" : 10090, # species NCBI identifier \n",
    "    \"limit\" : 1, # only one (best) identifier per input protein\n",
    "    \"echo_query\" : 1, # see your input identifiers in the output\n",
    "    \"caller_identity\" : \"www.awesome_app.org\" # your app name\n",
    "}\n",
    "\n",
    "## Construct URL\n",
    "request_url = \"/\".join([string_api_url, output_format, method])\n",
    "\n",
    "## Call STRING\n",
    "results = requests.post(request_url, data=params)\n",
    "\n",
    "identifiers_list = []\n",
    "\n",
    "# Read and parse the results\n",
    "for line in results.text.strip().split(\"\\n\"):\n",
    "    l = line.split(\"\\t\")\n",
    "    print(l)\n",
    "    actual_identifier = l[6]\n",
    "    \n",
    "    # Append string_identifier to the list\n",
    "    identifiers_list.append(actual_identifier)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6f164093-f3dd-46b3-8b42-6d3e42f6023b",
   "metadata": {},
   "outputs": [],
   "source": [
    "interaction = set(identifiers_list) & set(list(selected_cluster['Gene']))\n",
    "print(len(interaction) )\n",
    "print(len(identifiers_list) )\n",
    "print(len(selected_cluster['Gene'])) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4a0fdc2c-15b3-4688-a466-b07ce1696a46",
   "metadata": {},
   "outputs": [],
   "source": [
    "import requests ## python -m pip install requests\n",
    "\n",
    "string_api_url = \"https://version-11-5.string-db.org/api\"\n",
    "output_format = \"tsv-no-header\"\n",
    "method = \"interaction_partners\"\n",
    "\n",
    "\n",
    "my_genes = string_identifiers_list\n",
    "\n",
    "## Construct the request\n",
    "request_url = \"/\".join([string_api_url, output_format, method])\n",
    "\n",
    "\n",
    "## Set parameters\n",
    "params = {\n",
    "    \"identifiers\" : \"%0d\".join(my_genes), # your protein\n",
    "    \"species\" : 10090, # species NCBI identifier \n",
    "    \"limit\" : 10,\n",
    "    \"caller_identity\" : \"www.awesome_app.org\" # your app name\n",
    "}\n",
    "\n",
    "# Request\n",
    "response = requests.post(request_url, data=params)\n",
    "\n",
    "\n",
    "\n",
    "## Read and parse the results\n",
    "# Iterate over each line and append a dictionary to the list\n",
    "data_list = []\n",
    "for line in response.text.strip().split(\"\\n\"):\n",
    "    l = line.strip().split(\"\\t\")\n",
    "    row_dict = {\n",
    "        'query_ensp': l[0],\n",
    "        'query_name': l[2],\n",
    "        'partner_ensp': l[1],\n",
    "        'partner_name': l[3],\n",
    "        'combined_score': l[5]\n",
    "    }\n",
    "    data_list.append(row_dict)\n",
    "\n",
    "# Convert the list of dictionaries to a DataFrame\n",
    "df = pd.DataFrame(data_list)\n",
    "\n",
    "# Print the DataFrame\n",
    "df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "07f6c9b6-3a7b-46e0-afa1-e70a4ba0cc70",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "interaction = set(list(selected_cluster['Gene'])) & set(list(df[\"partner_name\"]))\n",
    "print(interaction)  \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "36c58659-312e-4304-9b4c-469d3bdd8f9f",
   "metadata": {},
   "outputs": [],
   "source": [
    "selected_cluster['Gene']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "399a3948-92bb-45e2-875c-16171ab8028a",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "transcriptome_env",
   "language": "python",
   "name": "transcriptome_env"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}