vcf-processing-and-clustering.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 179,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pysam import VariantFile\n",
    "from sklearn.cluster import KMeans\n",
    "\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import math"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 196,
   "metadata": {},
   "outputs": [],
   "source": [
    "vcf_file = 'DRR131565.variants.HC_init.wAnnot.vcf'\n",
    "mca_file = 'MVP_MCA_results.csv'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Parse VCF File and get variants"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 197,
   "metadata": {},
   "outputs": [],
   "source": [
    "def getGenesVariantList(vcf_file):\n",
    "    \n",
    "    vcf = VariantFile(vcf_file)\n",
    "    \n",
    "    variants = set()\n",
    "        \n",
    "    for rec in vcf.fetch():\n",
    "         for gene_name in rec.info[\"GENE\"]:\n",
    "                variants.add(gene_name)\n",
    "\n",
    "    return variants"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Apply clustering to MCA to identify genes closest to the closest experiment"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 189,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Uses a method to dynamically determine number of clusters, \n",
    "# and then runs kmeans\n",
    "def getKMeanLabels(X):\n",
    "    \n",
    "    ##################\n",
    "    # Code borrowed from Jessica Temporal \n",
    "    # https://jtemporal.com/kmeans-and-elbow-method/\n",
    "    def calculate_wcss(data):\n",
    "            wcss = []\n",
    "            for n in range(2, 21):\n",
    "                kmeans = KMeans(n_clusters=n)\n",
    "                kmeans.fit(X=data)\n",
    "                wcss.append(kmeans.inertia_)\n",
    "\n",
    "            return wcss\n",
    "\n",
    "    def optimal_number_of_clusters(wcss):\n",
    "        x1, y1 = 2, wcss[0]\n",
    "        x2, y2 = 20, wcss[len(wcss)-1]\n",
    "\n",
    "        distances = []\n",
    "        for i in range(len(wcss)):\n",
    "            x0 = i+2\n",
    "            y0 = wcss[i]\n",
    "            numerator = abs((y2-y1)*x0 - (x2-x1)*y0 + x2*y1 - y2*x1)\n",
    "            denominator = math.sqrt((y2 - y1)**2 + (x2 - x1)**2)\n",
    "            distances.append(numerator/denominator)\n",
    "\n",
    "        return distances.index(max(distances)) + 2\n",
    "    \n",
    "    # calculating the within clusters sum-of-squares for 19 cluster amounts\n",
    "    sum_of_squares = calculate_wcss(X)\n",
    "\n",
    "    # calculating the optimal number of clusters\n",
    "    n = optimal_number_of_clusters(sum_of_squares)\n",
    "    ###################\n",
    "    \n",
    "    kmeans = KMeans(n_clusters=n)\n",
    "    kmeans.fit(X)\n",
    "    y_kmeans = kmeans.predict(X)\n",
    "    \n",
    "    return y_kmeans, kmeans.cluster_centers_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 206,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get 2D distance between A and B\n",
    "def getDist(A,B):\n",
    "    return math.sqrt((A[0]-B[0])**2 + (A[1]-B[1])**2)\n",
    "\n",
    "# gets the minimum distance from point to \n",
    "# a center in centers\n",
    "def getMinDist(point, centers):\n",
    "    \n",
    "    min_dist = getDist(point,centers[0])\n",
    "    min_idx = 0\n",
    "    \n",
    "    for idx, center in enumerate(centers):\n",
    "        dist = getDist(point, center)\n",
    "        if dist < min_dist: \n",
    "            min_dist = dist\n",
    "            min_idx = idx\n",
    "            \n",
    "    return (min_idx, min_dist)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 218,
   "metadata": {},
   "outputs": [],
   "source": [
    "def getGenesOfInterestFromMCA(mca_file):\n",
    "    \n",
    "    mca_df = pd.read_csv('MVP_MCA_results.csv')\n",
    "\n",
    "    gene_df = mca_df[mca_df['AssociationType'].isin(['Gene' , 'DrugTerm'])]\n",
    "    no_gene_df = mca_df[~mca_df['AssociationType'].isin(['Gene', 'DrugTerm'])]\n",
    "    \n",
    "    X = np.array(list([list(a) for a in zip(gene_df['x'], gene_df['y'])]))\n",
    "    \n",
    "    # cluster genes\n",
    "    y_kmeans, cluster_centers = getKMeanLabels(X)\n",
    "    gene_df['cluster_center'] = y_kmeans\n",
    "    \n",
    "    # Display clustering\n",
    "    # plt.scatter(X[:, 0], X[:, 1], c=y_kmeans, s=50, cmap='viridis')\n",
    "    # plt.show()\n",
    "    \n",
    "    # add the distance from the closest center to each \n",
    "    no_gene_df['closestClusterCenter'] = \\\n",
    "        no_gene_df.apply(lambda row: getMinDist([row.x, row.y], cluster_centers)[0], axis=1)\n",
    "    \n",
    "    no_gene_df['distToClusterCenter'] = \\\n",
    "        no_gene_df.apply(lambda row: getMinDist([row.x, row.y], cluster_centers)[1], axis=1)\n",
    "    \n",
    "    no_gene_df = no_gene_df.sort_values(by=['distToClusterCenter'])\n",
    "    \n",
    "    cluster_of_interest = no_gene_df['closestClusterCenter'].iloc[0]\n",
    "    \n",
    "    return list(gene_df[gene_df['cluster_center']==cluster_of_interest]['Value'].unique())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 220,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "<ipython-input-218-45beca6bd270>:12: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame.\n",
      "Try using .loc[row_indexer,col_indexer] = value instead\n",
      "\n",
      "See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy\n",
      "  gene_df['cluster_center'] = y_kmeans\n",
      "<ipython-input-218-45beca6bd270>:19: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame.\n",
      "Try using .loc[row_indexer,col_indexer] = value instead\n",
      "\n",
      "See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy\n",
      "  no_gene_df['closestClusterCenter'] = \\\n",
      "<ipython-input-218-45beca6bd270>:22: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame.\n",
      "Try using .loc[row_indexer,col_indexer] = value instead\n",
      "\n",
      "See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy\n",
      "  no_gene_df['distToClusterCenter'] = \\\n"
     ]
    }
   ],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 224,
   "metadata": {},
   "outputs": [],
   "source": [
    "def getGenesOfInterest(gene_list, variants):\n",
    "    \n",
    "    # currently just taking the intersection, but\n",
    "    # can be made more sophisticated\n",
    "    return set(gene_list).intersection(variants)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 225,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "<ipython-input-218-45beca6bd270>:12: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame.\n",
      "Try using .loc[row_indexer,col_indexer] = value instead\n",
      "\n",
      "See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy\n",
      "  gene_df['cluster_center'] = y_kmeans\n",
      "<ipython-input-218-45beca6bd270>:19: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame.\n",
      "Try using .loc[row_indexer,col_indexer] = value instead\n",
      "\n",
      "See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy\n",
      "  no_gene_df['closestClusterCenter'] = \\\n",
      "<ipython-input-218-45beca6bd270>:22: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame.\n",
      "Try using .loc[row_indexer,col_indexer] = value instead\n",
      "\n",
      "See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy\n",
      "  no_gene_df['distToClusterCenter'] = \\\n"
     ]
    }
   ],
   "source": [
    "variants = getGenesVariantList(filename)\n",
    "gene_list = getGenesOfInterestFromMCA(mca_file)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 226,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'ATL2',\n",
       " 'EFCAB11',\n",
       " 'FBXO34',\n",
       " 'KIF22',\n",
       " 'KLHL29',\n",
       " 'MRPS5',\n",
       " 'PAQR5',\n",
       " 'RGS20',\n",
       " 'SACM1L',\n",
       " 'ZFAND5'}"
      ]
     },
     "execution_count": 226,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "getGenesOfInterest(gene_list, variants)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}