Seurat_EDA_e_integracion_de_datos.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "provenance": [],
      "authorship_tag": "ABX9TyPMzOPqHennnl7cGf+0yLZD",
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "ir",
      "display_name": "R"
    },
    "language_info": {
      "name": "R"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/martings/My_own_single_cell_toolkit/blob/main/Seurat_EDA_e_integracion_de_datos.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "#En colab tarda unos 25 minutos en instalarse el Seurat\n",
        "install.packages('Seurat')\n",
        "#install.packages('tidyverse')"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "GqN0LugyVHru",
        "outputId": "43192de0-8ad1-4a92-e12d-dbfe6f5c4c86"
      },
      "execution_count": 2,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stderr",
          "text": [
            "Installing package into ‘/usr/local/lib/R/site-library’\n",
            "(as ‘lib’ is unspecified)\n",
            "\n",
            "also installing the dependencies ‘bitops’, ‘gtools’, ‘caTools’, ‘plyr’, ‘tensor’, ‘BH’, ‘sitmo’, ‘globals’, ‘listenv’, ‘parallelly’, ‘zoo’, ‘htmlwidgets’, ‘lazyeval’, ‘crosstalk’, ‘promises’, ‘RcppTOML’, ‘here’, ‘gplots’, ‘reshape2’, ‘gridExtra’, ‘RcppArmadillo’, ‘sp’, ‘httpuv’, ‘xtable’, ‘sourcetools’, ‘later’, ‘spatstat.data’, ‘spatstat.random’, ‘spatstat.utils’, ‘spatstat.sparse’, ‘goftest’, ‘abind’, ‘deldir’, ‘polyclip’, ‘FNN’, ‘dqrng’, ‘cowplot’, ‘fitdistrplus’, ‘future’, ‘future.apply’, ‘ggrepel’, ‘ggridges’, ‘ica’, ‘igraph’, ‘irlba’, ‘leiden’, ‘lmtest’, ‘matrixStats’, ‘miniUI’, ‘patchwork’, ‘pbapply’, ‘plotly’, ‘png’, ‘progressr’, ‘RANN’, ‘Rcpp’, ‘RcppAnnoy’, ‘reticulate’, ‘ROCR’, ‘Rtsne’, ‘scattermore’, ‘sctransform’, ‘SeuratObject’, ‘shiny’, ‘spatstat.explore’, ‘spatstat.geom’, ‘uwot’, ‘RcppEigen’, ‘RcppProgress’\n",
            "\n",
            "\n",
            "The legacy packages maptools, rgdal, and rgeos, underpinning this package\n",
            "will retire shortly. Please refer to R-spatial evolution reports on\n",
            "https://r-spatial.org/r/2023/05/15/evolution4.html for details.\n",
            "This package is now running under evolution status 0 \n",
            "\n",
            "Attaching SeuratObject\n",
            "\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "library(Seurat)\n",
        "suppressMessages({\n",
        "    library(tidyverse)\n",
        "})"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "AVH1Uo53igaP",
        "outputId": "6a42b58f-6b7d-43c6-fc8b-628967ed00f4"
      },
      "execution_count": 2,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stderr",
          "text": [
            "The legacy packages maptools, rgdal, and rgeos, underpinning this package\n",
            "will retire shortly. Please refer to R-spatial evolution reports on\n",
            "https://r-spatial.org/r/2023/05/15/evolution4.html for details.\n",
            "This package is now running under evolution status 0 \n",
            "\n",
            "Attaching SeuratObject\n",
            "\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "#Analizo uno de los archivos a integrar para ver que pinta tienen\n",
        "AdultMG_Visvader<-readRDS(\"/content/SC_MG_Data/Visvader.10x.combined.BackUP.dim30.rds\")"
      ],
      "metadata": {
        "id": "EYe-zQgEZAlA"
      },
      "execution_count": 8,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "pbmc <- SCTransform(AdultMG_Visvader)\n",
        "rm(AdultMG_Visvader)\n",
        "pbmc <- RunPCA(pbmc)\n",
        "pbmc <- RunUMAP(pbmc, dims = 1:20)\n",
        "pbmc <- FindNeighbors(pbmc, dims = 1:20)\n",
        "pbmc <- FindClusters(pbmc, resolution = 0.1 )\n",
        "DimPlot(pbmc, reduction = \"umap\")"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "nxuFRPYBf5Im",
        "outputId": "b8fa4dfd-8319-41d9-b278-cf76089724d0"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stderr",
          "text": [
            "Calculating cell attributes from input UMI matrix: log_umi\n",
            "\n",
            "Variance stabilizing transformation of count matrix of size 13224 by 22287\n",
            "\n",
            "Model formula is y ~ log_umi\n",
            "\n",
            "Get Negative Binomial regression parameters per gene\n",
            "\n",
            "Using 2000 genes, 5000 cells\n",
            "\n"
          ]
        },
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "  |======================================================================| 100%\n"
          ]
        },
        {
          "output_type": "stream",
          "name": "stderr",
          "text": [
            "Found 26 outliers - those will be ignored in fitting/regularization step\n",
            "\n",
            "\n",
            "Second step: Get residuals using fitted parameters for 13224 genes\n",
            "\n"
          ]
        },
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "  |======================================================================| 100%\n"
          ]
        },
        {
          "output_type": "stream",
          "name": "stderr",
          "text": [
            "Computing corrected count matrix for 13224 genes\n",
            "\n"
          ]
        },
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "  |======================================================================| 100%\n"
          ]
        },
        {
          "output_type": "stream",
          "name": "stderr",
          "text": [
            "Calculating gene attributes\n",
            "\n",
            "Wall clock passed: Time difference of 2.256748 mins\n",
            "\n",
            "Determine variable features\n",
            "\n",
            "Place corrected count matrix in counts slot\n",
            "\n",
            "Centering data matrix\n",
            "\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 7,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 135
        },
        "id": "tjmEEpU-U944",
        "outputId": "2a44743e-afe4-46c5-ab70-6bdaae91e61a"
      },
      "outputs": [
        {
          "output_type": "error",
          "ename": "ERROR",
          "evalue": "ignored",
          "traceback": [
            "Error in readRDS(\"/content/SC_MG_Data/AdultMG_combined_BackUp.rds\"): error reading from connection\nTraceback:\n",
            "1. readRDS(\"/content/SC_MG_Data/AdultMG_combined_BackUp.rds\")"
          ]
        }
      ],
      "source": [
        "#Cargo los dataset de entrada para realizar la integracion\n",
        "\n",
        "Khaled_total<-readRDS(\"/content/SC_MG_Data/Khaled.10x.combined.estricto.features.BackUP.dim30.rds\")\n",
        "AdultMG_Spike<-readRDS(\"/content/SC_MG_Data/AdultMG_combined_BackUp.rds\")\n",
        "AdultMG_Visvader<-readRDS(\"/content/SC_MG_Data/Visvader.10x.combined.BackUP.dim30.rds\")\n",
        "\n",
        "#Armo una lista con todos los objetos a integrar\n",
        "seurat_object_list<-c(\"Khaled_total\",\"AdultMG_Spike\",\"AdultMG_Visvader\")\n",
        "pbmc.list <- list()\n",
        "for (i in seurat_object_list) {\n",
        "  pbmc.list[i] <- list(get(i))\n",
        "}\n",
        "#Borro los archivos originales de memoria RAM\n",
        "rm(Khaled_total)\n",
        "rm(AdultMG_Spike)\n",
        "rm(AdultMG_Visvader)\n",
        "\n",
        "#Transformo y normalizo los datos a integrar\n",
        "for (i in names(pbmc.list)) {\n",
        "  pbmc.list[[i]] <- SCTransform(pbmc.list[[i]], verbose = T)\n",
        "}\n",
        "\n",
        "#Aca busco los features (genes) que van a ser usados para realizar la integracion. \n",
        "#Busca los genes que tienen buena expresion y buen peso estadistico entre los tres data set como para usarlos de featues principales.\n",
        "#El numero de features, va a mejorar la calidad/resolucion de la integracion, pero va a requerir mas RAM para realizar la integracion\n",
        "pbmc.features <- SelectIntegrationFeatures(object.list = pbmc.list, nfeatures = 3000)\n",
        "\n",
        "#Realiza la normalizacion SCT a los datos a integrar\n",
        "pbmc.list <- PrepSCTIntegration(object.list = pbmc.list, anchor.features = pbmc.features)\n",
        "\n",
        "#Se puede usar un dataset de referencia, de mapa principal al cual el resto de los dataset van a tratar de amolgarse. \n",
        "#Es util por si queremos conservar las separaciones de alguno de los dataset antes que otros.\n",
        "#El programa puede correr sin referencia puntual o mismo puede hacerlo con varios. \n",
        "reference_dataset <- which(names(pbmc.list) == \"Khaled_total\")\n",
        "\n",
        "#\"\"Centra\"\" los datos de los tres dataset normalizados o tambien se puede pensar que los limita a un conjunto de features que contempla lo mejor posible los 3.\n",
        "pbmc.anchors <- FindIntegrationAnchors(object.list = pbmc.list, normalization.method = \"SCT\", \n",
        "                                       anchor.features = pbmc.features, reference = reference_dataset)\n",
        "#Borro por cuestiones de RAM la lista original\n",
        "rm(pbmc.list)\n",
        "#Con todos los datos centados y normalizados, los integro todos a un nuevo dataset\n",
        "pbmc.integrated <- IntegrateData(anchorset = pbmc.anchors, normalization.method = \"SCT\")\n",
        "#Es conveniente guardar el archivo crudo para tener a disposicion si \n",
        "saveRDS(pbmc.integrated, \"pbmc.integrated.backup.rds\")\n",
        "#Siempre es bueno ahorrar RAM\n",
        "rm(pbmc.anchors)\n",
        "\n",
        "#\n",
        "pbmc.integrated <- SCTransform(pbmc.integrated,nfeatures = 3000 )\n",
        "pbmc.integrated <- RunPCA(pbmc.integrated)\n",
        "pbmc.integrated <- RunUMAP(pbmc.integrated, dims = 1:20)\n",
        "pbmc.integrated <- FindNeighbors(pbmc.integrated, dims = 1:20)\n",
        "pbmc.integrated <- FindClusters(pbmc.integrated, resolution = 0.1 )\n",
        "DimPlot(pbmc.integrated, reduction = \"umap\", group.by = c(\"orig.ident\", \"Cell_Type\"))\n",
        "\n",
        "saveRDS(pbmc.integrated, \"pbmc.integrated.post.rds\")"
      ]
    }
  ]
}