spatialRNAseq_post_integration.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "360d6bb4-a72b-49b3-a9fe-e306e089d1fb",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "── \u001b[1mAttaching packages\u001b[22m ─────────────────────────────────────── tidyverse 1.3.1 ──\n",
      "\n",
      "\u001b[32m✔\u001b[39m \u001b[34mggplot2\u001b[39m 3.4.2     \u001b[32m✔\u001b[39m \u001b[34mpurrr  \u001b[39m 1.0.1\n",
      "\u001b[32m✔\u001b[39m \u001b[34mtibble \u001b[39m 3.2.1     \u001b[32m✔\u001b[39m \u001b[34mdplyr  \u001b[39m 1.1.1\n",
      "\u001b[32m✔\u001b[39m \u001b[34mtidyr  \u001b[39m 1.3.0     \u001b[32m✔\u001b[39m \u001b[34mstringr\u001b[39m 1.5.0\n",
      "\u001b[32m✔\u001b[39m \u001b[34mreadr  \u001b[39m 2.1.4     \u001b[32m✔\u001b[39m \u001b[34mforcats\u001b[39m 1.0.0\n",
      "\n",
      "── \u001b[1mConflicts\u001b[22m ────────────────────────────────────────── tidyverse_conflicts() ──\n",
      "\u001b[31m✖\u001b[39m \u001b[34mdplyr\u001b[39m::\u001b[32mfilter()\u001b[39m masks \u001b[34mstats\u001b[39m::filter()\n",
      "\u001b[31m✖\u001b[39m \u001b[34mdplyr\u001b[39m::\u001b[32mlag()\u001b[39m    masks \u001b[34mstats\u001b[39m::lag()\n",
      "\n",
      "Loading required package: Matrix\n",
      "\n",
      "\n",
      "Attaching package: ‘Matrix’\n",
      "\n",
      "\n",
      "The following objects are masked from ‘package:tidyr’:\n",
      "\n",
      "    expand, pack, unpack\n",
      "\n",
      "\n",
      "Loading required package: limma\n",
      "\n",
      "\n",
      "Attaching package: ‘reshape2’\n",
      "\n",
      "\n",
      "The following object is masked from ‘package:tidyr’:\n",
      "\n",
      "    smiths\n",
      "\n",
      "\n",
      "Loading required package: stats4\n",
      "\n",
      "Loading required package: BiocGenerics\n",
      "\n",
      "\n",
      "Attaching package: ‘BiocGenerics’\n",
      "\n",
      "\n",
      "The following object is masked from ‘package:limma’:\n",
      "\n",
      "    plotMA\n",
      "\n",
      "\n",
      "The following objects are masked from ‘package:dplyr’:\n",
      "\n",
      "    combine, intersect, setdiff, union\n",
      "\n",
      "\n",
      "The following objects are masked from ‘package:stats’:\n",
      "\n",
      "    IQR, mad, sd, var, xtabs\n",
      "\n",
      "\n",
      "The following objects are masked from ‘package:base’:\n",
      "\n",
      "    anyDuplicated, aperm, append, as.data.frame, basename, cbind,\n",
      "    colnames, dirname, do.call, duplicated, eval, evalq, Filter, Find,\n",
      "    get, grep, grepl, intersect, is.unsorted, lapply, Map, mapply,\n",
      "    match, mget, order, paste, pmax, pmax.int, pmin, pmin.int,\n",
      "    Position, rank, rbind, Reduce, rownames, sapply, setdiff, sort,\n",
      "    table, tapply, union, unique, unsplit, which.max, which.min\n",
      "\n",
      "\n",
      "\n",
      "Attaching package: ‘S4Vectors’\n",
      "\n",
      "\n",
      "The following objects are masked from ‘package:Matrix’:\n",
      "\n",
      "    expand, unname\n",
      "\n",
      "\n",
      "The following objects are masked from ‘package:dplyr’:\n",
      "\n",
      "    first, rename\n",
      "\n",
      "\n",
      "The following object is masked from ‘package:tidyr’:\n",
      "\n",
      "    expand\n",
      "\n",
      "\n",
      "The following objects are masked from ‘package:base’:\n",
      "\n",
      "    expand.grid, I, unname\n",
      "\n",
      "\n",
      "Loading required package: SummarizedExperiment\n",
      "\n",
      "Loading required package: MatrixGenerics\n",
      "\n",
      "Loading required package: matrixStats\n",
      "\n",
      "\n",
      "Attaching package: ‘matrixStats’\n",
      "\n",
      "\n",
      "The following object is masked from ‘package:dplyr’:\n",
      "\n",
      "    count\n",
      "\n",
      "\n",
      "\n",
      "Attaching package: ‘MatrixGenerics’\n",
      "\n",
      "\n",
      "The following objects are masked from ‘package:matrixStats’:\n",
      "\n",
      "    colAlls, colAnyNAs, colAnys, colAvgsPerRowSet, colCollapse,\n",
      "    colCounts, colCummaxs, colCummins, colCumprods, colCumsums,\n",
      "    colDiffs, colIQRDiffs, colIQRs, colLogSumExps, colMadDiffs,\n",
      "    colMads, colMaxs, colMeans2, colMedians, colMins, colOrderStats,\n",
      "    colProds, colQuantiles, colRanges, colRanks, colSdDiffs, colSds,\n",
      "    colSums2, colTabulates, colVarDiffs, colVars, colWeightedMads,\n",
      "    colWeightedMeans, colWeightedMedians, colWeightedSds,\n",
      "    colWeightedVars, rowAlls, rowAnyNAs, rowAnys, rowAvgsPerColSet,\n",
      "    rowCollapse, rowCounts, rowCummaxs, rowCummins, rowCumprods,\n",
      "    rowCumsums, rowDiffs, rowIQRDiffs, rowIQRs, rowLogSumExps,\n",
      "    rowMadDiffs, rowMads, rowMaxs, rowMeans2, rowMedians, rowMins,\n",
      "    rowOrderStats, rowProds, rowQuantiles, rowRanges, rowRanks,\n",
      "    rowSdDiffs, rowSds, rowSums2, rowTabulates, rowVarDiffs, rowVars,\n",
      "    rowWeightedMads, rowWeightedMeans, rowWeightedMedians,\n",
      "    rowWeightedSds, rowWeightedVars\n",
      "\n",
      "\n",
      "Loading required package: GenomicRanges\n",
      "\n",
      "Loading required package: IRanges\n",
      "\n",
      "\n",
      "Attaching package: ‘IRanges’\n",
      "\n",
      "\n",
      "The following objects are masked from ‘package:dplyr’:\n",
      "\n",
      "    collapse, desc, slice\n",
      "\n",
      "\n",
      "The following object is masked from ‘package:purrr’:\n",
      "\n",
      "    reduce\n",
      "\n",
      "\n",
      "Loading required package: GenomeInfoDb\n",
      "\n",
      "Loading required package: Biobase\n",
      "\n",
      "Welcome to Bioconductor\n",
      "\n",
      "    Vignettes contain introductory material; view with\n",
      "    'browseVignettes()'. To cite Bioconductor, see\n",
      "    'citation(\"Biobase\")', and for packages 'citation(\"pkgname\")'.\n",
      "\n",
      "\n",
      "\n",
      "Attaching package: ‘Biobase’\n",
      "\n",
      "\n",
      "The following object is masked from ‘package:MatrixGenerics’:\n",
      "\n",
      "    rowMedians\n",
      "\n",
      "\n",
      "The following objects are masked from ‘package:matrixStats’:\n",
      "\n",
      "    anyMissing, rowMedians\n",
      "\n",
      "\n",
      "\n",
      "Attaching package: ‘SingleCellExperiment’\n",
      "\n",
      "\n",
      "The following object is masked from ‘package:edgeR’:\n",
      "\n",
      "    cpm\n",
      "\n",
      "\n",
      "\n",
      "Attaching package: ‘data.table’\n",
      "\n",
      "\n",
      "The following object is masked from ‘package:SummarizedExperiment’:\n",
      "\n",
      "    shift\n",
      "\n",
      "\n",
      "The following object is masked from ‘package:GenomicRanges’:\n",
      "\n",
      "    shift\n",
      "\n",
      "\n",
      "The following object is masked from ‘package:IRanges’:\n",
      "\n",
      "    shift\n",
      "\n",
      "\n",
      "The following objects are masked from ‘package:S4Vectors’:\n",
      "\n",
      "    first, second\n",
      "\n",
      "\n",
      "The following objects are masked from ‘package:reshape2’:\n",
      "\n",
      "    dcast, melt\n",
      "\n",
      "\n",
      "The following objects are masked from ‘package:dplyr’:\n",
      "\n",
      "    between, first, last\n",
      "\n",
      "\n",
      "The following object is masked from ‘package:purrr’:\n",
      "\n",
      "    transpose\n",
      "\n",
      "\n",
      "Attaching SeuratObject\n",
      "\n",
      "\n",
      "Attaching package: ‘Seurat’\n",
      "\n",
      "\n",
      "The following object is masked from ‘package:SummarizedExperiment’:\n",
      "\n",
      "    Assays\n",
      "\n",
      "\n",
      "\n",
      "Attaching package: ‘patchwork’\n",
      "\n",
      "\n",
      "The following object is masked from ‘package:cowplot’:\n",
      "\n",
      "    align_plots\n",
      "\n",
      "\n",
      "here() starts at /global/u1/k/kserrano/Desktop\n",
      "\n",
      "Loading required package: viridisLite\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# Load libraries\n",
    "library(tidyverse)\n",
    "library(cowplot)\n",
    "library(Matrix.utils)\n",
    "library(edgeR)\n",
    "library(Matrix)\n",
    "library(reshape2)\n",
    "library(S4Vectors)\n",
    "library(SingleCellExperiment)\n",
    "library(pheatmap)\n",
    "library(apeglm)\n",
    "library(png)\n",
    "library(DESeq2)\n",
    "library(RColorBrewer)\n",
    "library(data.table)\n",
    "library(Seurat)\n",
    "library(ggplot2)\n",
    "library(patchwork)\n",
    "library(dplyr)\n",
    "library(here)\n",
    "library(tidyverse)\n",
    "library(viridis)\n",
    "library(RColorBrewer)\n",
    "options (future.globals.maxSize = 4000 * 1024^5)\n",
    "\n",
    "#color scheme\n",
    "farmer_colors = c(\"#FF9999\",\n",
    "\"#FF6666\",\n",
    "\"#FFCC66\",\n",
    "\"#FF6699\",\n",
    "\"#CC99FF\",\n",
    "\"#CC6699\",\n",
    "\"#99CCFF\",\n",
    "\"#ff9933\",\n",
    "\"#FF99FF\",\n",
    "\"#FF99CC\",\n",
    "\"#99CCCC\",\n",
    "\"#66CC99\",\n",
    "\"#CC6600\",\n",
    "\"#FFCCCC\",\n",
    "\"#339966\",\n",
    "\"#CCCCCC\",\n",
    "\"#CC9999\",\n",
    "\"#996666\",\n",
    "\"#FFFF99\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "64135e62-41a9-4f22-9d9a-63a43c6aee5c",
   "metadata": {},
   "outputs": [],
   "source": [
    "#loading in all integrated objects\n",
    "wt.integrated <- readRDS('/global/homes/k/kserrano/Desktop/Integrated_Objects/wt.integrated2.rds')\n",
    "ctrl.integrated <- readRDS('/global/homes/k/kserrano/Desktop/Integrated_Objects/ctrl.integrated1.rds')\n",
    "myc.integrated <- readRDS('/global/homes/k/kserrano/Desktop/Integrated_Objects/myc.integrated.rds')\n",
    "wtnuc <- readRDS(\"/global/cfs/projectdirs/m342/plant-microbe_shared/rObjects/med_all_april5.RDS\")\n",
    "myc.nuc <- readRDS('/global/cfs/projectdirs/m342/plant-microbe_shared/rObjects/med.RDS')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4e560b36-9cbe-4fce-92d9-f520c1f1f9a4",
   "metadata": {},
   "outputs": [],
   "source": [
    "# specify that we will perform downstream analysis on the corrected data note that the\n",
    "# original unmodified data still resides in the 'Spatial' assay\n",
    "DefaultAssay(wt.integrated) <- \"integrated\"\n",
    "DefaultAssay(ctrl.integrated) <- \"integrated\"\n",
    "DefaultAssay(myc.integrated) <- \"integrated\"\n",
    "DefaultAssay(wt.nuc) <- \"integrated\"\n",
    "DefaultAssay(myc.nuc) <- \"integrated\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d298a7cd-6d24-4611-895f-6e4c90ae5b38",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Figure 2 panel B code\n",
    "##First, let's make our list of AM marker genes from Table S2\n",
    "am.markers <-c( 'A17-----------MTR-6g021805' #MtEF1a,\n",
    "'A17-----------CCAMK',\n",
    "'A17-----------MTR-3g109610',\n",
    "'A17-----------CNGC15C',\n",
    "'A17-----------DMI1',\n",
    "'A17-----------ENOD11',\n",
    "'A17-----------MTR-5g026850',\n",
    "'A17-----------LYK3',\n",
    "'A17-----------MTR-7g100110',\n",
    "'A17-----------MTR-3g093270',\n",
    "'A17-----------MTR-3g097560',\n",
    "'A17-----------MTR-1g112940',\n",
    "'A17-----------MTR-8g038210',\n",
    "'A17-----------MTR-0041s0030',\n",
    "'A17-----------MTR-8g095040',\n",
    "'A17-----------MTR-1g471050',\n",
    "'A17-----------MTR-3g065980',\n",
    "'A17-----------MTR-8g068265',\n",
    "'A17-----------MTR-1g078400',\n",
    "'A17-----------MTR-1g017910',\n",
    "'A17-----------MTR-1g109110',\n",
    "'A17-----------MTR-1g090920',\n",
    "'A17-----------MTR-8g006790',\n",
    "'A17-----------MTR-7g116650',\n",
    "'A17-----------MTR-5g031030',\n",
    "'A17-----------MTR-1g028600',\n",
    "'A17-----------MTR-4g104020',\n",
    "'A17-----------MTR-7g027190',\n",
    "'A17-----------MTR-1g040500',\n",
    "'A17-----------MTR-3g118160',\n",
    "'A17-----------MTR-5g030910',\n",
    "'A17-----------MTR-3g089125',\n",
    "'A17-----------MTR-6g027840',\n",
    "'A17-----------MTR-8g012805',\n",
    "'A17-----------MTR-8g012835',\n",
    "'A17-----------MTR-7g068600',\n",
    "'R.irregularis-RIR-0168300',\n",
    "'R.irregularis-RIR-2656400',\n",
    "'R.irregularis-RIR-3143500')\n",
    "               \n",
    "#now let's make the dotplot using the list of AM marker genes  \n",
    "plot00 <- DotPlot(myc.integrated, assay = \"SCT\", features = am.markers, cols = c(\"darkblue\", \"red\"), col.min = -3, col.max = 3, dot.min = 0, dot.scale = 10, idents = NULL, group.by = NULL, split.by = NULL, cluster.idents = TRUE, scale = TRUE, scale.by = \"radius\", scale.min = NA, scale.max = NA) + RotatedAxis()          \n",
    "plot00"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a25436e9-4c10-4d55-b9cd-a443eb806831",
   "metadata": {},
   "outputs": [],
   "source": [
    "##code for Figure 2C\n",
    "#Integrating single nuclei myc and spatial myc datasets\n",
    "anchors <- FindTransferAnchors(reference = myc.nuc, query = myc.integrated, normalization.method = \"SCT\")\n",
    "imputation.assay <- TransferData(anchorset = anchors, refdata = GetAssayData(myc.nuc[['RNA']]),\n",
    "weight.reduction = myc.integrated[[\"pca\"]], dims = 1:30)\n",
    "myc.integrated[[\"imputation\"]] <- imputation.assay\n",
    "DefaultAssay(myc.integrated) <- \"imputation\"\n",
    "\n",
    "#for example, these genes were not detected in the spatial assay, but their location can now be predicted based on concurrently expressed genes in the single cell/single nucleus data\n",
    "p04 <- SpatialPlot(myc.integrated, features = \"A17-----------MTR-0041s0030\", images = \"slice1.2\", alpha =0.5) + scale_fill_viridis_c() #Skl1\n",
    "p05 <- SpatialPlot(myc.integrated, features = \"A17-----------DMI1\", images = \"slice1.2\", alpha =0.5) + scale_fill_viridis_c() #DMI1\n",
    "\n",
    "#lets now save it as svg\n",
    "ggsave(file=\"fig2panelc1.svg\", plot=p04, width=12, height=10)\n",
    "ggsave(file=\"fig2panelc2.svg\", plot=p05, width=12, height=10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0ab0b30f-c7fb-4141-9dd4-cf4d8613132d",
   "metadata": {},
   "outputs": [],
   "source": [
    "#proceeding to integration for prediction of cell types within spatial clusters\n",
    "##simplifying cluster identities into broad categories, pre-annotated object!\n",
    "###beginning of code for Figure 4b!\n",
    "Idents(wt.nuc) <- \"cell_type\"\n",
    "wt.nuc <- RenameIdents(wt.nuc,\"Cortex 1\"=\"Cortex\")\n",
    "wt.nuc<- RenameIdents(object = wt.nuc, `Cortex 2` = \"Cortex\")\n",
    "wt.nuc <- RenameIdents(object = wt.nuc, `Cortex 3` = \"Cortex\")\n",
    "wt.nuc <- RenameIdents(object = wt.nuc, `Cortex 4` = \"Cortex\")\n",
    "wt.nuc <- RenameIdents(object = wt.nuc, `Cortex 5` = \"Cortex\")\n",
    "wt.nuc <- RenameIdents(object = wt.nuc, `Cortex 6` = \"Cortex\")\n",
    "wt.nuc <- RenameIdents(object = wt.nuc, `Cortex 7` = \"Cortex\")\n",
    "wt.nuc <- RenameIdents(object = wt.nuc, `Cortex 8` = \"Cortex\")\n",
    "wt.nuc<- RenameIdents(object = wt.nuc, `LRP 1` = \"LRP\")\n",
    "wt.nuc <- RenameIdents(object = wt.nuc, `LRP 2` = \"LRP\")\n",
    "wt.nuc <- RenameIdents(object = wt.nuc, `Stele 1` = \"Stele\")\n",
    "wt.nuc <- RenameIdents(object = wt.nuc, `Stele 2` = \"Stele\")\n",
    "wt.nuc <- RenameIdents(object = wt.nuc, `Stele 3` = \"Stele\")\n",
    "wt.nuc  <- RenameIdents(object = wt.nuc, `Vascular Pericycle 1` = \"Vascular Pericycle\")\n",
    "wt.nuc <- RenameIdents(object = wt.nuc, `Vascular Pericycle 2` = \"Vascular Pericycle\")\n",
    "wt.nuc$cell_type<- Idents(wt.nuc)\n",
    "#check metadata\n",
    "wtnuc[[]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "09f7e7ab-d0a7-4baf-8585-4107bbfdce3c",
   "metadata": {},
   "outputs": [],
   "source": [
    "#beginning to transfer data from reference (single-nuclei) to query (spatial)\n",
    "wt.query <- wt.integrated\n",
    "wt.anchors <- FindTransferAnchors(reference = wt.nuc, reference.assay = \"integrated\", query = wt.query,\n",
    "    dims = 1:30, reference.reduction = \"pca\")\n",
    "\n",
    "#storing cell type prediction information as metadata into original spatial object\n",
    "predictions <- TransferData(anchorset = wt.anchors, refdata = wt.nuc$cell_type,\n",
    "    dims = 1:30)\n",
    "wt.integrated <- AddMetaData(wt.integrated, metadata = predictions)\n",
    "wt.integrated <- MapQuery(anchorset = wt.anchors, reference = wt.nuc, query = wt.query,\n",
    "    refdata = list(celltype = \"cell_type\"), reference.reduction = \"pca\", reduction.model = \"umap\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b32ec33a-e1bc-4fc0-a673-f60a0326fc1f",
   "metadata": {},
   "outputs": [],
   "source": [
    "#comparing cell type labels now between reference and query\n",
    "p1 <- DimPlot(wt.nuc, reduction = \"umap\", group.by = \"cell_type\", label = TRUE, label.size = 3,\n",
    "    repel = TRUE) + NoLegend() + ggtitle(\"Reference annotations\")\n",
    "p2 <- DimPlot(wt.integrated, reduction = \"ref.umap\", group.by = \"predicted.celltype\", label = TRUE,\n",
    "    label.size = 3, repel = TRUE)  + ggtitle(\"Query transferred labels\")\n",
    "p1 + p2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "5ebbb0b3-c74a-4e8a-a9d2-5e933b1d0871",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "ename": "ERROR",
     "evalue": "Error in FetchData(wt.integrated, vars = c(\"predicted.celltype\", \"seurat_clusters\")): could not find function \"FetchData\"\n",
     "output_type": "error",
     "traceback": [
      "Error in FetchData(wt.integrated, vars = c(\"predicted.celltype\", \"seurat_clusters\")): could not find function \"FetchData\"\nTraceback:\n"
     ]
    }
   ],
   "source": [
    "#now lets make a barplot of predicted cell types within each cluster\n",
    "FetchData(wt.integrated, vars = c(\"predicted.celltype\", \"seurat_clusters\"))\n",
    "ggplot(wt.integrated@meta.data, aes(x = seurat_clusters, fill = predicted.id)) + geom_bar() +\n",
    "    theme_classic()\n",
    "\n",
    "# SVG graphics device\n",
    "svg(\"/global/homes/k/kserrano/Desktop/plot2.svg\")\n",
    "\n",
    "# Code of the plot\n",
    "ggplot(wt@meta.data, aes(x = seurat_clusters, fill = predicted.id)) + geom_bar() +\n",
    "    theme_classic()\n",
    "\n",
    "# Close the graphics device\n",
    "dev.off() \n",
    "\n",
    "###end of code for Figure 4b!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "84235fa7-777f-4ce7-85cd-9abbf5702508",
   "metadata": {},
   "outputs": [],
   "source": [
    "#beginning of differential gene expression analysis for spatial dataset and data for Figure 4a\n",
    "#first done for all cells, ignoring clustering, changing cluster identities to a single cluster called \"all\"\n",
    "wt$ident<- \"all\"\n",
    "head(wt[[]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "23b7b22e-bb6d-46bc-9b2f-6fd5906d6e29",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Bring in Seurat object\n",
    "seurat <- wt\n",
    "\n",
    "# Extract raw counts and metadata to create SingleCellExperiment object\n",
    "counts <- seurat@assays$Spatial@counts \n",
    "\n",
    "metadata <- seurat@meta.data\n",
    "\n",
    "# Set up metadata as desired for aggregation and DE analysis\n",
    "metadata$ident <- factor(seurat$ident)\n",
    "\n",
    "# Create single cell experiment object\n",
    "sce <- SingleCellExperiment(assays = list(counts = counts), \n",
    "                           colData = metadata)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c99abd3b-160b-473c-9eb2-86c2173f5f57",
   "metadata": {},
   "outputs": [],
   "source": [
    "#now on to single cell experiment\n",
    "assays(sce)\n",
    "dim(counts(sce))\n",
    "counts(sce)[1:6, 1:6]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0b48b46d-fee4-4604-a8b4-7623d57d6d84",
   "metadata": {},
   "outputs": [],
   "source": [
    "#exploring cellular metadata\n",
    "dim(colData(sce))\n",
    "head(colData(sce))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b230abdc-c2fd-4dee-88a4-6ac21d13ea6f",
   "metadata": {},
   "outputs": [],
   "source": [
    "#preparing single cell dataset for pseudobulk\n",
    "# Extract unique names of clusters (= levels of cluster_id factor variable)\n",
    "cluster_names <- levels(colData(sce)$ident)\n",
    "cluster_names\n",
    "\n",
    "# Total number of clusters\n",
    "length(cluster_names)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ce7e521e-3cc5-4571-82e0-f967370e98c6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Extract unique names of samples (= levels of sample_id factor variable)\n",
    "sample_names <- levels(colData(sce)$sample)\n",
    "sample_names\n",
    "\n",
    "# Total number of clusters\n",
    "length(sample_names)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "66bd3c7f-5e88-4799-a448-c3960ff47c0a",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Subset metadata to include only the variables you want to aggregate across \n",
    "#(here, we want to aggregate by sample and by cluster)\n",
    "groups <- colData(sce)[, c(\"ident\", \"sample\")]\n",
    "head(groups)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "900e766a-f121-40c5-82c0-5b24ef1ecd45",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Aggregate across cluster-sample groups\n",
    "# transposing row/columns to have cell_ids as row names matching those of groups\n",
    "aggr_counts <- aggregate.Matrix(t(counts(sce)), \n",
    "                                groupings = groups, fun = \"sum\") \n",
    "\n",
    "# Explore output matrix\n",
    "class(aggr_counts)\n",
    "dim(aggr_counts)\n",
    "aggr_counts[1:6, 1:6]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6b500bb2-8c09-4f1b-b5d8-909e6d879cc9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Transpose aggregated matrix to have genes as rows and samples as columns\n",
    "aggr_counts <- t(aggr_counts)\n",
    "aggr_counts[1:6, 1:6]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3889182c-5bf6-4cfa-9405-ceebc5fbb414",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Understanding tstrsplit()\n",
    "\n",
    "## Exploring structure of function output (list)\n",
    "tstrsplit(colnames(aggr_counts), \"_\") %>% str()\n",
    "\n",
    "## Comparing the first 10 elements of our input and output strings\n",
    "head(colnames(aggr_counts), n = 10)\n",
    "head(tstrsplit(colnames(aggr_counts), \"_\")[[1]], n = 10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ebbd9c03-886a-4ee0-b549-0b4368d17e95",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Using which() to look up tstrsplit() output\n",
    "all_idx <- which(tstrsplit(colnames(aggr_counts), \"_\")[[1]] == \"all\")\n",
    "all_idx\n",
    "\n",
    "colnames(aggr_counts)[all_idx]\n",
    "aggr_counts[1:10, all_idx]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "abe95676-ea92-413e-a27a-8cb8f11df29f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# As a reminder, we stored our cell types in a vector called cluster_names\n",
    "cluster_names\n",
    "\n",
    "\n",
    "# Loop over all cell types to extract corresponding counts, and store information in a list\n",
    "\n",
    "## Initiate empty list\n",
    "counts_ls <- list()\n",
    "\n",
    "for (i in 1:length(cluster_names)) {\n",
    "\n",
    "  ## Extract indexes of columns in the global matrix that match a given cluster\n",
    "  column_idx <- which(tstrsplit(colnames(aggr_counts), \"_\")[[1]] == cluster_names[i])\n",
    "  \n",
    "  ## Store corresponding sub-matrix as one element of a list\n",
    "  counts_ls[[i]] <- aggr_counts[, column_idx]\n",
    "  names(counts_ls)[i] <- cluster_names[i]\n",
    "\n",
    "}\n",
    "\n",
    "# Explore the different components of the list\n",
    "str(counts_ls)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1ce271f9-bea7-499b-8b6b-6ab260b589cc",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Reminder: explore structure of metadata\n",
    "head(colData(sce))\n",
    "\n",
    "# Extract sample-level variables\n",
    "metadata <- colData(sce) %>% \n",
    "  as.data.frame() %>% \n",
    "  dplyr::select(sample, treatment)\n",
    "\n",
    "dim(metadata)\n",
    "head(metadata)\n",
    "\n",
    "# Exclude duplicated rows\n",
    "metadata <- metadata[!duplicated(metadata), ]\n",
    "\n",
    "dim(metadata)\n",
    "head(metadata)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ed2dbb3b-1519-41db-a0ac-45c77b3e6690",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Rename rows\n",
    "rownames(metadata) <- metadata$sample\n",
    "head(metadata)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f781f55a-6484-4a63-aeb3-ed721667a5aa",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Number of cells per sample and cluster\n",
    "t <- table(colData(sce)$sample,\n",
    "           colData(sce)$ident)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4a469916-c2cc-4b76-9fcd-632464ccd472",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Creating metadata list\n",
    "\n",
    "## Initiate empty list\n",
    "metadata_ls <- list()\n",
    "\n",
    "for (i in 1:length(counts_ls)) {\n",
    "  \n",
    "    ## Initiate a data frame for cluster i with one row per sample (matching column names in the counts matrix)\n",
    "    df <- data.frame(ident_sample = colnames(counts_ls[[i]]))\n",
    "    \n",
    "    ## Use tstrsplit() to separate cluster (cell type) and sample IDs\n",
    "    df$ident <- tstrsplit(df$ident_sample, \"_\")[[1]]\n",
    "    df$sample  <- tstrsplit(df$ident_sample, \"_\")[[2]]\n",
    "    \n",
    "    \n",
    "    ## Retrieve cell count information for this cluster from global cell count table\n",
    "    idx <- which(colnames(t) == unique(df$ident))\n",
    "    cell_counts <- t[, idx]\n",
    "    \n",
    "    ## Remove samples with zero cell contributing to the cluster\n",
    "    cell_counts <- cell_counts[cell_counts > 0]\n",
    "    \n",
    "    ## Match order of cell_counts and sample_ids\n",
    "    sample_order <- match(df$sample, names(cell_counts))\n",
    "    cell_counts <- cell_counts[sample_order]\n",
    "    \n",
    "    ## Append cell_counts to data frame\n",
    "    df$cell_count <- cell_counts\n",
    "    \n",
    "    \n",
    "    ## Join data frame (capturing metadata specific to cluster) to generic metadata\n",
    "    df <- plyr::join(df, metadata, \n",
    "                     by = intersect(names(df), names(metadata)))\n",
    "    \n",
    "    ## Update rownames of metadata to match colnames of count matrix, as needed later for DE\n",
    "    rownames(df) <- df$ident_sample\n",
    "    \n",
    "    ## Store complete metadata for cluster i in list\n",
    "    metadata_ls[[i]] <- df\n",
    "    names(metadata_ls)[i] <- unique(df$ident)\n",
    "\n",
    "}\n",
    "\n",
    "# Explore the different components of the list\n",
    "str(metadata_ls)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1cb42c41-910d-4d75-8cc5-d98c5765fdec",
   "metadata": {},
   "outputs": [],
   "source": [
    "#now let's actually look at differential expression!!\n",
    "# Select cell type of interest\n",
    "cluster_names\n",
    "\n",
    "# Double-check that both lists have same names\n",
    "all(names(counts_ls) == names(metadata_ls))\n",
    "idx <- which(names(counts_ls) == \"all\")\n",
    "cluster_counts <- counts_ls[[idx]]\n",
    "cluster_metadata <- metadata_ls[[idx]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "37e830b9-f58b-4fd4-a62c-07b725348b60",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Check contents of extracted objects\n",
    "cluster_counts[1:6, 1:6]\n",
    "head(cluster_metadata)\n",
    "\n",
    "# Check matching of matrix columns and metadata rows\n",
    "all(colnames(cluster_counts) == rownames(cluster_metadata))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4670aee9-cfb7-4836-91a2-47d6de01036e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create DESeq2 object        \n",
    "dds <- DESeqDataSetFromMatrix(cluster_counts, \n",
    "                              colData = cluster_metadata, \n",
    "                              design = ~ treatment)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e89f1cc2-d468-4086-8763-5131f8dd9656",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Transform counts for data visualization\n",
    "rld <- rlog(dds, blind=TRUE)\n",
    "\n",
    "# Plot PCA\n",
    "DESeq2::plotPCA(rld, ntop = 500, intgroup = \"treatment\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f9932e10-161e-4ab3-a924-adc9ef8a0efa",
   "metadata": {},
   "outputs": [],
   "source": [
    "DESeq2::plotPCA(rld, ntop = 500, intgroup = \"cell_count\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0195e483-c4eb-4279-a035-a8ff69a44700",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Extract the rlog matrix from the object and compute pairwise correlation values\n",
    "rld_mat <- assay(rld)\n",
    "rld_cor <- cor(rld_mat)\n",
    "\n",
    "# Plot heatmap\n",
    "pheatmap(rld_cor, annotation = cluster_metadata[, c(\"treatment\"), drop=F])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0143b9b1-f66d-47ce-8c78-172c3ade2403",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Run DESeq2 differential expression analysis\n",
    "dds <- DESeq(dds)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e68ab5d0-250a-489c-933f-95ee162d45ef",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot dispersion estimates\n",
    "plotDispEsts(dds)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1a67cd02-f15c-4e8c-addb-307edc1677bb",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Check the coefficients for the comparison\n",
    "resultsNames(dds)\n",
    "\n",
    "# Generate results object\n",
    "res <- results(dds, \n",
    "               name = \"treatment_Myc_vs_Ctrl\",\n",
    "               alpha = 0.05)\n",
    "\n",
    "# Shrink the log2 fold changes to be more appropriate using the apeglm method - should cite [paper]() when using this method\n",
    "res <- lfcShrink(dds, \n",
    "                 coef = \"treatment_Myc_vs_Ctrl\",\n",
    "                 res=res,\n",
    "                 type = \"apeglm\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6473e337-b20d-4cd4-a5d7-da236271c5f2",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Turn the DESeq2 results object into a tibble for use with tidyverse functions\n",
    "res_tbl <- res %>%\n",
    "  data.frame() %>%\n",
    "  rownames_to_column(var = \"gene\") %>%\n",
    "  as_tibble() %>%\n",
    "  arrange(padj)\n",
    "\n",
    "# Check results output\n",
    "head(res_tbl)\n",
    "\n",
    "# Write all results to file\n",
    "write.csv(res_tbl, \"/global/homes/k/kserrano/Desktop/resultsall.csv\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6f5a15f1-ba7b-42cb-9fc0-2493556244a9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set thresholds\n",
    "padj_cutoff <- 0.05\n",
    "\n",
    "# Subset the significant results\n",
    "sig_res <- dplyr::filter(res_tbl, padj < padj_cutoff) %>%\n",
    "  dplyr::arrange(padj)\n",
    "\n",
    "# Check significant genes output\n",
    "sig_res\n",
    "\n",
    "# Write significant results to file\n",
    "write.csv(sig_res, \"/global/homes/k/kserrano/Desktop/sig.resultsall.csv\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8051b3fc-2469-4811-8978-015ec9efedf0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set thresholds\n",
    "log2fc_cutoff <- 1.0\n",
    "\n",
    "# Count significantly up/down genes above threshold\n",
    "n_sig_up <- dplyr::filter(sig_res, log2FoldChange >= log2fc_cutoff) %>% \n",
    "  nrow()\n",
    "n_sig_dn <- dplyr::filter(sig_res, log2FoldChange <= -log2fc_cutoff) %>% \n",
    "  nrow()\n",
    "\n",
    "n_sig_up\n",
    "n_sig_down"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ae7dbac4-d693-4cc7-97d1-e3ee9f1f9c02",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Scatterplot\n",
    "\n",
    "## Extract normalized counts from dds object\n",
    "normalized_counts <- counts(dds, normalized = TRUE)\n",
    "\n",
    "## Extract top 20 DEG from resLFC (make sure to order by padj)\n",
    "top20_sig_genes <- sig_res %>%\n",
    "  dplyr::arrange(padj) %>%\n",
    "  dplyr::pull(gene) %>%\n",
    "  head(n = 20)\n",
    "\n",
    "## Extract matching normalized count values from matrix\n",
    "top20_sig_counts <- normalized_counts[rownames(normalized_counts) %in% top20_sig_genes, ]\n",
    "top20_sig_counts\n",
    "\n",
    "## Convert wide matrix to long data frame for ggplot2\n",
    "top20_sig_df <- data.frame(top20_sig_counts)\n",
    "top20_sig_df$gene <- rownames(top20_sig_counts)\n",
    "\n",
    "top20_sig_df <- melt(setDT(top20_sig_df), \n",
    "                     id.vars = c(\"gene\"),\n",
    "                     variable.name = \"ident_sample\") %>% \n",
    "  data.frame()\n",
    "top20_sig_df\n",
    "\n",
    "## Replace \".\" by \" \" in cluster_sample_id variable (melt() introduced the \".\")\n",
    "top20_sig_df$ident_sample <- gsub(\"\\\\X\", \"\", top20_sig_df$ident_sample)\n",
    "top20_sig_df\n",
    "\n",
    "##make a new dataframe for dds metadata \n",
    "df2 <- as.data.frame(colData(dds))\n",
    "df2\n",
    "\n",
    "\n",
    "## Join counts data frame with metadata\n",
    "top20_sig_df <- plyr::join(top20_sig_df, df2,\n",
    "                           by = \"ident_sample\")\n",
    "top20_sig_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ba76c684-d8e7-483d-ba29-1fdccf016424",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Generate plot\n",
    "plot1 <- ggplot(top20_sig_df, aes(y = value, x = treatment, col = treatment)) +\n",
    "  geom_jitter(height = 0, width = 0.15) +\n",
    "  scale_y_continuous(trans = 'log10') +\n",
    "  ylab(\"log10 of normalized expression level\") +\n",
    "  xlab(\"condition\") +\n",
    "  ggtitle(\"Top 20 Significant DE Genes\") +\n",
    "  theme(plot.title = element_text(hjust = 0.5)) +\n",
    "  facet_wrap(~ gene)\n",
    "\n",
    "svg(\"top20spatial.svg\")\n",
    "\n",
    "# Code of the plot\n",
    "plot1\n",
    "\n",
    "# Close the graphics device\n",
    "dev.off() \n",
    "\n",
    "top20_sig_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a6cd8eb5-a83e-45c3-a673-4ad111963719",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Heatmap\n",
    "\n",
    "## Extract normalized counts for significant genes only\n",
    "sig_counts <- normalized_counts[rownames(normalized_counts) %in% sig_res$gene, ]\n",
    "\n",
    "## Set a color-blind friendly palette\n",
    "heat_colors <- rev(brewer.pal(11, \"PuOr\"))\n",
    "\n",
    "## Run pheatmap using the metadata data frame for the annotation\n",
    "pheatmap(sig_counts, \n",
    "         color = heat_colors, \n",
    "         cluster_rows = TRUE, \n",
    "         show_rownames = FALSE,\n",
    "         annotation = cluster_metadata[, c(\"treatment\", \"ident\")], \n",
    "         border_color = NA, \n",
    "         fontsize = 10, \n",
    "         scale = \"row\", \n",
    "         fontsize_row = 10, \n",
    "         height = 20) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "aca0ebd9-4c08-41d6-8f92-1c5ad2d3a421",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Volcano plot\n",
    "res_table_thres <- res_tbl[!is.na(res_tbl$padj), ] %>% \n",
    "  mutate(threshold = padj < padj_cutoff & abs(log2FoldChange) >= log2fc_cutoff)\n",
    "min(log10(res_table_thres$padj))\n",
    "\n",
    "## Generate plot\n",
    "ggplot(res_table_thres) +\n",
    "  geom_point(aes(x = log2FoldChange, y = -log10(padj), colour = threshold)) +\n",
    "  ggtitle(\"Volcano plot of myc cluster 0 cells relative to ctrl\") +\n",
    "  xlab(\"log2 fold change\") +\n",
    "  xlim(-4.5, 12) +\n",
    "  ylab(\"-log10 adjusted p-value\") +\n",
    "  scale_y_continuous(limits = c(0, 100)) +\n",
    "  scale_color_manual(values = c(\"grey60\", \"red3\")) +\n",
    "  theme(legend.position = \"none\",\n",
    "        plot.title = element_text(size = rel(1.3), hjust = 0.5),\n",
    "        axis.title = element_text(size = rel(1.15)))     "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "04ecb255-e069-43e1-9418-3204c038990b",
   "metadata": {},
   "outputs": [],
   "source": [
    "#end of differential gene expression analysis for all voxels at once\n",
    "#moving onto doing this analysis for each individual cluster\n",
    "#code for figure 4c counts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f27c1e16-9f40-41c7-b6d6-63c39e14be08",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Bring in Seurat object\n",
    "seurat <- wt\n",
    "\n",
    "# Extract raw counts and metadata to create SingleCellExperiment object\n",
    "counts <- seurat@assays$Spatial@counts \n",
    "\n",
    "metadata <- seurat@meta.data\n",
    "\n",
    "# Set up metadata as desired for aggregation and DE analysis\n",
    "#instead of using ident, we now use the seurat_clusters metadata column\n",
    "metadata$seurat_clusters <- factor(seurat$seurat_clusters)\n",
    "\n",
    "# Create single cell experiment object\n",
    "sce <- SingleCellExperiment(assays = list(counts = counts), \n",
    "                           colData = metadata)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fe2756ab-3947-4e2e-b83b-c962334d9726",
   "metadata": {},
   "outputs": [],
   "source": [
    "#now on to single cell experiment\n",
    "assays(sce)\n",
    "dim(counts(sce))\n",
    "counts(sce)[1:6, 1:6]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "411db1df-4071-4707-ae38-8361405539f5",
   "metadata": {},
   "outputs": [],
   "source": [
    "#exploring cellular metadata\n",
    "dim(colData(sce))\n",
    "head(colData(sce))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6a19f18d-abde-482a-b4ef-da9bb9ee8767",
   "metadata": {},
   "outputs": [],
   "source": [
    "#preparing single cell dataset for pseudobulk\n",
    "# Extract unique names of clusters (= levels of cluster_id factor variable)\n",
    "cluster_names <- levels(colData(sce)$seurat_clusters)\n",
    "cluster_names\n",
    "\n",
    "# Total number of clusters\n",
    "length(cluster_names)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "52bafb1a-3e81-4d0f-b3b8-77d4f8e54dff",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Extract unique names of samples (= levels of sample_id factor variable)\n",
    "sample_names <- levels(colData(sce)$sample)\n",
    "sample_names\n",
    "\n",
    "# Total number of clusters\n",
    "length(sample_names)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f9296573-af9b-4516-870b-f406c152d420",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Subset metadata to include only the variables you want to aggregate across \n",
    "#(here, we want to aggregate by sample and by cluster)\n",
    "groups <- colData(sce)[, c(\"seurat_clusters\", \"sample\")]\n",
    "head(groups)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c2b41d73-b94b-46bc-9840-aabd442d6ae3",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Aggregate across cluster-sample groups\n",
    "# transposing row/columns to have cell_ids as row names matching those of groups\n",
    "aggr_counts <- aggregate.Matrix(t(counts(sce)), \n",
    "                                groupings = groups, fun = \"sum\") \n",
    "\n",
    "# Explore output matrix\n",
    "class(aggr_counts)\n",
    "dim(aggr_counts)\n",
    "aggr_counts[1:6, 1:6]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "85eedec4-4694-4483-874a-d1412b33e981",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Transpose aggregated matrix to have genes as rows and samples as columns\n",
    "aggr_counts <- t(aggr_counts)\n",
    "aggr_counts[1:6, 1:6]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e031c5d8-bc58-441b-abaf-92a1fd278596",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Understanding tstrsplit()\n",
    "\n",
    "## Exploring structure of function output (list)\n",
    "tstrsplit(colnames(aggr_counts), \"_\") %>% str()\n",
    "\n",
    "## Comparing the first 10 elements of our input and output strings\n",
    "head(colnames(aggr_counts), n = 10)\n",
    "head(tstrsplit(colnames(aggr_counts), \"_\")[[1]], n = 10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d12dbe43-d655-400d-a90b-128136eb6709",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Using which() to look up tstrsplit() output\n",
    "#starting with cluster number 0 here\n",
    "zero_idx <- which(tstrsplit(colnames(aggr_counts), \"_\")[[1]] == \"0\")\n",
    "zero_idx\n",
    "\n",
    "colnames(aggr_counts)[zero_idx]\n",
    "aggr_counts[1:10, zero_idx]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "920c6cb0-879c-47fe-94e9-9fae43dda947",
   "metadata": {},
   "outputs": [],
   "source": [
    "# As a reminder, we stored our cell types in a vector called cluster_names\n",
    "cluster_names\n",
    "\n",
    "\n",
    "# Loop over all cell types to extract corresponding counts, and store information in a list\n",
    "\n",
    "## Initiate empty list\n",
    "counts_ls <- list()\n",
    "\n",
    "for (i in 1:length(cluster_names)) {\n",
    "\n",
    "  ## Extract indexes of columns in the global matrix that match a given cluster\n",
    "  column_idx <- which(tstrsplit(colnames(aggr_counts), \"_\")[[1]] == cluster_names[i])\n",
    "  \n",
    "  ## Store corresponding sub-matrix as one element of a list\n",
    "  counts_ls[[i]] <- aggr_counts[, column_idx]\n",
    "  names(counts_ls)[i] <- cluster_names[i]\n",
    "\n",
    "}\n",
    "\n",
    "# Explore the different components of the list\n",
    "str(counts_ls)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1604945a-2611-400e-80af-b260cf709e14",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Reminder: explore structure of metadata\n",
    "head(colData(sce))\n",
    "\n",
    "# Extract sample-level variables\n",
    "metadata <- colData(sce) %>% \n",
    "  as.data.frame() %>% \n",
    "  dplyr::select(sample, treatment)\n",
    "\n",
    "dim(metadata)\n",
    "head(metadata)\n",
    "\n",
    "# Exclude duplicated rows\n",
    "metadata <- metadata[!duplicated(metadata), ]\n",
    "\n",
    "dim(metadata)\n",
    "head(metadata)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5c172487-551f-484f-ae62-b9143d538ea7",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Rename rows\n",
    "rownames(metadata) <- metadata$sample\n",
    "head(metadata)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6f29c907-ec54-4389-b558-cd74227d8d80",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Number of cells per sample and cluster\n",
    "t <- table(colData(sce)$sample,\n",
    "           colData(sce)$seurat_clusters)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ffbbf71c-5fdb-46d7-832b-5fc87d62c2ee",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Creating metadata list\n",
    "\n",
    "## Initiate empty list\n",
    "metadata_ls <- list()\n",
    "\n",
    "for (i in 1:length(counts_ls)) {\n",
    "  \n",
    "    ## Initiate a data frame for cluster i with one row per sample (matching column names in the counts matrix)\n",
    "    df <- data.frame(seurat_clusters_sample = colnames(counts_ls[[i]]))\n",
    "    \n",
    "    ## Use tstrsplit() to separate cluster (cell type) and sample IDs\n",
    "    df$seurat_clusters <- tstrsplit(df$seurat_clusters_sample, \"_\")[[1]]\n",
    "    df$sample  <- tstrsplit(df$seurat_clusters_sample, \"_\")[[2]]\n",
    "    \n",
    "    \n",
    "    ## Retrieve cell count information for this cluster from global cell count table\n",
    "    idx <- which(colnames(t) == unique(df$seurat_clusters))\n",
    "    cell_counts <- t[, idx]\n",
    "    \n",
    "    ## Remove samples with zero cell contributing to the cluster\n",
    "    cell_counts <- cell_counts[cell_counts > 0]\n",
    "    \n",
    "    ## Match order of cell_counts and sample_ids\n",
    "    sample_order <- match(df$sample, names(cell_counts))\n",
    "    cell_counts <- cell_counts[sample_order]\n",
    "    \n",
    "    ## Append cell_counts to data frame\n",
    "    df$cell_count <- cell_counts\n",
    "    \n",
    "    \n",
    "    ## Join data frame (capturing metadata specific to cluster) to generic metadata\n",
    "    df <- plyr::join(df, metadata, \n",
    "                     by = intersect(names(df), names(metadata)))\n",
    "    \n",
    "    ## Update rownames of metadata to match colnames of count matrix, as needed later for DE\n",
    "    rownames(df) <- df$seurat_clusters_sample\n",
    "    \n",
    "    ## Store complete metadata for cluster i in list\n",
    "    metadata_ls[[i]] <- df\n",
    "    names(metadata_ls)[i] <- unique(df$seurat_clusters)\n",
    "\n",
    "}\n",
    "\n",
    "# Explore the different components of the list\n",
    "str(metadata_ls)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3f9244cb-bc55-41c9-b334-e59470485e8d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Select cell type of interest\n",
    "cluster_names\n",
    "\n",
    "# Double-check that both lists have same names\n",
    "all(names(counts_ls) == names(metadata_ls))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e3ed92a3-14e2-411b-9db2-0d4c848ceb15",
   "metadata": {},
   "outputs": [],
   "source": [
    "idx <- which(names(counts_ls) == \"0\")\n",
    "cluster_counts <- counts_ls[[idx]]\n",
    "cluster_metadata <- metadata_ls[[idx]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e3c00b62-239c-498d-845e-e04f74357b94",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Check contents of extracted objects\n",
    "cluster_counts[1:6, 1:6]\n",
    "head(cluster_metadata)\n",
    "\n",
    "# Check matching of matrix columns and metadata rows\n",
    "all(colnames(cluster_counts) == rownames(cluster_metadata))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "da22f5d4-58f7-4c0c-8899-f0fb0a8d3194",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create DESeq2 object        \n",
    "dds0 <- DESeqDataSetFromMatrix(cluster_counts, \n",
    "                              colData = cluster_metadata, \n",
    "                              design = ~ treatment)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ceaba76d-17e9-49e1-b863-b33c128433c6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Transform counts for data visualization\n",
    "rld <- rlog(dds0, blind=TRUE)\n",
    "\n",
    "# Plot PCA\n",
    "DESeq2::plotPCA(rld, ntop = 500, intgroup = \"treatment\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8bb2c56c-45bf-4a17-af30-bdb95b56397a",
   "metadata": {},
   "outputs": [],
   "source": [
    "DESeq2::plotPCA(rld, ntop = 500, intgroup = \"cell_count\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0515d943-16e0-4098-bff8-ff83dad0f49e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Extract the rlog matrix from the object and compute pairwise correlation values\n",
    "rld_mat <- assay(rld)\n",
    "rld_cor <- cor(rld_mat)\n",
    "\n",
    "# Plot heatmap\n",
    "pheatmap(rld_cor, annotation = cluster_metadata[, c(\"treatment\"), drop=F])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5f324e5a-f8b4-4bf9-a150-76c2d7d26920",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Run DESeq2 differential expression analysis\n",
    "dds <- DESeq(dds0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "db495043-a21b-40c9-a8b5-1a5dbe43c708",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot dispersion estimates\n",
    "plotDispEsts(dds)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d503d79d-6e37-4fe8-93e3-8ef6e9cb34b0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Check the coefficients for the comparison\n",
    "resultsNames(dds)\n",
    "\n",
    "# Generate results object\n",
    "res0 <- results(dds, \n",
    "               name = \"treatment_Myc_vs_Ctrl\",\n",
    "               alpha = 0.05)\n",
    "\n",
    "# Shrink the log2 fold changes to be more appropriate using the apeglm method - should cite [paper]() when using this method\n",
    "res0 <- lfcShrink(dds, \n",
    "                 coef = \"treatment_Myc_vs_Ctrl\",\n",
    "                 res=res0,\n",
    "                 type = \"apeglm\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fae7680c-ffd9-4e04-b25a-50bbb6308e37",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Turn the DESeq2 results object into a tibble for use with tidyverse functions\n",
    "res_tbl <- res0 %>%\n",
    "  data.frame() %>%\n",
    "  rownames_to_column(var = \"gene\") %>%\n",
    "  as_tibble() %>%\n",
    "  arrange(padj)\n",
    "\n",
    "# Check results output\n",
    "head(res_tbl)\n",
    "\n",
    "# Write all results to file\n",
    "write.csv(res_tbl, \"/global/homes/k/kserrano/Desktop/results0.csv\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2ae3670f-6e6b-4d72-8ae3-157e74fe0150",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set thresholds\n",
    "padj_cutoff <- 0.05\n",
    "\n",
    "# Subset the significant results\n",
    "sig_res0 <- dplyr::filter(res_tbl, padj < padj_cutoff) %>%\n",
    "  dplyr::arrange(padj)\n",
    "\n",
    "# Check significant genes output\n",
    "sig_res0\n",
    "\n",
    "# Write significant results to file\n",
    "write.csv(sig_res0, \"/global/homes/k/kserrano/Desktop/sig.results0.csv\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "56e7da6b-1d89-4d32-a6a3-908cdcb4c93f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set thresholds\n",
    "log2fc_cutoff <- 1.0\n",
    "\n",
    "# Count significantly up/down genes above threshold\n",
    "n_sig_up0 <- dplyr::filter(sig_res0, log2FoldChange >= log2fc_cutoff) %>% \n",
    "  nrow()\n",
    "n_sig_dn0 <- dplyr::filter(sig_res0, log2FoldChange <= -log2fc_cutoff) %>% \n",
    "  nrow()\n",
    "\n",
    "n_sig_up0\n",
    "n_sig_dn0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4604774d-1261-4b34-93c4-674ff68fb74a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Scatterplot\n",
    "\n",
    "## Extract normalized counts from dds object\n",
    "normalized_counts <- counts(dds, normalized = TRUE)\n",
    "\n",
    "## Extract top 20 DEG from resLFC (make sure to order by padj)\n",
    "top20_sig_genes_0 <- sig_res0 %>%\n",
    "  dplyr::arrange(padj) %>%\n",
    "  dplyr::pull(gene) %>%\n",
    "  head(n = 20)\n",
    "\n",
    "## Extract matching normalized count values from matrix\n",
    "top20_sig_counts_0 <- normalized_counts[rownames(normalized_counts) %in% top20_sig_genes_0, ]\n",
    "top20_sig_counts_0\n",
    "\n",
    "## Convert wide matrix to long data frame for ggplot2\n",
    "top20_sig_df_0 <- data.frame(top20_sig_counts_0)\n",
    "top20_sig_df_0$gene <- rownames(top20_sig_counts_0)\n",
    "\n",
    "top20_sig_df_0 <- melt(setDT(top20_sig_df_0), \n",
    "                     id.vars = c(\"gene\"),\n",
    "                     variable.name = \"seurat_clusters_sample\") %>% \n",
    "  data.frame()\n",
    "top20_sig_df_0\n",
    "\n",
    "## Replace \".\" by \" \" in cluster_sample_id variable (melt() introduced the \".\")\n",
    "top20_sig_df_0$seurat_clusters_sample <- gsub(\"\\\\X\", \"\", top20_sig_df_0$seurat_clusters_sample)\n",
    "top20_sig_df_0\n",
    "\n",
    "##make a new dataframe for dds metadata \n",
    "df2 <- as.data.frame(colData(dds))\n",
    "df2\n",
    "\n",
    "\n",
    "## Join counts data frame with metadata\n",
    "top20_sig_df <- plyr::join(top20_sig_df_0, df2,\n",
    "                           by = \"seurat_clusters_sample\")\n",
    "top20_sig_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "68777568-25dd-4cbc-aa7e-3cdbb2d8fdef",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Generate plot\n",
    "plot1 <- ggplot(top20_sig_df_0, aes(y = value, x = \"treatment\", col = \"treatment\")) +\n",
    "  geom_jitter(height = 0, width = 0.15) +\n",
    "  scale_y_continuous(trans = 'log10') +\n",
    "  ylab(\"log10 of normalized expression level\") +\n",
    "  xlab(\"condition\") +\n",
    "  ggtitle(\"Top 20 Significant DE Genes\") +\n",
    "  theme(plot.title = element_text(hjust = 0.5)) +\n",
    "  facet_wrap(~ gene)\n",
    "\n",
    "svg(\"top20spatial_0.svg\")\n",
    "\n",
    "# Code of the plot\n",
    "plot1\n",
    "\n",
    "# Close the graphics device\n",
    "dev.off() \n",
    "\n",
    "top20_sig_df_0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f7382bae-7712-48c0-91a4-d7cb7c094dcc",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Heatmap\n",
    "\n",
    "## Extract normalized counts for significant genes only\n",
    "sig_counts <- normalized_counts[rownames(normalized_counts) %in% sig_res0$gene, ]\n",
    "\n",
    "## Set a color-blind friendly palette\n",
    "heat_colors <- rev(brewer.pal(11, \"PuOr\"))\n",
    "\n",
    "## Run pheatmap using the metadata data frame for the annotation\n",
    "pheatmap(sig_counts, \n",
    "         color = heat_colors, \n",
    "         cluster_rows = TRUE, \n",
    "         show_rownames = FALSE,\n",
    "         annotation = cluster_metadata[, c(\"treatment\", \"seurat_clusters\")], \n",
    "         border_color = NA, \n",
    "         fontsize = 10, \n",
    "         scale = \"row\", \n",
    "         fontsize_row = 10, \n",
    "         height = 20) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8157d430-3f6e-4ea7-8d36-17f0184566ce",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Volcano plot\n",
    "res_table_thres0 <- res_tbl[!is.na(res_tbl$padj), ] %>% \n",
    "  mutate(threshold = padj < padj_cutoff & abs(log2FoldChange) >= log2fc_cutoff)\n",
    "min(log10(res_table_thres0$padj))\n",
    "\n",
    "## Generate plot\n",
    "ggplot(res_table_thres0) +\n",
    "  geom_point(aes(x = log2FoldChange, y = -log10(padj), colour = threshold)) +\n",
    "  ggtitle(\"Volcano plot of myc cluster 0 cells relative to ctrl\") +\n",
    "  xlab(\"log2 fold change\") +\n",
    "  xlim(-4.5, 12) +\n",
    "  ylab(\"-log10 adjusted p-value\") +\n",
    "  scale_y_continuous(limits = c(0, 100)) +\n",
    "  scale_color_manual(values = c(\"grey60\", \"red3\")) +\n",
    "  theme(legend.position = \"none\",\n",
    "        plot.title = element_text(size = rel(1.3), hjust = 0.5),\n",
    "        axis.title = element_text(size = rel(1.15)))     "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "be73cd5a-240c-46a1-bb61-d609f1ee2832",
   "metadata": {},
   "outputs": [],
   "source": [
    "#trying to do all clusters at once, to avoid doing them individually\n",
    "get_dds_LRTresults <- function(clustx){\n",
    "  \n",
    "  print(clustx) # useful for debugging\n",
    "  \n",
    "  # Extract counts matrix and metadata for cluster x\n",
    "  idx <- which(names(counts_ls) == clustx)\n",
    "  cluster_counts <- counts_ls[[idx]]\n",
    "  cluster_metadata <- metadata_ls[[idx]]\n",
    "  \n",
    "  # Print error message if sample names do not match\n",
    "  if ( all(colnames(cluster_counts) != rownames(cluster_metadata)) ) {\n",
    "    print(\"ERROR: sample names in counts matrix columns and metadata rows do not match!\")\n",
    "  }\n",
    "  \n",
    "  # Run DESeq2\n",
    "  dds <- DESeqDataSetFromMatrix(cluster_counts, \n",
    "                                colData = cluster_metadata, \n",
    "                                design = ~ treatment)\n",
    "  dds_lrt <- DESeq(dds, test = \"LRT\", reduced = ~ 1)\n",
    "  \n",
    "  # Extract results\n",
    "  res_LRT <- results(dds_lrt)\n",
    "  \n",
    "  # Create a tibble for LRT results\n",
    "  res_LRT_tb <- res_LRT %>%\n",
    "    data.frame() %>%\n",
    "    rownames_to_column(var = \"gene\") %>% \n",
    "    as_tibble()\n",
    "  \n",
    "  # Save all results\n",
    "  if (!dir.exists(\"results\")) { dir.create(\"results\") }\n",
    "  write.csv(res_LRT_tb,\n",
    "            paste0(\"results/\", clustx, \"_LRT_all_genes.csv\"),\n",
    "            quote = FALSE, \n",
    "            row.names = FALSE)\n",
    "  \n",
    "  # Subset to return genes with padj < 0.05\n",
    "  sigLRT_genes <- res_LRT_tb %>% \n",
    "    filter(padj < 0.05)\n",
    "  \n",
    "  # Save significant results\n",
    "  write.csv(sigLRT_genes,\n",
    "            paste0(\"results/\", clustx, \"_LRT_signif_genes.csv\"),\n",
    "            quote = FALSE, \n",
    "            row.names = FALSE)\n",
    "  \n",
    "  # Transform counts for data visualization\n",
    "  rld <- rlog(dds_lrt, blind = TRUE)\n",
    "  \n",
    "  # Extract the rlog matrix from the object and compute pairwise correlation values\n",
    "  rld_mat <- assay(rld)\n",
    "  rld_cor <- cor(rld_mat)\n",
    "  \n",
    "  # Obtain rlog values for those significant genes\n",
    "  cluster_rlog <- rld_mat[sigLRT_genes$gene, ]\n",
    "  cluster_meta_sig <- cluster_metadata[which(rownames(cluster_metadata) %in% colnames(cluster_rlog)), ]\n",
    "  \n",
    "  save(dds_lrt, res_LRT, sigLRT_genes, \n",
    "       file = paste0(\"results/\", clustx, \"_all_LRTresults.Rdata\"))\n",
    "  \n",
    "}\n",
    "\n",
    "map(cluster_names, get_dds_LRTresults)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8451505c-9070-46ae-857c-e5180fe0ffd8",
   "metadata": {},
   "outputs": [],
   "source": [
    "#end of differential gene expression!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b7420ff2-000f-4fab-a1fd-19f7f1e5118e",
   "metadata": {},
   "outputs": [],
   "source": [
    "#looking at all rhizophagus gene expression\n",
    "#loading in mycorrhizal spatial objects individually\n",
    "B1 <- readRDS(\"/global/homes/k/kserrano/Desktop/Spatial_Objects/B1.rds\")\n",
    "B2 <- readRDS(\"/global/homes/k/kserrano/Desktop/Spatial_Objects/B2.rds\")\n",
    "B3 <- readRDS(\"/global/homes/k/kserrano/Desktop/Spatial_Objects/B3.rds\")\n",
    "A2 <- readRDS(\"/global/homes/k/kserrano/Desktop/Spatial_Objects/A2.rds\")\n",
    "A3 <- readRDS(\"/global/homes/k/kserrano/Desktop/Spatial_Objects/A3.rds\")\n",
    "C2 <- readRDS(\"/global/homes/k/kserrano/Desktop/Spatial_Objects/C2.rds\")\n",
    "C3 <- readRDS(\"/global/homes/k/kserrano/Desktop/Spatial_Objects/C3.rds\")\n",
    "D2 <- readRDS(\"/global/homes/k/kserrano/Desktop/Spatial_Objects/D2.rds\")\n",
    "D3 <- readRDS(\"/global/homes/k/kserrano/Desktop/Spatial_Objects/D3.rds\")\n",
    "\n",
    "#merging into a single object\n",
    "myc.merge <- merge(B1, y = c(B2, B3, A2, A3, C2, C3, D2, D3), add.cell.ids = c(\"B1\", \"B2\", \"B3\", \"A2\", \"A3\", \"C2\", \"C3\", \"D2\", \"D3\"))\n",
    "myc.merge\n",
    "\n",
    "#getting a list of all genes expressed in this merged object and their expression values\n",
    "#figure 5 panel a - getting list to compare to previous datasets\n",
    "apply(myc.merge@assays$Spatial@data,1,mean) -> gene.expression.myc\n",
    "sort(gene.expression.myc, decreasing = TRUE) -> gene.expression.myc\n",
    "as.data.frame(gene.expression.myc) -> mycexp\n",
    "write.csv(mycexp, \"/global/homes/k/kserrano/Desktop/mycexp.csv\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0b4bce46-e24c-4150-81f9-f40f4b0983e4",
   "metadata": {},
   "outputs": [],
   "source": [
    "#figure 5 panel c\n",
    "#plotting dimplot for a particular gene\n",
    "sobj <- myc.integrated\n",
    "#ggplot function\n",
    "get_reduction <- function(sobj, reduction_name) {\n",
    "  rd <- Embeddings(sobj, reduction = reduction_name) %>% \n",
    "    as_tibble(rownames = \"Cell\")\n",
    "}\n",
    "\n",
    "get_cell_data <- function(sobj, features = NULL, melt = TRUE) {\n",
    "  md <- as_tibble(sobj@meta.data, rownames = \"Cell\")\n",
    "  reducs <- names(sobj@reductions)\n",
    "  rd <- map(reducs, get_reduction, sobj = sobj) %>% reduce(left_join)\n",
    "  combined_data <- left_join(md, rd)\n",
    "  if(!is.null(features)) {\n",
    "    exp_data <- GetAssayData(sobj[features,], slot = \"counts\", assay = \"Spatial\") %>% \n",
    "      Matrix::t() %>% \n",
    "      as_tibble(rownames = \"Cell\")\n",
    "    if(melt) {\n",
    "      exp_data <- gather(exp_data, \"Locus\", \"Counts\", -Cell)\n",
    "    }\n",
    "    combined_data <- left_join(combined_data, exp_data)\n",
    "  }\n",
    "  return(combined_data)\n",
    "}\n",
    "\n",
    "## Plot by expression\n",
    "features <- c(\"R.irregularis-RIR-2656400\") #RiFTR1\n",
    "cell_data <- get_cell_data(sobj, features = features)\n",
    "ggplot(arrange(cell_data, Counts), aes(x = UMAP_1, y = UMAP_2)) +\n",
    "  geom_point(aes(color = log1p(Counts)),size=0.5) +\n",
    "  scale_color_viridis_c(name = \"Log(Counts + 1)\", option = \"plasma\") +\n",
    "  labs(x = \"UMAP 1\", y = \"UMAP 2\") +\n",
    "  facet_wrap(\"Locus\") +\n",
    "  theme_bw() +\n",
    "  theme(panel.grid = element_blank(),\n",
    "        axis.line = element_blank(),\n",
    "        axis.ticks = element_blank(),\n",
    "        axis.text = element_blank())\n",
    "\n",
    "ggsave(\"RiFTR1.pdf\", plot=last_plot(), device=\"pdf\", width = 5.5, height = 4, units=\"in\", dpi = 300)\n",
    "\n",
    "## Plot by expression\n",
    "features <- c(\"R.irregularis-RIR-0168300\") #RiEF1a\n",
    "cell_data <- get_cell_data(sobj, features = features)\n",
    "ggplot(arrange(cell_data, Counts), aes(x = UMAP_1, y = UMAP_2)) +\n",
    "  geom_point(aes(color = log1p(Counts)),size=0.5) +\n",
    "  scale_color_viridis_c(name = \"Log(Counts + 1)\", option = \"plasma\") +\n",
    "  labs(x = \"UMAP 1\", y = \"UMAP 2\") +\n",
    "  facet_wrap(\"Locus\") +\n",
    "  theme_bw() +\n",
    "  theme(panel.grid = element_blank(),\n",
    "        axis.line = element_blank(),\n",
    "        axis.ticks = element_blank(),\n",
    "        axis.text = element_blank())\n",
    "\n",
    "ggsave(\"RiEF1a.pdf\", plot=last_plot(), device=\"pdf\", width = 5.5, height = 4, units=\"in\", dpi = 300)\n",
    "\n",
    "## Plot by expression\n",
    "features <- c(\"R.irregularis-RIR-1036900\") #Rihypothetical\n",
    "cell_data <- get_cell_data(sobj, features = features)\n",
    "ggplot(arrange(cell_data, Counts), aes(x = UMAP_1, y = UMAP_2)) +\n",
    "  geom_point(aes(color = log1p(Counts)),size=0.5) +\n",
    "  scale_color_viridis_c(name = \"Log(Counts + 1)\", option = \"plasma\") +\n",
    "  labs(x = \"UMAP 1\", y = \"UMAP 2\") +\n",
    "  facet_wrap(\"Locus\") +\n",
    "  theme_bw() +\n",
    "  theme(panel.grid = element_blank(),\n",
    "        axis.line = element_blank(),\n",
    "        axis.ticks = element_blank(),\n",
    "        axis.text = element_blank())\n",
    "\n",
    "ggsave(\"Rihyp.pdf\", plot=last_plot(), device=\"pdf\", width = 5.5, height = 4, units=\"in\", dpi = 300)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "seurat_shifter",
   "language": "R",
   "name": "seurat-shifter"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "4.2.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}