FigS1and2.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "de0a9802",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import libraries. \n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import seaborn as sns\n",
    "import matplotlib.pyplot as plt\n",
    "from kneed import KneeLocator\n",
    "\n",
    "# Dimension reduction tools\n",
    "from sklearn.decomposition import PCA as PCA\n",
    "from sklearn.manifold import TSNE\n",
    "import umap "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b71509be",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Creating Color Palette for Different Populations\n",
    "\n",
    "#for Oryx datasets:\n",
    "custom_palette = {\n",
    "    \"EAD_A\": \"#cada45\", #green\n",
    "    \"EAD_B\": \"#d4a2e1\", #purple\n",
    "    \"EEP\": \"#55e0c6\", #blue\n",
    "    \"USA\": \"#f0b13c\", #orange\n",
    "}\n",
    "\n",
    "#for Galapagos Datasets:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f8f2e145",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#Loading the population data and covariance matrix\n",
    "\n",
    "#load population data\n",
    "population_names = pd.read_csv('input_files/oryx_pop_info_sorted_46_final.txt', sep='\\t', header=0)\n",
    "#load the covariance matrix\n",
    "filename='input_files/oryx_6xyh_1K.cov'\n",
    "cov_mat= pd.read_csv(filename, sep=' ', header=None)\n",
    "#Generating the pandas dataframe called Data_Struct\n",
    "Data_Struct=population_names\n",
    "filenameforplot='SO_6x'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0f1a7937",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Function to plot the scree plot\n",
    "def plot_scree(explained_variance,filename_title,elbow_point):\n",
    "    plt.figure(figsize=(8, 5))\n",
    "    # Convert to a simple list if it's not already\n",
    "    explained_variance = list(explained_variance)\n",
    "    plt.plot(range(1, len(explained_variance) + 1), explained_variance, marker='o', linestyle='--')\n",
    "    plt.axvline(x=elbow_point, color='r', linestyle='--')\n",
    "    plt.text(elbow_point + 0.1, max(explained_variance) * 0.9, f'Elbow point: {elbow_point}', color='red', verticalalignment='center', fontsize=18)\n",
    "    plt.title(f'{filename_title}', fontsize=18)\n",
    "    plt.xlabel('Number of Components', fontsize=18)\n",
    "    plt.ylabel('Variance Explained', fontsize=18)\n",
    "    plt.grid()\n",
    "    #plt.show()\n",
    "     # Adjust tick label size\n",
    "    plt.xticks(fontsize=18)\n",
    "    plt.yticks(fontsize=18)\n",
    "     # Save the plot as a PDF file\n",
    "    pdf_filename = f'{filename_title}.pdf'  # Define the PDF file name\n",
    "    plt.savefig(pdf_filename, format='pdf')  # Save the figure as a PDF\n",
    "    plt.close()  # Close the figure to avoid displaying it in the notebook\n",
    "\n",
    "# Function to find the elbow point\n",
    "def find_elbow_point(explained_variance, sensitivity=1.0):\n",
    "    explained_variance = list(explained_variance)\n",
    "    kneedle = KneeLocator(range(1, len(explained_variance) + 1), explained_variance, \n",
    "                          curve='convex', direction='decreasing', \n",
    "                          S=sensitivity, interp_method='polynomial')\n",
    "    return kneedle.elbow\n",
    "\n",
    "\n",
    "#Calculate PCA\n",
    "#convert covariance matrix to numpy array\n",
    "cov_mat_np=cov_mat.to_numpy()\n",
    "\n",
    "# calculate eigen vectors and eigen values from the initial covariance matrix\n",
    "eigen_vals, eigen_vecs = np.linalg.eig(cov_mat_np)\n",
    "eigen_pairs = [(np.abs(eigen_vals[i]), eigen_vecs[:, i]) for i in range(len(eigen_vals))]\n",
    "eigen_pairs.sort(key=lambda k: k[0], reverse=True)\n",
    "feature_vector = np.hstack([eigen_pairs[i][1][:, np.newaxis] for i in range(len(eigen_vals))])\n",
    "principal_components = cov_mat_np.dot(feature_vector) \n",
    "\n",
    "# sorting them from largest to smallest\n",
    "idx = eigen_vals.argsort()[::-1]   \n",
    "eigenValues = eigen_vals[idx]\n",
    "eigenVectors = eigen_vecs[:,idx]\n",
    "\n",
    "# calculating the total explained variance\n",
    "expl_pre=eigenValues/sum(eigenValues)\n",
    "expl=np.cumsum(expl_pre)\n",
    "\n",
    "expl_df=pd.DataFrame(expl_pre*100,columns=['explained_variance'])\n",
    "expl_df['cumulative_expl']=expl*100\n",
    "expl_df.set_index(np.arange(1, eigenVectors.shape[0] + 1), inplace=True)\n",
    "\n",
    "\n",
    "# Plot the scree plot\n",
    "#plot_filename = f'scree_plot_{filename_title}.png'\n",
    "\n",
    " # Find the elbow point\n",
    "elbow_point = find_elbow_point(expl_df['explained_variance'])\n",
    "print(\"Optimal number of principal components):\", elbow_point)\n",
    "\n",
    "plot_scree(expl_df['explained_variance'],filenameforplot,elbow_point)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bbf4cb09",
   "metadata": {},
   "outputs": [],
   "source": [
    "#to generate FigS1\n",
    "from pdf2image import convert_from_path\n",
    "\n",
    "images = convert_from_path('panelA.pdf')\n",
    "for i, image in enumerate(images):\n",
    "    image.save(f'panelA_page{i}.png', 'PNG')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "beyondpca",
   "language": "python",
   "name": "python3"
  },
  "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.10.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}