Real_Data/Cusanovich_2018/run_clustering_cusanovich2018.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import scanpy as sc\n",
    "import os\n",
    "from sklearn.cluster import KMeans\n",
    "from sklearn.cluster import AgglomerativeClustering\n",
    "from sklearn.metrics.cluster import adjusted_rand_score\n",
    "from sklearn.metrics.cluster import adjusted_mutual_info_score\n",
    "from sklearn.metrics.cluster import homogeneity_score\n",
    "import rpy2.robjects as robjects\n",
    "from rpy2.robjects import pandas2ri"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_metrics = pd.DataFrame(columns=['ARI_Louvain','ARI_kmeans','ARI_HC',\n",
    "                                   'AMI_Louvain','AMI_kmeans','AMI_HC',\n",
    "                                   'Homogeneity_Louvain','Homogeneity_kmeans','Homogeneity_HC'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "workdir = './output/'\n",
    "path_fm = os.path.join(workdir,'feature_matrices/')\n",
    "path_clusters = os.path.join(workdir,'clusters/')\n",
    "path_metrics = os.path.join(workdir,'metrics/')\n",
    "os.system('mkdir -p '+path_clusters)\n",
    "os.system('mkdir -p '+path_metrics)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "13\n"
     ]
    }
   ],
   "source": [
    "metadata = pd.read_csv('./input/metadata.tsv',sep='\\t',index_col=0)\n",
    "num_clusters = len(np.unique(metadata['label']))\n",
    "print(num_clusters)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "files = [x for x in os.listdir(path_fm) if x.startswith('FM')]\n",
    "len(files)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['FM_SnapATAC_cusanovich2018.rds']"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def getNClusters(adata,n_cluster,range_min=0,range_max=3,max_steps=20):\n",
    "    this_step = 0\n",
    "    this_min = float(range_min)\n",
    "    this_max = float(range_max)\n",
    "    while this_step < max_steps:\n",
    "        print('step ' + str(this_step))\n",
    "        this_resolution = this_min + ((this_max-this_min)/2)\n",
    "        sc.tl.louvain(adata,resolution=this_resolution)\n",
    "        this_clusters = adata.obs['louvain'].nunique()\n",
    "        \n",
    "        print('got ' + str(this_clusters) + ' at resolution ' + str(this_resolution))\n",
    "        \n",
    "        if this_clusters > n_cluster:\n",
    "            this_max = this_resolution\n",
    "        elif this_clusters < n_cluster:\n",
    "            this_min = this_resolution\n",
    "        else:\n",
    "            return(this_resolution, adata)\n",
    "        this_step += 1\n",
    "    \n",
    "    print('Cannot find the number of clusters')\n",
    "    print('Clustering solution from last iteration is used:' + str(this_clusters) + ' at resolution ' + str(this_resolution))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "SnapATAC\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/data/pinello/SHARED_SOFTWARE/anaconda3/envs/ATACseq_clustering/lib/python3.7/site-packages/rpy2/robjects/pandas2ri.py:191: FutureWarning: from_items is deprecated. Please use DataFrame.from_dict(dict(items), ...) instead. DataFrame.from_dict(OrderedDict(items)) may be used to preserve the key order.\n",
      "  res = PandasDataFrame.from_items(items)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "step 0\n",
      "got 30 at resolution 1.5\n",
      "step 1\n",
      "got 21 at resolution 0.75\n",
      "step 2\n",
      "got 19 at resolution 0.375\n",
      "step 3\n",
      "got 17 at resolution 0.1875\n",
      "step 4\n",
      "got 12 at resolution 0.09375\n",
      "step 5\n",
      "got 15 at resolution 0.140625\n",
      "step 6\n",
      "got 13 at resolution 0.1171875\n"
     ]
    }
   ],
   "source": [
    "for file in files:\n",
    "    file_split = file.split('_')\n",
    "    method = file_split[1]\n",
    "    dataset = file_split[2].split('.')[0]\n",
    "    if(len(file_split)>3):\n",
    "        method = method + '_' + '_'.join(file_split[3:]).split('.')[0]\n",
    "    print(method)\n",
    "\n",
    "    pandas2ri.activate()\n",
    "    readRDS = robjects.r['readRDS']\n",
    "    df_rds = readRDS(os.path.join(path_fm,file))\n",
    "    fm_mat = pandas2ri.ri2py(robjects.r['data.frame'](robjects.r['as.matrix'](df_rds)))\n",
    "    fm_mat.fillna(0,inplace=True)\n",
    "    fm_mat.columns = metadata.index\n",
    "    \n",
    "    adata = sc.AnnData(fm_mat.T)\n",
    "    adata.var_names_make_unique()\n",
    "    adata.obs = metadata.loc[adata.obs.index,]\n",
    "    df_metrics.loc[method,] = \"\"\n",
    "    #Louvain\n",
    "    sc.pp.neighbors(adata, n_neighbors=15,use_rep='X')\n",
    "#     sc.tl.louvain(adata)\n",
    "    getNClusters(adata,n_cluster=num_clusters)\n",
    "    #kmeans\n",
    "    kmeans = KMeans(n_clusters=num_clusters, random_state=2019).fit(adata.X)\n",
    "    adata.obs['kmeans'] = pd.Series(kmeans.labels_,index=adata.obs.index).astype('category')\n",
    "    #hierachical clustering\n",
    "    hc = AgglomerativeClustering(n_clusters=num_clusters).fit(adata.X)\n",
    "    adata.obs['hc'] = pd.Series(hc.labels_,index=adata.obs.index).astype('category')\n",
    "    #clustering metrics\n",
    "    \n",
    "    #adjusted rank index\n",
    "    ari_louvain = adjusted_rand_score(adata.obs['label'], adata.obs['louvain'])\n",
    "    ari_kmeans = adjusted_rand_score(adata.obs['label'], adata.obs['kmeans'])\n",
    "    ari_hc = adjusted_rand_score(adata.obs['label'], adata.obs['hc'])\n",
    "    #adjusted mutual information\n",
    "    ami_louvain = adjusted_mutual_info_score(adata.obs['label'], adata.obs['louvain'],average_method='arithmetic')\n",
    "    ami_kmeans = adjusted_mutual_info_score(adata.obs['label'], adata.obs['kmeans'],average_method='arithmetic')   \n",
    "    ami_hc = adjusted_mutual_info_score(adata.obs['label'], adata.obs['hc'],average_method='arithmetic')\n",
    "    #homogeneity\n",
    "    homo_louvain = homogeneity_score(adata.obs['label'], adata.obs['louvain'])\n",
    "    homo_kmeans = homogeneity_score(adata.obs['label'], adata.obs['kmeans'])\n",
    "    homo_hc = homogeneity_score(adata.obs['label'], adata.obs['hc'])\n",
    "\n",
    "    df_metrics.loc[method,['ARI_Louvain','ARI_kmeans','ARI_HC']] = [ari_louvain,ari_kmeans,ari_hc]\n",
    "    df_metrics.loc[method,['AMI_Louvain','AMI_kmeans','AMI_HC']] = [ami_louvain,ami_kmeans,ami_hc]\n",
    "    df_metrics.loc[method,['Homogeneity_Louvain','Homogeneity_kmeans','Homogeneity_HC']] = [homo_louvain,homo_kmeans,homo_hc] \n",
    "    adata.obs[['louvain','kmeans','hc']].to_csv(os.path.join(path_clusters ,method + '_clusters.tsv'),sep='\\t')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_metrics.to_csv(path_metrics+'clustering_scores.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>ARI_Louvain</th>\n",
       "      <th>ARI_kmeans</th>\n",
       "      <th>ARI_HC</th>\n",
       "      <th>AMI_Louvain</th>\n",
       "      <th>AMI_kmeans</th>\n",
       "      <th>AMI_HC</th>\n",
       "      <th>Homogeneity_Louvain</th>\n",
       "      <th>Homogeneity_kmeans</th>\n",
       "      <th>Homogeneity_HC</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>SnapATAC</th>\n",
       "      <td>0.385046</td>\n",
       "      <td>0.342715</td>\n",
       "      <td>0.238682</td>\n",
       "      <td>0.619034</td>\n",
       "      <td>0.553361</td>\n",
       "      <td>0.554868</td>\n",
       "      <td>0.595526</td>\n",
       "      <td>0.539072</td>\n",
       "      <td>0.515852</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "         ARI_Louvain ARI_kmeans    ARI_HC AMI_Louvain AMI_kmeans    AMI_HC  \\\n",
       "SnapATAC    0.385046   0.342715  0.238682    0.619034   0.553361  0.554868   \n",
       "\n",
       "         Homogeneity_Louvain Homogeneity_kmeans Homogeneity_HC  \n",
       "SnapATAC            0.595526           0.539072       0.515852  "
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df_metrics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:ATACseq_clustering]",
   "language": "python",
   "name": "conda-env-ATACseq_clustering-py"
  },
  "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}