Method iClusterPlus.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "Copia di Method iClusterPlus.ipynb",
      "provenance": [],
      "collapsed_sections": [],
      "toc_visible": true
    },
    "kernelspec": {
      "display_name": "Python 3",
      "name": "python3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "cells": [
    {
      "cell_type": "code",
      "metadata": {
        "id": "a4n13KiPNOpw"
      },
      "source": [
        "import os\n",
        "import gdown\n",
        "import pandas as pd\n",
        "import matplotlib.pyplot as plt\n",
        "import seaborn as sns\n",
        "import numpy as np\n",
        "import scipy as sp\n",
        "import matplotlib.cm as cm\n",
        "\n",
        "#np.random.seed(0)\n",
        "\n",
        "from collections import Counter\n",
        "from sklearn.preprocessing import StandardScaler, MinMaxScaler\n",
        "from sklearn.model_selection import StratifiedShuffleSplit\n",
        "from sklearn.decomposition import PCA\n",
        "from sklearn.cluster import KMeans, DBSCAN, SpectralClustering, AgglomerativeClustering\n",
        "from sklearn.metrics import silhouette_score, confusion_matrix, adjusted_mutual_info_score, adjusted_rand_score, silhouette_samples\n",
        "from sklearn.neighbors import NearestNeighbors\n",
        "from scipy.spatial import distance_matrix\n",
        "\n",
        "import warnings\n",
        "warnings.simplefilter(action='ignore', category=FutureWarning)\n",
        "\n",
        "# Import utils.py\n",
        "!gdown 'https://drive.google.com/uc?id=13I5w4WajPg6MObtLPQjxznm8w5hKlEY0' -O ./utils.py\n",
        "from utils import *"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "B0md7FKGMY-S"
      },
      "source": [
        "# Retrieve datasets and reduce their dimension with PCA"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "tAX5P3cTMhMf"
      },
      "source": [
        "The dimensionality reduction is needed, because iClusterPlus is very computationally heavy and Colab RAM gets saturated by running the algorithm on the original omics (Lung and Kidney in particular)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "76gyEjGjNJxg"
      },
      "source": [
        "## Synthetic"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "LYrY3l7aNWTj"
      },
      "source": [
        "Load omics"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "ph_CtKVyMsji"
      },
      "source": [
        "# Load dataset\n",
        "!wget --no-check-certificate 'https://docs.google.com/uc?export=download&id=12PArkc1RsOm2437mbysxRF4hQMddZOsc' -O ./mRNA.txt\n",
        "!wget --no-check-certificate 'https://docs.google.com/uc?export=download&id=1aJkDF0ckxzY4vsnS53s-V89DdAnRVbPo' -O ./meth.txt\n",
        "!wget --no-check-certificate 'https://docs.google.com/uc?export=download&id=1iS4u1SZH6r_Dvs7qRSKC444kc_bqmGhJ' -O ./prot.txt\n",
        "!wget --no-check-certificate 'https://docs.google.com/uc?export=download&id=1UtHj4BzBx5hnQkkklJ9ugU5ERjLhjx4W' -O ./clusters.txt\n",
        "\n",
        "ds = {}     # this will contain each omic\n",
        "omics = ['mRNA','meth','prot']\n",
        "for omic_name in omics:\n",
        "    path = omic_name + \".txt\"\n",
        "    ds[omic_name] = pd.read_csv(path, sep='\\t', index_col=0)\n",
        "    ds[omic_name].index.name = None\n",
        "    ds[omic_name] = ds[omic_name].T\n",
        "    # N.B.: the matrices have been transposed so that now we have samples as rows and features as columns\n",
        "\n",
        "\n",
        "y = pd.read_csv('clusters.txt', sep='\\t').set_index('subjects')    # this will contain the true cluster label of each sample\n",
        "true_cluster_labels_synthetic = y.values.reshape(y.shape[0])-1"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ipWxcNChNf9B"
      },
      "source": [
        "Pre-process omics"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "l4lc2486NhLj"
      },
      "source": [
        "for omic in omics:\n",
        "    if omic == \"meth\":\n",
        "        mu = 0\n",
        "        sigma = 0.4 # 0.4\n",
        "\n",
        "    elif omic == \"prot\":\n",
        "        mu = 0\n",
        "        sigma = 4 # 4\n",
        "\n",
        "    else:\n",
        "        mu = 0\n",
        "        sigma = 4 # 4\n",
        "        \n",
        "    n_samples, n_features = ds[omic].shape\n",
        "\n",
        "    # Add salt & pepper noise\n",
        "    np.random.seed(42)   # fixed seed for reproducibility\n",
        "    noise = np.random.normal(mu, sigma, size=(n_samples, n_features))\n",
        "    p = 0.5 #0.8     # salt & pepper noise: add (gaussian) noise to a particolar feature of a particular sample with probability p\n",
        "    raveled_indices = np.random.choice(np.arange(n_samples*n_features), replace=False, size=int(n_samples*n_features*(1-p)))\n",
        "    indices = np.unravel_index(raveled_indices, (n_samples, n_features))\n",
        "    noise[indices] = 0\n",
        "\n",
        "    print(f\"Salt & pepper gaussian noise N({mu},{sigma**2}) is added to the {omic} dataset\")\n",
        "    ds[f'{omic}_noisy'] = ds[f'{omic}'] + noise\n",
        "\n",
        "\n",
        "\n",
        "    # Normalize omic\n",
        "    ds[f'{omic}_normalized'] = MinMaxScaler().fit_transform(ds[f'{omic}_noisy'])\n",
        "    ds[f'{omic}_normalized_no_noise'] = MinMaxScaler().fit_transform(ds[omic])  # version without noise"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "9_rSX5ZSNkEC"
      },
      "source": [
        "Reduce the n. of feature of each omic to 128 with PCA"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "nJigA9HFNjyy"
      },
      "source": [
        "for omic in omics:\n",
        "    pca = PCA(128)\n",
        "    # Reduce n. features\n",
        "    principalComponents_df = pd.DataFrame(pca.fit_transform(ds[f'{omic}_normalized']))\n",
        "    # Min-max normalize again\n",
        "    #principalComponents_df = pd.DataFrame(MinMaxScaler().fit_transform(principalComponents_df))\n",
        "    # Save reduced omic\n",
        "    principalComponents_df.to_csv(f'{omic}_synthetic_128.txt', sep=',', header=False)\n",
        "\n",
        "\n",
        "for omic in omics:\n",
        "    pca = PCA(128)\n",
        "    # Reduce n. features\n",
        "    principalComponents_df = pd.DataFrame(pca.fit_transform(ds[f'{omic}_normalized_no_noise']))\n",
        "    # Min-max normalize again\n",
        "    #principalComponents_df = pd.DataFrame(MinMaxScaler().fit_transform(principalComponents_df))\n",
        "    # Save reduced omic\n",
        "    principalComponents_df.to_csv(f'{omic}_synthetic_no_noise_128.txt', sep=',', header=False)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "3nOMxkyCNR07"
      },
      "source": [
        "## Lung"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "e3drZfGeNW6g"
      },
      "source": [
        "Load omics"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "7ZxnN1yJNS8L"
      },
      "source": [
        "# Load dataset\n",
        "!gdown 'https://drive.google.com/uc?id=1kNQxRoCs6TIGVzdlpEpUkAqx8FyoSrsi' -O ./mRNA.txt\n",
        "!gdown 'https://drive.google.com/uc?id=16V4tY8GeCUemOrf_KBqJwdixrQYZmFHv' -O ./miRNA.txt\n",
        "!gdown 'https://drive.google.com/uc?id=1UuCNcXbxHS1lW3bRkixsbo77xIFsWxIG' -O ./meth.txt\n",
        "!gdown 'https://drive.google.com/uc?id=1RT1EiQCj19gGD8bYr0UBjxvUPl8473FX' -O ./cnv.txt  # corrected for batch effects\n",
        "!gdown 'https://drive.google.com/uc?id=1t1rTrZNVHqZRo_F0YA0heD3WbvijHox7' -O ./clusters.txt\n",
        "\n",
        "ds = {}     # this will contain each omic\n",
        "omics = ['mRNA','miRNA','meth','cnv']\n",
        "for omic_name in omics:\n",
        "    path = omic_name + \".txt\"\n",
        "    if omic_name not in ds:\n",
        "        ds[omic_name] = pd.read_csv(path, sep='\\t', index_col=0)\n",
        "\n",
        "y = pd.read_csv('clusters.txt', sep='\\t', index_col=0)\n",
        "true_cluster_labels_lung = y.values.reshape(y.shape[0])"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "lA4R-m6YRMJV"
      },
      "source": [
        "Pre-process omics"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "G9jVCWyJRNq-"
      },
      "source": [
        "# mRNA\n",
        "# Keep only protein coding genes\n",
        "# NB: the file idx_mRNA_prot_cod.txt, containing the protein coding genes of the\n",
        "# human genome, has been obtained in the colab notebook 'Find protein coding genes'\n",
        "!gdown 'https://drive.google.com/uc?id=1i6lXPIh60BtCF-ujfaqf4WrvWM4XMvNF' -O ./idx_mRNA_prot_cod_lung.txt\n",
        "\n",
        "idx_mRNA_prot_cod = pd.read_csv('idx_mRNA_prot_cod_lung.txt')\n",
        "idx_mRNA_prot_cod = idx_mRNA_prot_cod['idx'].values\n",
        "\n",
        "ds['mRNA'] = ds['mRNA'].iloc[:, idx_mRNA_prot_cod]\n",
        "\n",
        "# Delete genes with a zero expression value across all the samples\n",
        "ds['mRNA'] = ds['mRNA'].loc[:, (ds['mRNA'] != 0).any(axis=0)]\n",
        "\n",
        "# Normalize mRNA with MinMax Scaler\n",
        "ds['mRNA_normalized'] = MinMaxScaler().fit_transform(ds['mRNA'].values)\n",
        "\n",
        "\n",
        "\n",
        "# miRNA\n",
        "# Delete sequences with a zero expression value across all the samples\n",
        "ds['miRNA'] = ds['miRNA'].loc[:, (ds['miRNA'] != 0).any(axis=0)]\n",
        "\n",
        "# Normalize with log2 normalization\n",
        "ds['miRNA'] = np.log(ds['miRNA'] + 1) / np.log(2)\n",
        "\n",
        "# Normalize with MinMaxScaler\n",
        "ds['miRNA_normalized'] = MinMaxScaler().fit_transform(ds['miRNA'].values)\n",
        "\n",
        "\n",
        "\n",
        "# meth\n",
        "# Delete sequences with a zero expression value across all the samples\n",
        "ds['meth'] = ds['meth'].loc[:, (ds['meth'] != 0).any(axis=0)]\n",
        "\n",
        "# Normalize with MinMaxScaler\n",
        "ds['meth_normalized'] = MinMaxScaler().fit_transform(ds['meth'].values)\n",
        "\n",
        "\n",
        "\n",
        "# cnv\n",
        "# Normalize with MinMaxScaler\n",
        "ds['cnv_normalized'] = MinMaxScaler().fit_transform(ds['cnv'].values)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "AukNl_GHROM2"
      },
      "source": [
        "Reduce the n. of feature of each omic to 256 with PCA"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "niDTEbgwRO8c"
      },
      "source": [
        "for omic in omics:\n",
        "    pca = PCA(256)\n",
        "    # Reduce n. features\n",
        "    principalComponents_df = pd.DataFrame(pca.fit_transform(ds[f'{omic}_normalized']))\n",
        "    # Min-max normalize again\n",
        "    #principalComponents_df = pd.DataFrame(MinMaxScaler().fit_transform(principalComponents_df))\n",
        "    # Save reduced omic\n",
        "    principalComponents_df.to_csv(f'{omic}_lung_256.txt', sep=',', header=False)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "PFTqqP3ENTUp"
      },
      "source": [
        "## Kidney"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "RKuRYG0eNXUs"
      },
      "source": [
        "Load omics"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "hLPX0rpyNUiA"
      },
      "source": [
        "# Load dataset\n",
        "!gdown 'https://drive.google.com/uc?id=1i1do_UTzwXzPVIDDmYSFJEholK2Mp8g_' -O ./mRNA.txt\n",
        "!gdown 'https://drive.google.com/uc?id=1liKeOBKjnbCi1CIjcOPA3Zxv2fRzCfa2' -O ./miRNA.txt\n",
        "!gdown 'https://drive.google.com/uc?id=1qr9joY0bAVDLvjWsKF5xf3CaRolBu-mP' -O ./meth.txt\n",
        "!gdown 'https://drive.google.com/uc?id=1R-U2iDgM4oEyzNRfBIA2kXMbKw_s0QtI' -O ./clusters.txt\n",
        "\n",
        "ds = {}     # this will contain each omic\n",
        "omics = ['mRNA','miRNA','meth']\n",
        "for omic_name in omics:\n",
        "    path = omic_name + \".txt\"\n",
        "    if omic_name not in ds:\n",
        "        ds[omic_name] = pd.read_csv(path, sep='\\t', index_col=0)\n",
        "\n",
        "y = pd.read_csv('clusters.txt', sep='\\t', index_col=0)\n",
        "true_cluster_labels_kidney = y.values.reshape(y.shape[0])"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "88MHWyhVRSTj"
      },
      "source": [
        "Pre-process omics"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "wbTfK0ZXRSqT"
      },
      "source": [
        "# mRNA\n",
        "# Keep only protein coding genes\n",
        "# NB: the file idx_mRNA_prot_cod.txt, containing the protein coding genes of the\n",
        "# human genome, has been obtained in the colab notebook 'Find protein coding genes'\n",
        "if not os.path.exists(\"./idx_mRNA_prot_cod_kidney.txt\"):\n",
        "    !gdown 'https://drive.google.com/uc?id=1Pi4u8y_YAc2tmOWZYaeLn9wGdzu4cFC5' -O ./idx_mRNA_prot_cod_kidney.txt\n",
        "\n",
        "idx_mRNA_prot_cod = pd.read_csv('idx_mRNA_prot_cod_kidney.txt')\n",
        "idx_mRNA_prot_cod = idx_mRNA_prot_cod['idx'].values\n",
        "\n",
        "ds['mRNA'] = ds['mRNA'].iloc[:, idx_mRNA_prot_cod]\n",
        "\n",
        "# Delete genes with a zero expression value across all the samples\n",
        "ds['mRNA'] = ds['mRNA'].loc[:, (ds['mRNA'] != 0).any(axis=0)]\n",
        "\n",
        "# Normalize mRNA with MinMax Scaler\n",
        "ds['mRNA_normalized'] = MinMaxScaler().fit_transform(ds['mRNA'].values)\n",
        "\n",
        "\n",
        "\n",
        "# miRNA\n",
        "# Delete sequences with a zero expression value across all the samples\n",
        "ds['miRNA'] = ds['miRNA'].loc[:, (ds['miRNA'] != 0).any(axis=0)]\n",
        "\n",
        "# Normalize with log2 normalization\n",
        "ds['miRNA'] = np.log(ds['miRNA'] + 1) / np.log(2)\n",
        "\n",
        "# Normalize with MinMaxScaler\n",
        "ds['miRNA_normalized'] = MinMaxScaler().fit_transform(ds['miRNA'].values)\n",
        "\n",
        "\n",
        "\n",
        "# meth\n",
        "# Delete sequences with a zero expression value across all the samples\n",
        "ds['meth'] = ds['meth'].loc[:, (ds['meth'] != 0).any(axis=0)]\n",
        "\n",
        "# Normalize with MinMaxScaler\n",
        "ds['meth_normalized'] = MinMaxScaler().fit_transform(ds['meth'].values)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "t39SXOjBRTDb"
      },
      "source": [
        "Reduce the n. of feature of each omic to 256 with PCA"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "K2GP8OdGRTVM"
      },
      "source": [
        "for omic in omics:\n",
        "    pca = PCA(256)\n",
        "    # Reduce n. features\n",
        "    principalComponents_df = pd.DataFrame(pca.fit_transform(ds[f'{omic}_normalized']))\n",
        "    # Min-max normalize again\n",
        "    #principalComponents_df = pd.DataFrame(MinMaxScaler().fit_transform(principalComponents_df))\n",
        "    # Save reduced omic\n",
        "    principalComponents_df.to_csv(f'{omic}_kidney_256.txt', sep=',', header=False)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "1DyZdamhPNtw"
      },
      "source": [
        "# Import R libraries and rpy"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "1QC5f4hZ1Eug"
      },
      "source": [
        "%load_ext rpy2.ipython"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "-L_OJRyD1URB"
      },
      "source": [
        "%%R\n",
        "if (!requireNamespace(\"BiocManager\", quietly = TRUE))\n",
        "    install.packages(\"BiocManager\")\n",
        "BiocManager::install(\"iClusterPlus\")\n",
        "library('iClusterPlus')\n",
        "\n",
        "#install.packages('googledrive')\n",
        "#library('googledrive')\n",
        "#drive_deauth()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "DUUuDZCWPUAA"
      },
      "source": [
        "# Apply iClusterPlus"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "z_-FvOKW2c8Z"
      },
      "source": [
        "## Synthetic (noisy)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "C_-0yPWpPdl_"
      },
      "source": [
        "Load omics and adapt them to iClusterPlus accepted format"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "rbBwn4gh5YC9"
      },
      "source": [
        "%%R\n",
        "print('mRNA syntethic 128 loading...')\n",
        "mRNA = read.csv(\"mRNA_synthetic_128.txt\", sep=\",\", header=FALSE)\n",
        "mRNA <- as.matrix(mRNA) \n",
        "mRNA <- mRNA[,-1]\n",
        "colnames(mRNA) <- NULL\n",
        "\n",
        "print('meth syntethic 128 loading...')\n",
        "meth = read.csv(\"meth_synthetic_128.txt\", sep=\",\", header=FALSE)\n",
        "meth <- as.matrix(meth)\n",
        "meth <- meth[,-1]\n",
        "colnames(meth) <- NULL\n",
        "\n",
        "print('prot syntethic 128 loading...')\n",
        "prot = read.csv(\"prot_synthetic_128.txt\", sep=\",\", header=FALSE)\n",
        "prot <- as.matrix(prot)\n",
        "prot <- prot[,-1]\n",
        "colnames(prot) <- NULL"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "thya9MW3Qynq"
      },
      "source": [
        "Apply iClusterPlus"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "UKvKl3HWQ0Ct"
      },
      "source": [
        "%%R\n",
        "# run iClusterPlus algorithm\n",
        "r.icluster <- iClusterPlus::iClusterPlus(\n",
        "  mRNA, # Providing each omics type\n",
        "  meth,\n",
        "  prot,\n",
        "  type=c(\"gaussian\", \"gaussian\", \"gaussian\"), # Providing the distributions\n",
        "  K=64, # provide the number of factors to learn\n",
        "  alpha=c(1,1,1), # as well as other model parameters\n",
        "  lambda=c(.03,.03,.03)\n",
        ")\n",
        "\n",
        "# extract the H and W matrices from the run result\n",
        "# here, we refer to H as z, to keep with iCluster terminology\n",
        "# NB: H is basically the encoded dataset: (500,dim_ls) where dim_ls = K = 64\n",
        "icluster.z <- r.icluster$meanZ\n",
        "#icluster.ws <- r.icluster$beta\n",
        "\n",
        "# Save H matrix\n",
        "write.table(icluster.z, file=\"icluster_z_mat_synthetic.txt\", row.names=FALSE, col.names=FALSE, sep=',')"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "I3rOEZBpXyWy"
      },
      "source": [
        "Results"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Yea_K5N8XxFE"
      },
      "source": [
        "# Load H matrix (iClusterPlus-encoded-integrated dataset)\n",
        "H_synthetic = pd.read_csv(\"icluster_z_mat_synthetic.txt\", sep=',', header=None)\n",
        "\n",
        "for k in range(2,11):\n",
        "    # Apply K-means\n",
        "    kmeans = KMeans(n_clusters=k, random_state=0)\n",
        "    cluster_labels = kmeans.fit_predict(H_synthetic)\n",
        "\n",
        "    # Perform a 2D PCA to visualize the dataset\n",
        "    pca = PCA(2)\n",
        "    principalComponents = pca.fit_transform(H_synthetic)\n",
        "    kmeans.cluster_centers_ = pca.transform(kmeans.cluster_centers_)\n",
        "\n",
        "    # Plot the clustered dataset with cluster assignments and true cluster labels\n",
        "    plot_2D_dataset(principalComponents, cluster_labels, title=f'Dataset visualization', caption=f'predicted {k} clusters')\n",
        "    plot_2D_dataset(principalComponents, true_cluster_labels_synthetic, title=f'Dataset visualization', caption='true clusters')\n",
        "\n",
        "    # Plot the confusion matrix\n",
        "    plot_confusion_matrix(true_cluster_labels_synthetic, cluster_labels)\n",
        "\n",
        "    # Compute silhouette on the iClusterPlus-encoded dataset with cluster assignments and true cluster labels\n",
        "    print(f\"Silhouette, predicted clusters: {silhouette_score(H_synthetic, cluster_labels)}\")\n",
        "    print(f\"Silhouette, true clusters: {silhouette_score(H_synthetic, true_cluster_labels_synthetic)}\")\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "QpH0nfEFqngN"
      },
      "source": [
        "We can clearly see from the visualized dataset that iClusterPlus is not robust to noise in the input data"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "qoi7sZ5EqDpF"
      },
      "source": [
        "## Synthetic (no noise)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "iyd_mFMAqGSV"
      },
      "source": [
        "Load omics and adapt them to iClusterPlus accepted format"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "EX4Lme_DqGll"
      },
      "source": [
        "%%R\n",
        "print('mRNA syntethic no noise 128 loading...')\n",
        "mRNA = read.csv(\"mRNA_synthetic_no_noise_128.txt\", sep=\",\", header=FALSE)\n",
        "mRNA <- as.matrix(mRNA) \n",
        "mRNA <- mRNA[,-1]\n",
        "colnames(mRNA) <- NULL\n",
        "\n",
        "print('meth syntethic no noise 128 loading...')\n",
        "meth = read.csv(\"meth_synthetic_no_noise_128.txt\", sep=\",\", header=FALSE)\n",
        "meth <- as.matrix(meth)\n",
        "meth <- meth[,-1]\n",
        "colnames(meth) <- NULL\n",
        "\n",
        "print('prot syntethic no noise 128 loading...')\n",
        "prot = read.csv(\"prot_synthetic_no_noise_128.txt\", sep=\",\", header=FALSE)\n",
        "prot <- as.matrix(prot)\n",
        "prot <- prot[,-1]\n",
        "colnames(prot) <- NULL"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "6aVQxKOxqGy8"
      },
      "source": [
        "Apply iClusterPlus"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "H2Ws4pr3qHFI"
      },
      "source": [
        "%%R\n",
        "# run iClusterPlus algorithm\n",
        "r.icluster <- iClusterPlus::iClusterPlus(\n",
        "  mRNA, # Providing each omics type\n",
        "  meth,\n",
        "  prot,\n",
        "  type=c(\"gaussian\", \"gaussian\", \"gaussian\"), # Providing the distributions\n",
        "  K=64, # provide the number of factors to learn\n",
        "  alpha=c(1,1,1), # as well as other model parameters\n",
        "  lambda=c(.03,.03,.03)\n",
        ")\n",
        "\n",
        "# extract the H and W matrices from the run result\n",
        "# here, we refer to H as z, to keep with iCluster terminology\n",
        "# NB: H is basically the encoded dataset: (500,dim_ls) where dim_ls = K = 64\n",
        "icluster.z <- r.icluster$meanZ\n",
        "#icluster.ws <- r.icluster$beta\n",
        "\n",
        "# Save H matrix\n",
        "write.table(icluster.z, file=\"icluster_z_mat_synthetic_no_noise.txt\", row.names=FALSE, col.names=FALSE, sep=',')"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "OLrMNrJoqHZX"
      },
      "source": [
        "Results"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "9t2N9iAlqHps"
      },
      "source": [
        "# Load H matrix (iClusterPlus-encoded-integrated dataset)\n",
        "H_synthetic = pd.read_csv(\"icluster_z_mat_synthetic_no_noise.txt\", sep=',', header=None)\n",
        "\n",
        "for k in range(2,11):\n",
        "    # Apply K-means\n",
        "    kmeans = KMeans(n_clusters=k, random_state=0)\n",
        "    cluster_labels = kmeans.fit_predict(H_synthetic)\n",
        "\n",
        "    # Perform a 2D PCA to visualize the dataset\n",
        "    pca = PCA(2)\n",
        "    principalComponents = pca.fit_transform(H_synthetic)\n",
        "    kmeans.cluster_centers_ = pca.transform(kmeans.cluster_centers_)\n",
        "\n",
        "    # Plot the clustered dataset with cluster assignments and true cluster labels\n",
        "    plot_2D_dataset(principalComponents, cluster_labels, title=f'Dataset visualization', caption=f'predicted {k} clusters')\n",
        "    plot_2D_dataset(principalComponents, true_cluster_labels_synthetic, title=f'Dataset visualization', caption='true clusters')\n",
        "\n",
        "    # Plot the confusion matrix\n",
        "    plot_confusion_matrix(true_cluster_labels_synthetic, cluster_labels)\n",
        "\n",
        "    # Compute silhouette on the iClusterPlus-encoded dataset with cluster assignments and true cluster labels\n",
        "    print(f\"Silhouette, predicted clusters: {silhouette_score(H_synthetic, cluster_labels)}\")\n",
        "    print(f\"Silhouette, true clusters: {silhouette_score(H_synthetic, true_cluster_labels_synthetic)}\")\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ouixXCzU2eBv"
      },
      "source": [
        "## Lung"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Pp2TElr3TyFb"
      },
      "source": [
        "Load omics and adapt them to iClusterPlus accepted format"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "CMjOWlGb1-oB"
      },
      "source": [
        "%%R\n",
        "print('mRNA lung 256 loading...')\n",
        "mRNA = read.csv(\"mRNA_lung_256.txt\", sep=\",\", header=FALSE)\n",
        "mRNA <- as.matrix(mRNA)\n",
        "mRNA <- mRNA[,-1]\n",
        "colnames(mRNA) <- NULL\n",
        "\n",
        "print('miRNA lung 256 loading...')\n",
        "miRNA = read.csv(\"miRNA_lung_256.txt\", sep=\",\", header=FALSE)\n",
        "miRNA <- as.matrix(miRNA) \n",
        "miRNA <- miRNA[,-1]\n",
        "colnames(miRNA) <- NULL\n",
        "\n",
        "print('meth lung 256 loading...')\n",
        "meth = read.csv(\"meth_lung_256.txt\", sep=\",\", header=FALSE)\n",
        "meth <- as.matrix(meth)\n",
        "meth <- meth[,-1]\n",
        "colnames(meth) <- NULL\n",
        "\n",
        "print('cnv lung 256 loading...')\n",
        "cnv = read.csv(\"cnv_lung_256.txt\", sep=\",\", header=FALSE)\n",
        "cnv <- as.matrix(cnv) \n",
        "cnv <- cnv[,-1]\n",
        "colnames(cnv) <- NULL"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "g4Qvxo69gK6j"
      },
      "source": [
        "Apply iClusterPlus"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "yFhIUKeq3J6k"
      },
      "source": [
        "%%R\n",
        "# run iClusterPlus algorithm\n",
        "r.icluster <- iClusterPlus::iClusterPlus(\n",
        "  mRNA, # Providing each omics type\n",
        "  miRNA,\n",
        "  meth,\n",
        "  cnv,\n",
        "  type=c(\"gaussian\", \"gaussian\", \"gaussian\", \"gaussian\"), # Providing the distributions\n",
        "  K=64, # provide the number of factors to learn\n",
        "  alpha=c(1,1,1,1), # as well as other model parameters\n",
        "  lambda=c(.03,.03,.03,.03)\n",
        ")\n",
        "\n",
        "# extract the H and W matrices from the run result\n",
        "# here, we refer to H as z, to keep with iCluster terminology\n",
        "# NB: H is basically the encoded dataset: (783,dim_ls) where dim_ls = K = 64\n",
        "icluster.z <- r.icluster$meanZ\n",
        "#icluster.ws <- r.icluster$beta\n",
        "\n",
        "# Save H matrix\n",
        "write.table(icluster.z, file=\"icluster_z_mat_lung.txt\", row.names=FALSE, col.names=FALSE, sep=',')"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "1OIn129ygO9W"
      },
      "source": [
        "Results"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "LlT9PzVV3cgh"
      },
      "source": [
        "# Load H matrix (iClusterPlus-encoded-integrated dataset)\n",
        "H_synthetic = pd.read_csv(\"icluster_z_mat_lung.txt\", sep=',', header=None)\n",
        "\n",
        "for k in range(2,11):\n",
        "    # Apply K-means\n",
        "    kmeans = KMeans(n_clusters=k, random_state=0)\n",
        "    cluster_labels = kmeans.fit_predict(H_synthetic)\n",
        "\n",
        "    # Perform a 2D PCA to visualize the dataset\n",
        "    pca = PCA(2)\n",
        "    principalComponents = pca.fit_transform(H_synthetic)\n",
        "    kmeans.cluster_centers_ = pca.transform(kmeans.cluster_centers_)\n",
        "\n",
        "    # Plot the clustered dataset with cluster assignments and true cluster labels\n",
        "    plot_2D_dataset(principalComponents, cluster_labels, title=f'Dataset visualization', caption=f'predicted {k} clusters')\n",
        "    plot_2D_dataset(principalComponents, true_cluster_labels_lung, title=f'Dataset visualization', caption='true clusters')\n",
        "\n",
        "    # Plot the confusion matrix\n",
        "    plot_confusion_matrix(true_cluster_labels_lung, cluster_labels)\n",
        "\n",
        "    # Compute silhouette on the iClusterPlus-encoded dataset with cluster assignments and true cluster labels\n",
        "    print(f\"Silhouette, predicted clusters: {silhouette_score(H_synthetic, cluster_labels)}\")\n",
        "    print(f\"Silhouette, true clusters: {silhouette_score(H_synthetic, true_cluster_labels_lung)}\")\n",
        "\n",
        "\n",
        "    print()\n",
        "    print()\n",
        "    print()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "O5snNwSX1PzR"
      },
      "source": [
        "## Kidney"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "pBpzSnjnTwRj"
      },
      "source": [
        "Load omics and adapt them to iClusterPlus accepted format"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "RYUV8wCSqxoT"
      },
      "source": [
        "%%R\n",
        "print('mRNA kidney 256 loading...')\n",
        "mRNA = read.csv(\"mRNA_kidney_256.txt\", sep=\",\", header=FALSE)\n",
        "mRNA <- as.matrix(mRNA) \n",
        "mRNA <- mRNA[,-1]\n",
        "colnames(mRNA) <- NULL\n",
        "\n",
        "print('miRNA kidney 256 loading...')\n",
        "miRNA = read.csv(\"miRNA_kidney_256.txt\", sep=\",\", header=FALSE)\n",
        "miRNA <- as.matrix(miRNA) \n",
        "miRNA <- miRNA[,-1]\n",
        "colnames(miRNA) <- NULL\n",
        "\n",
        "print('meth kidney 256 loading...')\n",
        "meth = read.csv(\"meth_kidney_256.txt\", sep=\",\", header=FALSE)\n",
        "meth <- as.matrix(meth)\n",
        "meth <- meth[,-1]\n",
        "colnames(meth) <- NULL"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "JefGc9-oUJ8N"
      },
      "source": [
        "Apply iClusterPlus"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "12-sNZN9UK_X"
      },
      "source": [
        "%%R\n",
        "# run iClusterPlus algorithm\n",
        "r.icluster <- iClusterPlus::iClusterPlus(\n",
        "  mRNA, # Providing each omics type\n",
        "  miRNA,\n",
        "  meth,\n",
        "  type=c(\"gaussian\", \"gaussian\", \"gaussian\"), # Providing the distributions\n",
        "  K=64, # provide the number of factors to learn\n",
        "  alpha=c(1,1,1), # as well as other model parameters\n",
        "  lambda=c(.03,.03,.03)\n",
        ")\n",
        "\n",
        "# extract the H and W matrices from the run result\n",
        "# here, we refer to H as z, to keep with iCluster terminology\n",
        "# NB: H is basically the encoded dataset: (783,dim_ls) where dim_ls = K = 64\n",
        "icluster.z <- r.icluster$meanZ\n",
        "#icluster.ws <- r.icluster$beta\n",
        "\n",
        "# Save H matrix\n",
        "write.table(icluster.z, file=\"icluster_z_mat_kidney.txt\", row.names=FALSE, col.names=FALSE, sep=',')"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "KEJXszVtgQm-"
      },
      "source": [
        "Results"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "J-LSONEEgRFO"
      },
      "source": [
        "# Load H matrix (iClusterPlus-encoded-integrated dataset)\n",
        "H_synthetic = pd.read_csv(\"icluster_z_mat_kidney.txt\", sep=',', header=None)\n",
        "\n",
        "# Apply K-means\n",
        "for k in range(2,11):\n",
        "    kmeans = KMeans(n_clusters=k, random_state=0)\n",
        "    cluster_labels = kmeans.fit_predict(H_synthetic)\n",
        "\n",
        "    # Perform a 2D PCA to visualize the dataset\n",
        "    pca = PCA(2)\n",
        "    principalComponents = pca.fit_transform(H_synthetic)\n",
        "    kmeans.cluster_centers_ = pca.transform(kmeans.cluster_centers_)\n",
        "\n",
        "    # Plot the clustered dataset with cluster assignments and true cluster labels\n",
        "    plot_2D_dataset(principalComponents, cluster_labels, title=f'Dataset visualization', caption=f'predicted {k} clusters')\n",
        "    plot_2D_dataset(principalComponents, true_cluster_labels_kidney, title=f'Dataset visualization', caption='true clusters')\n",
        "\n",
        "    # Plot the confusion matrix\n",
        "    plot_confusion_matrix(true_cluster_labels_kidney, cluster_labels)\n",
        "\n",
        "    # Compute silhouette on the iClusterPlus-encoded dataset with cluster assignments and true cluster labels\n",
        "    print(f\"Silhouette, predicted clusters: {silhouette_score(H_synthetic, cluster_labels)}\")\n",
        "    print(f\"Silhouette, true clusters: {silhouette_score(H_synthetic, true_cluster_labels_kidney)}\")"
      ],
      "execution_count": null,
      "outputs": []
    }
  ]
}