spatialRNAseq_pre_integration.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0af55ff2-b7f1-47ba-882c-b18c88542e12",
   "metadata": {},
   "outputs": [],
   "source": [
    "#loading in the necessary packages\n",
    "##the last line here is a workaround for memory requirements if using Mac\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",
    "library(Seurat.utils)\n",
    "options (future.globals.maxSize = 4000 * 1024^5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8408e97f-1417-4a54-a19e-cb914e7131b7",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Loading in all capture area spaceranger output files using the here package to specify paths on local computer.\n",
    "## Each VGE slide was assigned a number 1-5.\n",
    "##V11Y10-339 = Slide 1 [mycorrhizal]\n",
    "##V11Y10-341 = Slide 2 [mycorrhizal]\n",
    "##V11M22-399 = Slide 3 [mycorrhizal]\n",
    "##V11M22-398 = Slide 4 [control] \n",
    "##V11M22-401 = Slide 5 [control]\n",
    "##Thus, each spatial object for each capture area was named for the capture area letter code and the slide number\n",
    "##As an example, Capture Area \"B1\" from slide V11Y10-339 or slide 1, is named \"B1\", while Capture Area \"B1\" from slide V11Y10-341 or slide 2, is named \"B2\".\n",
    "\n",
    "#B1 from V11Y10-339 - Mycorrhizal\n",
    "B1dataDir1<-here(\"V11Y10-339/B1\")\n",
    "B1<-Load10X_Spatial(data.dir = B1dataDir1, assay='Spatial', filename = \"filtered_feature_bc_matrix.h5\")\n",
    "B1@meta.data['Section_ID'] = 'B1'\n",
    "\n",
    "#A1 from V11Y10-341 - Mycorrhizal \n",
    "A1dataDir2<-here(\"V11Y10-341/A1\")\n",
    "A2<-Load10X_Spatial(data.dir = A1dataDir2, assay='Spatial', filename = \"filtered_feature_bc_matrix.h5\")\n",
    "A2@meta.data['Section_ID'] = 'A1'\n",
    "#B1 from V11Y10-341 - Mycorrhizal \n",
    "B1dataDir2<-here(\"V11Y10-341/B1\")\n",
    "B2<-Load10X_Spatial(data.dir = B1dataDir2, assay='Spatial', filename = \"filtered_feature_bc_matrix.h5\")\n",
    "B2@meta.data['Section_ID'] = 'B1'\n",
    "#C1 from V11Y10-341 - Mycorrhizal \n",
    "C1dataDir2<-here(\"V11Y10-341/C1\")\n",
    "C2<-Load10X_Spatial(data.dir = C1dataDir2, assay='Spatial', filename = \"filtered_feature_bc_matrix.h5\")\n",
    "C2@meta.data['Section_ID'] = 'C1'\n",
    "#D1 from V11Y10-341 - Mycorrhizal \n",
    "D1dataDir2<-here(\"V11Y10-341/D1\")\n",
    "D2<-Load10X_Spatial(data.dir = D1dataDir2, assay='Spatial', filename = \"filtered_feature_bc_matrix.h5\")\n",
    "D2@meta.data['Section_ID'] = 'D1'\n",
    "\n",
    "#A1 from V11M22-399 - Mycorrhizal\n",
    "A1dataDir3<-here(\"V11M22-399_dont/A1\")\n",
    "A3<-Load10X_Spatial(data.dir = A1dataDir3, assay='Spatial', filename = \"filtered_feature_bc_matrix.h5\")\n",
    "A3@meta.data['Section_ID'] = 'A1'\n",
    "#B1 from V11M22-399 - Mycorrhizal\n",
    "B1dataDir3<-here(\"V11M22-399_dont/B1\")\n",
    "B3<-Load10X_Spatial(data.dir = B1dataDir2, assay='Spatial', filename = \"filtered_feature_bc_matrix.h5\")\n",
    "B3@meta.data['Section_ID'] = 'B1'\n",
    "#C1 from V11M22-399 - Mycorrhizal\n",
    "C1dataDir3<-here(\"V11M22-399_dont/C1\")\n",
    "C3<-Load10X_Spatial(data.dir = C1dataDir3, assay='Spatial', filename = \"filtered_feature_bc_matrix.h5\")\n",
    "C3@meta.data['Section_ID'] = 'C1'\n",
    "#C1 from V11M22-399 - Mycorrhizal\n",
    "D1dataDir3<-here(\"V11M22-399_dont/D1\")\n",
    "D3<-Load10X_Spatial(data.dir = D1dataDir2, assay='Spatial', filename = \"filtered_feature_bc_matrix.h5\")\n",
    "D3@meta.data['Section_ID'] = 'D1'\n",
    "\n",
    "#A1 from V11M22-398 - Control\n",
    "A1dataDir4<-here(\"V11M22-398/A1\")\n",
    "A4<-Load10X_Spatial(data.dir = A1dataDir4, assay='Spatial', filename = \"filtered_feature_bc_matrix.h5\")\n",
    "A4@meta.data['Section_ID'] = 'A1'\n",
    "#B1 from V11M22-398 - Control\n",
    "B1dataDir4<-here(\"V11M22-398/B1\")\n",
    "B4<-Load10X_Spatial(data.dir = B1dataDir4, assay='Spatial', filename = \"filtered_feature_bc_matrix.h5\")\n",
    "B4@meta.data['Section_ID'] = 'B1'\n",
    "#C1 from V11M22-398 - Control\n",
    "C1dataDir4<-here(\"V11M22-398/C1\")\n",
    "C4<-Load10X_Spatial(data.dir = C1dataDir4, assay='Spatial', filename = \"filtered_feature_bc_matrix.h5\")\n",
    "C4@meta.data['Section_ID'] = 'C1'\n",
    "#D1 from V11M22-398 - Control\n",
    "D1dataDir4<-here(\"V11M22-398/D1\")\n",
    "D4<-Load10X_Spatial(data.dir = D1dataDir4, assay='Spatial', filename = \"filtered_feature_bc_matrix.h5\")\n",
    "D4@meta.data['Section_ID'] = 'D1'\n",
    "\n",
    "#A1 from V11M22-401 - Control\n",
    "A1dataDir5<-here(\"V11M22-401/A1\")\n",
    "A5<-Load10X_Spatial(data.dir = A1dataDir5, assay='Spatial', filename = \"filtered_feature_bc_matrix.h5\")\n",
    "A5@meta.data['Section_ID'] = 'A1'\n",
    "#B1 from V11M22-401 - Control\n",
    "B1dataDir5<-here(\"V11M22-401/B1\")\n",
    "B5<-Load10X_Spatial(data.dir = B1dataDir5, assay='Spatial', filename = \"filtered_feature_bc_matrix.h5\")\n",
    "B5@meta.data['Section_ID'] = 'B1'\n",
    "#C1 from V11M22-401 - Control\n",
    "C1dataDir5<-here(\"V11M22-401/C1\")\n",
    "C5<-Load10X_Spatial(data.dir = C1dataDir5, assay='Spatial', filename = \"filtered_feature_bc_matrix.h5\")\n",
    "C5@meta.data['Section_ID'] = 'C1'\n",
    "#D1 from V11M22-401 - Control\n",
    "D1dataDir5<-here(\"V11M22-401/D1\")\n",
    "D5<-Load10X_Spatial(data.dir = D1dataDir5, assay='Spatial', filename = \"filtered_feature_bc_matrix.h5\")\n",
    "D5@meta.data['Section_ID'] = 'D1'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4a262c39-c3cd-4812-9880-c9f91a6f1e18",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Now that our capture areas are Seurat objects, we can proceed to filtering and normalization.\n",
    "##From the raw expression values, we selected 13 outlier ribosomal features to remove from all of the datasets prior to merged filtering\n",
    "###removing these features from the counts for each capture area prior to full filtering steps\n",
    "counts <- GetAssayData(B1, assay = \"Spatial\")\n",
    "  counts <- counts[-(which(rownames(counts) %in% c('A17-----------LSU-rRNA-eukarya.17','A17-----------LSU-rRNA-eukarya.16','A17-----------LSU-rRNA-eukarya.7','A17-----------SSU-rRNA-eukarya.5','A17-----------SSU-rRNA-eukarya.4','A17-----------LSU-rRNA-eukarya.20', \"A17-----------SSU-rRNA-eukarya\", 'R.irregularis-RIR-2468470','R.irregularis-RIR-1548030','R.irregularis-RIR-2275770','R.irregularis-RIR-0708530','R.irregularis-RIR-0708570', \"R.irregularis-RIR-1508070\"))),]\n",
    "  B1 <- subset(B1, features = rownames(counts))\n",
    "  \n",
    "counts <- GetAssayData(A2, assay = \"Spatial\")\n",
    "  counts <- counts[-(which(rownames(counts) %in% c('A17-----------LSU-rRNA-eukarya.17','A17-----------LSU-rRNA-eukarya.16','A17-----------LSU-rRNA-eukarya.7','A17-----------SSU-rRNA-eukarya.5','A17-----------SSU-rRNA-eukarya.4','A17-----------LSU-rRNA-eukarya.20', \"A17-----------SSU-rRNA-eukarya\", 'R.irregularis-RIR-2468470','R.irregularis-RIR-1548030','R.irregularis-RIR-2275770','R.irregularis-RIR-0708530','R.irregularis-RIR-0708570', \"R.irregularis-RIR-1508070\"))),]\n",
    "  A2 <- subset(A2, features = rownames(counts))\n",
    "  \n",
    "counts <- GetAssayData(B2, assay = \"Spatial\")\n",
    "  counts <- counts[-(which(rownames(counts) %in% c('A17-----------LSU-rRNA-eukarya.17','A17-----------LSU-rRNA-eukarya.16','A17-----------LSU-rRNA-eukarya.7','A17-----------SSU-rRNA-eukarya.5','A17-----------SSU-rRNA-eukarya.4','A17-----------LSU-rRNA-eukarya.20', \"A17-----------SSU-rRNA-eukarya\", 'R.irregularis-RIR-2468470','R.irregularis-RIR-1548030','R.irregularis-RIR-2275770','R.irregularis-RIR-0708530','R.irregularis-RIR-0708570', \"R.irregularis-RIR-1508070\"))),]\n",
    "  B2 <- subset(B2, features = rownames(counts))\n",
    "  \n",
    "counts <- GetAssayData(C2, assay = \"Spatial\")\n",
    "  counts <- counts[-(which(rownames(counts) %in% c('A17-----------LSU-rRNA-eukarya.17','A17-----------LSU-rRNA-eukarya.16','A17-----------LSU-rRNA-eukarya.7','A17-----------SSU-rRNA-eukarya.5','A17-----------SSU-rRNA-eukarya.4','A17-----------LSU-rRNA-eukarya.20', \"A17-----------SSU-rRNA-eukarya\", 'R.irregularis-RIR-2468470','R.irregularis-RIR-1548030','R.irregularis-RIR-2275770','R.irregularis-RIR-0708530','R.irregularis-RIR-0708570', \"R.irregularis-RIR-1508070\"))),]\n",
    "  C2 <- subset(C2, features = rownames(counts))  \n",
    "  \n",
    "counts <- GetAssayData(D2, assay = \"Spatial\")\n",
    "  counts <- counts[-(which(rownames(counts) %in% c('A17-----------LSU-rRNA-eukarya.17','A17-----------LSU-rRNA-eukarya.16','A17-----------LSU-rRNA-eukarya.7','A17-----------SSU-rRNA-eukarya.5','A17-----------SSU-rRNA-eukarya.4','A17-----------LSU-rRNA-eukarya.20', \"A17-----------SSU-rRNA-eukarya\", 'R.irregularis-RIR-2468470','R.irregularis-RIR-1548030','R.irregularis-RIR-2275770','R.irregularis-RIR-0708530','R.irregularis-RIR-0708570', \"R.irregularis-RIR-1508070\"))),]\n",
    "  D2 <- subset(D2, features = rownames(counts))\n",
    "  \n",
    "counts <- GetAssayData(A3, assay = \"Spatial\")\n",
    "  counts <- counts[-(which(rownames(counts) %in% c('A17-----------LSU-rRNA-eukarya.17','A17-----------LSU-rRNA-eukarya.16','A17-----------LSU-rRNA-eukarya.7','A17-----------SSU-rRNA-eukarya.5','A17-----------SSU-rRNA-eukarya.4','A17-----------LSU-rRNA-eukarya.20', \"A17-----------SSU-rRNA-eukarya\", 'R.irregularis-RIR-2468470','R.irregularis-RIR-1548030','R.irregularis-RIR-2275770','R.irregularis-RIR-0708530','R.irregularis-RIR-0708570', \"R.irregularis-RIR-1508070\"))),]\n",
    "  A3 <- subset(A3, features = rownames(counts))\n",
    "  \n",
    "counts <- GetAssayData(B3, assay = \"Spatial\")\n",
    "  counts <- counts[-(which(rownames(counts) %in% c('A17-----------LSU-rRNA-eukarya.17','A17-----------LSU-rRNA-eukarya.16','A17-----------LSU-rRNA-eukarya.7','A17-----------SSU-rRNA-eukarya.5','A17-----------SSU-rRNA-eukarya.4','A17-----------LSU-rRNA-eukarya.20', \"A17-----------SSU-rRNA-eukarya\", 'R.irregularis-RIR-2468470','R.irregularis-RIR-1548030','R.irregularis-RIR-2275770','R.irregularis-RIR-0708530','R.irregularis-RIR-0708570', \"R.irregularis-RIR-1508070\"))),]\n",
    "  B3 <- subset(B3, features = rownames(counts))\n",
    "  \n",
    "counts <- GetAssayData(C3, assay = \"Spatial\")\n",
    "  counts <- counts[-(which(rownames(counts) %in% c('A17-----------LSU-rRNA-eukarya.17','A17-----------LSU-rRNA-eukarya.16','A17-----------LSU-rRNA-eukarya.7','A17-----------SSU-rRNA-eukarya.5','A17-----------SSU-rRNA-eukarya.4','A17-----------LSU-rRNA-eukarya.20', \"A17-----------SSU-rRNA-eukarya\", 'R.irregularis-RIR-2468470','R.irregularis-RIR-1548030','R.irregularis-RIR-2275770','R.irregularis-RIR-0708530','R.irregularis-RIR-0708570', \"R.irregularis-RIR-1508070\"))),]\n",
    "  C3 <- subset(C3, features = rownames(counts))  \n",
    "  \n",
    "counts <- GetAssayData(D3, assay = \"Spatial\")\n",
    "  counts <- counts[-(which(rownames(counts) %in% c('A17-----------LSU-rRNA-eukarya.17','A17-----------LSU-rRNA-eukarya.16','A17-----------LSU-rRNA-eukarya.7','A17-----------SSU-rRNA-eukarya.5','A17-----------SSU-rRNA-eukarya.4','A17-----------LSU-rRNA-eukarya.20', \"A17-----------SSU-rRNA-eukarya\", 'R.irregularis-RIR-2468470','R.irregularis-RIR-1548030','R.irregularis-RIR-2275770','R.irregularis-RIR-0708530','R.irregularis-RIR-0708570', \"R.irregularis-RIR-1508070\"))),]\n",
    "  D3 <- subset(D3, features = rownames(counts))\n",
    "  \n",
    "counts <- GetAssayData(A4, assay = \"Spatial\")\n",
    "  counts <- counts[-(which(rownames(counts) %in% c('A17-----------LSU-rRNA-eukarya.17','A17-----------LSU-rRNA-eukarya.16','A17-----------LSU-rRNA-eukarya.7','A17-----------SSU-rRNA-eukarya.5','A17-----------SSU-rRNA-eukarya.4','A17-----------LSU-rRNA-eukarya.20', \"A17-----------SSU-rRNA-eukarya\", 'R.irregularis-RIR-2468470','R.irregularis-RIR-1548030','R.irregularis-RIR-2275770','R.irregularis-RIR-0708530','R.irregularis-RIR-0708570', \"R.irregularis-RIR-1508070\"))),]\n",
    "  A4 <- subset(A4, features = rownames(counts))\n",
    "  \n",
    "counts <- GetAssayData(B4, assay = \"Spatial\")\n",
    "  counts <- counts[-(which(rownames(counts) %in% c('A17-----------LSU-rRNA-eukarya.17','A17-----------LSU-rRNA-eukarya.16','A17-----------LSU-rRNA-eukarya.7','A17-----------SSU-rRNA-eukarya.5','A17-----------SSU-rRNA-eukarya.4','A17-----------LSU-rRNA-eukarya.20', \"A17-----------SSU-rRNA-eukarya\", 'R.irregularis-RIR-2468470','R.irregularis-RIR-1548030','R.irregularis-RIR-2275770','R.irregularis-RIR-0708530','R.irregularis-RIR-0708570', \"R.irregularis-RIR-1508070\"))),]\n",
    "  B4 <- subset(B4, features = rownames(counts))\n",
    "  \n",
    "counts <- GetAssayData(C4, assay = \"Spatial\")\n",
    "  counts <- counts[-(which(rownames(counts) %in% c('A17-----------LSU-rRNA-eukarya.17','A17-----------LSU-rRNA-eukarya.16','A17-----------LSU-rRNA-eukarya.7','A17-----------SSU-rRNA-eukarya.5','A17-----------SSU-rRNA-eukarya.4','A17-----------LSU-rRNA-eukarya.20', \"A17-----------SSU-rRNA-eukarya\", 'R.irregularis-RIR-2468470','R.irregularis-RIR-1548030','R.irregularis-RIR-2275770','R.irregularis-RIR-0708530','R.irregularis-RIR-0708570', \"R.irregularis-RIR-1508070\"))),]\n",
    "  C4 <- subset(C4, features = rownames(counts))  \n",
    "  \n",
    "counts <- GetAssayData(D4, assay = \"Spatial\")\n",
    "  counts <- counts[-(which(rownames(counts) %in% c('A17-----------LSU-rRNA-eukarya.17','A17-----------LSU-rRNA-eukarya.16','A17-----------LSU-rRNA-eukarya.7','A17-----------SSU-rRNA-eukarya.5','A17-----------SSU-rRNA-eukarya.4','A17-----------LSU-rRNA-eukarya.20', \"A17-----------SSU-rRNA-eukarya\", 'R.irregularis-RIR-2468470','R.irregularis-RIR-1548030','R.irregularis-RIR-2275770','R.irregularis-RIR-0708530','R.irregularis-RIR-0708570', \"R.irregularis-RIR-1508070\"))),]\n",
    "  D4 <- subset(D4, features = rownames(counts))\n",
    "  \n",
    "counts <- GetAssayData(A5, assay = \"Spatial\")\n",
    "  counts <- counts[-(which(rownames(counts) %in% c('A17-----------LSU-rRNA-eukarya.17','A17-----------LSU-rRNA-eukarya.16','A17-----------LSU-rRNA-eukarya.7','A17-----------SSU-rRNA-eukarya.5','A17-----------SSU-rRNA-eukarya.4','A17-----------LSU-rRNA-eukarya.20', \"A17-----------SSU-rRNA-eukarya\", 'R.irregularis-RIR-2468470','R.irregularis-RIR-1548030','R.irregularis-RIR-2275770','R.irregularis-RIR-0708530','R.irregularis-RIR-0708570', \"R.irregularis-RIR-1508070\"))),]\n",
    "  A5 <- subset(A5, features = rownames(counts))\n",
    "  \n",
    "counts <- GetAssayData(B5, assay = \"Spatial\")\n",
    "  counts <- counts[-(which(rownames(counts) %in% c('A17-----------LSU-rRNA-eukarya.17','A17-----------LSU-rRNA-eukarya.16','A17-----------LSU-rRNA-eukarya.7','A17-----------SSU-rRNA-eukarya.5','A17-----------SSU-rRNA-eukarya.4','A17-----------LSU-rRNA-eukarya.20', \"A17-----------SSU-rRNA-eukarya\", 'R.irregularis-RIR-2468470','R.irregularis-RIR-1548030','R.irregularis-RIR-2275770','R.irregularis-RIR-0708530','R.irregularis-RIR-0708570', \"R.irregularis-RIR-1508070\"))),]\n",
    "  B5 <- subset(B5, features = rownames(counts))\n",
    "  \n",
    "counts <- GetAssayData(C5, assay = \"Spatial\")\n",
    "  counts <- counts[-(which(rownames(counts) %in% c('A17-----------LSU-rRNA-eukarya.17','A17-----------LSU-rRNA-eukarya.16','A17-----------LSU-rRNA-eukarya.7','A17-----------SSU-rRNA-eukarya.5','A17-----------SSU-rRNA-eukarya.4','A17-----------LSU-rRNA-eukarya.20', \"A17-----------SSU-rRNA-eukarya\", 'R.irregularis-RIR-2468470','R.irregularis-RIR-1548030','R.irregularis-RIR-2275770','R.irregularis-RIR-0708530','R.irregularis-RIR-0708570', \"R.irregularis-RIR-1508070\"))),]\n",
    "  C5 <- subset(C5, features = rownames(counts))  \n",
    "  \n",
    "counts <- GetAssayData(D5, assay = \"Spatial\")\n",
    "  counts <- counts[-(which(rownames(counts) %in% c('A17-----------LSU-rRNA-eukarya.17','A17-----------LSU-rRNA-eukarya.16','A17-----------LSU-rRNA-eukarya.7','A17-----------SSU-rRNA-eukarya.5','A17-----------SSU-rRNA-eukarya.4','A17-----------LSU-rRNA-eukarya.20', \"A17-----------SSU-rRNA-eukarya\", 'R.irregularis-RIR-2468470','R.irregularis-RIR-1548030','R.irregularis-RIR-2275770','R.irregularis-RIR-0708530','R.irregularis-RIR-0708570', \"R.irregularis-RIR-1508070\"))),]\n",
    "  D5 <- subset(D5, features = rownames(counts))  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1c5276a8-8eb0-410b-bc13-54989cf8b3e5",
   "metadata": {},
   "outputs": [],
   "source": [
    "#First We want to merge all the capture areas into a single merged object since working with each one is laborious, called that object \"alldata\"\n",
    "##Visium experiments tend to display a high capture rate of ribosomal RNA. \n",
    "###To account for this, we first need to assess the percentage of ribosomal RNA expression across the capture areas.\n",
    "####We used the ribosomal RNA identifier \"[LS]SU-rRNA\" to get the percentage of transcripts with that identifer per voxel and added that percentage to the metadata.\n",
    "#####We were also curious about the breakdown of R. irregularis vs M. truncatula transcripts in each voxel, so we used the \"A17-----------MTR\" and \"R.irregularis-RIR\" identifiers for that\n",
    "alldata <- merge(B1, c(A2, B2, C2, D2, A3, B3, C3, D3, A4, B4, C4, D4, A5, B5, C5, D5), add.cell.ids = c(\"B1\", \"A2\", \"B2\", \"C2\", \"D2\", \"A3\", \"B3\", \"C3\", \"D3\", \"A4\",\"B4\",\"C4\",\"D4\", \"A5\", \"B5\", \"C5\", \"D5\"))\n",
    "#The calculation here is simply the column sum of the matrix present in the counts slot for features belonging to the set divided by the column sum for all features times 100.\n",
    "alldata <- PercentageFeatureSet(\n",
    "  alldata,\n",
    "  pattern = \"R.irregularis-RIR\",\n",
    "  features = NULL,\n",
    "  col.name = \"percent.fungal\",\n",
    "  assay = NULL)\n",
    "\n",
    "alldata <- PercentageFeatureSet(\n",
    "  alldata,\n",
    "  pattern = \"A17-----------MTR\",\n",
    "  features = NULL,\n",
    "  col.name = \"percent.plant\",\n",
    "  assay = NULL)\n",
    "\n",
    "alldata <- PercentageFeatureSet(\n",
    "  alldata,\n",
    "  pattern = \"A17-----------[LS]SU-rRNA\",\n",
    "  features = NULL,\n",
    "  col.name = \"percent.rb\",\n",
    "  assay = NULL)\n",
    "```\n",
    "```{r}\n",
    "#Let's make all of the statistics we're interested in into a list and visualize each of these stats in a volcano plot\n",
    "stats <- c(\"nFeature_Spatial\", \"nCount_Spatial\", \"percent.plant\", \"percent.fungal\", \"percent.rb\")\n",
    "VlnPlot(alldata, group.by = \"orig.ident\", features = stats, pt.size = 0.1, ncol = 3) + NoLegend()\n",
    "```\n",
    "```{r}\n",
    "#let's filter out some of the outliers based on these volcano plots\n",
    "##this may be different depending on capture area which is why we are using the inidividual objects rather than the merged\n",
    "##First filtering for feature counts\n",
    "B1 <- subset(B1, subset = nFeature_Spatial > 100 & nFeature_Spatial < 10000)\n",
    "A2 <- subset(A2, subset = nFeature_Spatial > 100 & nFeature_Spatial < 10000)\n",
    "B2 <- subset(B2, subset = nFeature_Spatial > 100 & nFeature_Spatial < 10000)\n",
    "C2 <- subset(C2, subset = nFeature_Spatial > 100 & nFeature_Spatial < 10000)\n",
    "D2 <- subset(D2, subset = nFeature_Spatial > 100 & nFeature_Spatial < 10000)\n",
    "A3 <- subset(A3, subset = nFeature_Spatial > 100 & nFeature_Spatial < 10000)\n",
    "B3 <- subset(B3, subset = nFeature_Spatial > 100 & nFeature_Spatial < 10000)\n",
    "C3 <- subset(C3, subset = nFeature_Spatial > 100 & nFeature_Spatial < 10000)\n",
    "D3 <- subset(D3, subset = nFeature_Spatial > 100 & nFeature_Spatial < 10000)\n",
    "A4 <- subset(A4, subset = nFeature_Spatial > 100 & nFeature_Spatial < 10000)\n",
    "B4 <- subset(B4, subset = nFeature_Spatial > 100 & nFeature_Spatial < 10000)\n",
    "D4 <- subset(C4, subset = nFeature_Spatial > 100 & nFeature_Spatial < 10000)\n",
    "A4 <- subset(D4, subset = nFeature_Spatial > 100 & nFeature_Spatial < 10000)\n",
    "A5 <- subset(A5, subset = nFeature_Spatial > 100 & nFeature_Spatial < 10000)\n",
    "B5 <- subset(B5, subset = nFeature_Spatial > 100 & nFeature_Spatial < 10000)\n",
    "C5 <- subset(C5, subset = nFeature_Spatial > 100 & nFeature_Spatial < 10000)\n",
    "D5 <- subset(D5, subset = nFeature_Spatial > 100 & nFeature_Spatial < 10000)\n",
    "\n",
    "#now umi counts\n",
    "B1 <- subset(B1, subset = nCount_Spatial > 100 & nCount_Spatial < 60000)\n",
    "A2 <- subset(A2, subset = nCount_Spatial > 100 & nCount_Spatial < 60000)\n",
    "B2 <- subset(B2, subset = nCount_Spatial > 100 & nCount_Spatial < 60000)\n",
    "C2 <- subset(C2, subset = nCount_Spatial > 100 & nCount_Spatial < 60000)\n",
    "D2 <- subset(D2, subset = nCount_Spatial > 100 & nCount_Spatial < 60000)\n",
    "A3 <- subset(A3, subset = nCount_Spatial > 100 & nCount_Spatial < 60000)\n",
    "B3 <- subset(B3, subset = nCount_Spatial > 100 & nCount_Spatial < 60000)\n",
    "C3 <- subset(C3, subset = nCount_Spatial > 100 & nCount_Spatial < 60000)\n",
    "D3 <- subset(D3, subset = nCount_Spatial > 100 & nCount_Spatial < 60000)\n",
    "A4 <- subset(A4, subset = nCount_Spatial > 100 & nCount_Spatial < 60000)\n",
    "B4 <- subset(B4, subset = nCount_Spatial > 100 & nCount_Spatial < 60000)\n",
    "D4 <- subset(C4, subset = nCount_Spatial > 100 & nCount_Spatial < 60000)\n",
    "A4 <- subset(D4, subset = nCount_Spatial > 100 & nCount_Spatial < 60000)\n",
    "A5 <- subset(A5, subset = nCount_Spatial > 100 & nCount_Spatial < 60000)\n",
    "B5 <- subset(B5, subset = nCount_Spatial > 100 & nCount_Spatial < 60000)\n",
    "C5 <- subset(C5, subset = nCount_Spatial > 100 & nCount_Spatial < 60000)\n",
    "D5 <- subset(D5, subset = nCount_Spatial > 100 & nCount_Spatial < 60000)\n",
    "\n",
    "#now ribosomal RNA\n",
    "##for this we have to add the PercentFeatureSet to each capture area first since we only added it to the merged dataset\n",
    "B1 <- PercentageFeatureSet(\n",
    "  B1,\n",
    "  pattern = \"A17-----------[LS]SU-rRNA\",\n",
    "  features = NULL,\n",
    "  col.name = \"percent.rb\",\n",
    "  assay = NULL)\n",
    "\n",
    "A2 <- PercentageFeatureSet(\n",
    "  A2,\n",
    "  pattern = \"A17-----------[LS]SU-rRNA\",\n",
    "  features = NULL,\n",
    "  col.name = \"percent.rb\",\n",
    "  assay = NULL)\n",
    "\n",
    "B2 <- PercentageFeatureSet(\n",
    "  B2,\n",
    "  pattern = \"A17-----------[LS]SU-rRNA\",\n",
    "  features = NULL,\n",
    "  col.name = \"percent.rb\",\n",
    "  assay = NULL)\n",
    "\n",
    "C2 <- PercentageFeatureSet(\n",
    "  C2,\n",
    "  pattern = \"A17-----------[LS]SU-rRNA\",\n",
    "  features = NULL,\n",
    "  col.name = \"percent.rb\",\n",
    "  assay = NULL)\n",
    "\n",
    "D2 <- PercentageFeatureSet(\n",
    "  D2,\n",
    "  pattern = \"A17-----------[LS]SU-rRNA\",\n",
    "  features = NULL,\n",
    "  col.name = \"percent.rb\",\n",
    "  assay = NULL)\n",
    "\n",
    "A3 <- PercentageFeatureSet(\n",
    "  A3,\n",
    "  pattern = \"A17-----------[LS]SU-rRNA\",\n",
    "  features = NULL,\n",
    "  col.name = \"percent.rb\",\n",
    "  assay = NULL)\n",
    "\n",
    "B3 <- PercentageFeatureSet(\n",
    "  B3,\n",
    "  pattern = \"A17-----------[LS]SU-rRNA\",\n",
    "  features = NULL,\n",
    "  col.name = \"percent.rb\",\n",
    "  assay = NULL)\n",
    "\n",
    "C3 <- PercentageFeatureSet(\n",
    "  C3,\n",
    "  pattern = \"A17-----------[LS]SU-rRNA\",\n",
    "  features = NULL,\n",
    "  col.name = \"percent.rb\",\n",
    "  assay = NULL)\n",
    "\n",
    "D3 <- PercentageFeatureSet(\n",
    "  D3,\n",
    "  pattern = \"A17-----------[LS]SU-rRNA\",\n",
    "  features = NULL,\n",
    "  col.name = \"percent.rb\",\n",
    "  assay = NULL)\n",
    "\n",
    "A4 <- PercentageFeatureSet(\n",
    "  A4,\n",
    "  pattern = \"A17-----------[LS]SU-rRNA\",\n",
    "  features = NULL,\n",
    "  col.name = \"percent.rb\",\n",
    "  assay = NULL)\n",
    "\n",
    "B4 <- PercentageFeatureSet(\n",
    "  B4,\n",
    "  pattern = \"A17-----------[LS]SU-rRNA\",\n",
    "  features = NULL,\n",
    "  col.name = \"percent.rb\",\n",
    "  assay = NULL)\n",
    "\n",
    "C4 <- PercentageFeatureSet(\n",
    "  C4,\n",
    "  pattern = \"A17-----------[LS]SU-rRNA\",\n",
    "  features = NULL,\n",
    "  col.name = \"percent.rb\",\n",
    "  assay = NULL)\n",
    "\n",
    "D4 <- PercentageFeatureSet(\n",
    "  D4,\n",
    "  pattern = \"A17-----------[LS]SU-rRNA\",\n",
    "  features = NULL,\n",
    "  col.name = \"percent.rb\",\n",
    "  assay = NULL)\n",
    "\n",
    "A5 <- PercentageFeatureSet(\n",
    "  A5,\n",
    "  pattern = \"A17-----------[LS]SU-rRNA\",\n",
    "  features = NULL,\n",
    "  col.name = \"percent.rb\",\n",
    "  assay = NULL)\n",
    "\n",
    "B5 <- PercentageFeatureSet(\n",
    "  B5,\n",
    "  pattern = \"A17-----------[LS]SU-rRNA\",\n",
    "  features = NULL,\n",
    "  col.name = \"percent.rb\",\n",
    "  assay = NULL)\n",
    "\n",
    "C5 <- PercentageFeatureSet(\n",
    "  D5,\n",
    "  pattern = \"A17-----------[LS]SU-rRNA\",\n",
    "  features = NULL,\n",
    "  col.name = \"percent.rb\",\n",
    "  assay = NULL)\n",
    "\n",
    "D5 <- PercentageFeatureSet(\n",
    "  D5,\n",
    "  pattern = \"A17-----------[LS]SU-rRNA\",\n",
    "  features = NULL,\n",
    "  col.name = \"percent.rb\",\n",
    "  assay = NULL)\n",
    "\n",
    "#adding in metadata per capture area\n",
    "B1 <- PercentageFeatureSet(\n",
    "  B1,\n",
    "  pattern = \"R.irregularis-RIR\",\n",
    "  features = NULL,\n",
    "  col.name = \"percent.ri\",\n",
    "  assay = NULL)\n",
    "\n",
    "A2 <- PercentageFeatureSet(\n",
    "  A2,\n",
    "  pattern = \"R.irregularis-RIR\",\n",
    "  features = NULL,\n",
    "  col.name = \"percent.ri\",\n",
    "  assay = NULL)\n",
    "\n",
    "B2 <- PercentageFeatureSet(\n",
    "  B2,\n",
    "  pattern = \"R.irregularis-RIR\",\n",
    "  features = NULL,\n",
    "  col.name = \"percent.ri\",\n",
    "  assay = NULL)\n",
    "\n",
    "C2 <- PercentageFeatureSet(\n",
    "  C2,\n",
    "  pattern = \"R.irregularis-RIR\",\n",
    "  features = NULL,\n",
    "  col.name = \"percent.ri\",\n",
    "  assay = NULL)\n",
    "\n",
    "D2 <- PercentageFeatureSet(\n",
    "  D2,\n",
    "  pattern = \"R.irregularis-RIR\",\n",
    "  features = NULL,\n",
    "  col.name = \"percent.ri\",\n",
    "  assay = NULL)\n",
    "\n",
    "A3 <- PercentageFeatureSet(\n",
    "  A3,\n",
    "  pattern = \"R.irregularis-RIR\",\n",
    "  features = NULL,\n",
    "  col.name = \"percent.ri\",\n",
    "  assay = NULL)\n",
    "\n",
    "B3 <- PercentageFeatureSet(\n",
    "  B3,\n",
    "  pattern = \"R.irregularis-RIR\",\n",
    "  features = NULL,\n",
    "  col.name = \"percent.ri\",\n",
    "  assay = NULL)\n",
    "\n",
    "C3 <- PercentageFeatureSet(\n",
    "  C3,\n",
    "  pattern = \"R.irregularis-RIR\",\n",
    "  features = NULL,\n",
    "  col.name = \"percent.ri\",\n",
    "  assay = NULL)\n",
    "\n",
    "D3 <- PercentageFeatureSet(\n",
    "  D3,\n",
    "  pattern = \"R.irregularis-RIR\",\n",
    "  features = NULL,\n",
    "  col.name = \"percent.ri\",\n",
    "  assay = NULL)\n",
    "\n",
    "A4 <- PercentageFeatureSet(\n",
    "  A4,\n",
    "  pattern = \"R.irregularis-RIR\",\n",
    "  features = NULL,\n",
    "  col.name = \"percent.ri\",\n",
    "  assay = NULL)\n",
    "\n",
    "B4 <- PercentageFeatureSet(\n",
    "  B4,\n",
    "  pattern = \"R.irregularis-RIR\",\n",
    "  features = NULL,\n",
    "  col.name = \"percent.ri\",\n",
    "  assay = NULL)\n",
    "\n",
    "C4 <- PercentageFeatureSet(\n",
    "  C4,\n",
    "  pattern = \"R.irregularis-RIR\",\n",
    "  features = NULL,\n",
    "  col.name = \"percent.ri\",\n",
    "  assay = NULL)\n",
    "\n",
    "D4 <- PercentageFeatureSet(\n",
    "  D4,\n",
    "  pattern = \"R.irregularis-RIR\",\n",
    "  features = NULL,\n",
    "  col.name = \"percent.ri\",\n",
    "  assay = NULL)\n",
    "\n",
    "A5 <- PercentageFeatureSet(\n",
    "  A5,\n",
    "  pattern = \"R.irregularis-RIR\",\n",
    "  features = NULL,\n",
    "  col.name = \"percent.ri\",\n",
    "  assay = NULL)\n",
    "\n",
    "B5 <- PercentageFeatureSet(\n",
    "  B5,\n",
    "  pattern = \"R.irregularis-RIR\",\n",
    "  features = NULL,\n",
    "  col.name = \"percent.ri\",\n",
    "  assay = NULL)\n",
    "\n",
    "C5 <- PercentageFeatureSet(\n",
    "  D5,\n",
    "  pattern = \"R.irregularis-RIR\",\n",
    "  features = NULL,\n",
    "  col.name = \"percent.ri\",\n",
    "  assay = NULL)\n",
    "\n",
    "D5 <- PercentageFeatureSet(\n",
    "  D5,\n",
    "  pattern = \"R.irregularis-RIR\",\n",
    "  features = NULL,\n",
    "  col.name = \"percent.ri\",\n",
    "  assay = NULL)\n",
    "\n",
    "#filtering to keep only cells with below 30% rRNA content, which according to the plot kept the majority of cells and got rid of outliers\n",
    "selected_ribo <- WhichCells(B1, expression = percent.rb < 30)\n",
    "B1 <- subset(B1, cells = selected_ribo)\n",
    "selected_ribo <- WhichCells(A2, expression = percent.rb < 30)\n",
    "A2 <- subset(A2, cells = selected_ribo)\n",
    "selected_ribo <- WhichCells(B2, expression = percent.rb < 30)\n",
    "B2 <- subset(B2, cells = selected_ribo)\n",
    "selected_ribo <- WhichCells(C2, expression = percent.rb < 30)\n",
    "C2 <- subset(C2, cells = selected_ribo)\n",
    "selected_ribo <- WhichCells(D2, expression = percent.rb < 30)\n",
    "D2 <- subset(D2, cells = selected_ribo)\n",
    "selected_ribo <- WhichCells(A3, expression = percent.rb < 30)\n",
    "A3 <- subset(A3, cells = selected_ribo)\n",
    "selected_ribo <- WhichCells(B3, expression = percent.rb < 30)\n",
    "B3 <- subset(B3, cells = selected_ribo)\n",
    "selected_ribo <- WhichCells(C3, expression = percent.rb < 30)\n",
    "C3 <- subset(C3, cells = selected_ribo)\n",
    "selected_ribo <- WhichCells(D3, expression = percent.rb < 30)\n",
    "D3 <- subset(D3, cells = selected_ribo)\n",
    "selected_ribo <- WhichCells(A4, expression = percent.rb < 30)\n",
    "A4 <- subset(A4, cells = selected_ribo)\n",
    "selected_ribo <- WhichCells(B4, expression = percent.rb < 30)\n",
    "B4 <- subset(B4, cells = selected_ribo)\n",
    "selected_ribo <- WhichCells(C4, expression = percent.rb < 30)\n",
    "C4 <- subset(C4, cells = selected_ribo)\n",
    "selected_ribo <- WhichCells(D4, expression = percent.rb < 30)\n",
    "D4 <- subset(D4, cells = selected_ribo)\n",
    "selected_ribo <- WhichCells(A5, expression = percent.rb < 30)\n",
    "A5 <- subset(A5, cells = selected_ribo)\n",
    "selected_ribo <- WhichCells(B5, expression = percent.rb < 30)\n",
    "B5 <- subset(B5, cells = selected_ribo)\n",
    "selected_ribo <- WhichCells(C5, expression = percent.rb < 30)\n",
    "C5 <- subset(C5, cells = selected_ribo)\n",
    "selected_ribo <- WhichCells(D5, expression = percent.rb < 30)\n",
    "D5 <- subset(D5, cells = selected_ribo)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "83f8987c-6085-452d-9495-987e5580e987",
   "metadata": {},
   "outputs": [],
   "source": [
    "#After filtering is done, let's normalize the data using SCT Transform\n",
    "#SCT builds regularized negative binomial models of gene expression in order to account for technical artifacts while preserving biological variance. \n",
    "#SCTransform normalizes the data, detects high-variance features, and stores the data in the SCT assay. No need to scale before clustering when using.\n",
    "B1 <- SCTransform(B1, assay = \"Spatial\", verbose = FALSE)\n",
    "A2 <- SCTransform(A2, assay = \"Spatial\", verbose = FALSE)\n",
    "B2 <- SCTransform(B2, assay = \"Spatial\", verbose = FALSE)\n",
    "C2 <- SCTransform(C2, assay = \"Spatial\", verbose = FALSE)\n",
    "D2 <- SCTransform(D2, assay = \"Spatial\", verbose = FALSE)\n",
    "A3 <- SCTransform(A3, assay = \"Spatial\", verbose = FALSE)\n",
    "B3 <- SCTransform(B3, assay = \"Spatial\", verbose = FALSE)\n",
    "C3 <- SCTransform(C3, assay = \"Spatial\", verbose = FALSE)\n",
    "D3 <- SCTransform(D3, assay = \"Spatial\", verbose = FALSE)\n",
    "A4 <- SCTransform(A4, assay = \"Spatial\", verbose = FALSE)\n",
    "B4 <- SCTransform(B4, assay = \"Spatial\", verbose = FALSE)\n",
    "C4 <- SCTransform(C4, assay = \"Spatial\", verbose = FALSE)\n",
    "D4 <- SCTransform(D4, assay = \"Spatial\", verbose = FALSE)\n",
    "A5 <- SCTransform(A5, assay = \"Spatial\", verbose = FALSE)\n",
    "B5 <- SCTransform(B5, assay = \"Spatial\", verbose = FALSE)\n",
    "C5 <- SCTransform(C5, assay = \"Spatial\", verbose = FALSE)\n",
    "D5 <- SCTransform(D5, assay = \"Spatial\", verbose = FALSE)\n",
    "\n",
    "#Now we can perform dimensionality reduction and clustering for each capture area.\n",
    "B1 <- RunPCA(B1, assay = \"SCT\", verbose = FALSE)\n",
    "A2 <- RunPCA(A2, assay = \"SCT\", verbose = FALSE)\n",
    "B2 <- RunPCA(B2, assay = \"SCT\", verbose = FALSE)\n",
    "C2 <- RunPCA(C2, assay = \"SCT\", verbose = FALSE)\n",
    "D2 <- RunPCA(D2, assay = \"SCT\", verbose = FALSE)\n",
    "A3 <- RunPCA(A3, assay = \"SCT\", verbose = FALSE)\n",
    "B3 <- RunPCA(B3, assay = \"SCT\", verbose = FALSE)\n",
    "C3 <- RunPCA(C3, assay = \"SCT\", verbose = FALSE)\n",
    "D3 <- RunPCA(D3, assay = \"SCT\", verbose = FALSE)\n",
    "A4 <- RunPCA(A4, assay = \"SCT\", verbose = FALSE)\n",
    "B4 <- RunPCA(B4, assay = \"SCT\", verbose = FALSE)\n",
    "C4 <- RunPCA(C4, assay = \"SCT\", verbose = FALSE)\n",
    "D4 <- RunPCA(D4, assay = \"SCT\", verbose = FALSE)\n",
    "A5 <- RunPCA(A5, assay = \"SCT\", verbose = FALSE)\n",
    "B5 <- RunPCA(B5, assay = \"SCT\", verbose = FALSE)\n",
    "C5 <- RunPCA(C5, assay = \"SCT\", verbose = FALSE)\n",
    "D5 <- RunPCA(D5, assay = \"SCT\", verbose = FALSE)\n",
    "\n",
    "#Let's see how well the data is captured by the PCs\n",
    "ElbowPlot(B1)\n",
    "ElbowPlot(A2)\n",
    "ElbowPlot(B2)\n",
    "ElbowPlot(C2)\n",
    "ElbowPlot(D2)\n",
    "ElbowPlot(A3)\n",
    "ElbowPlot(B3)\n",
    "ElbowPlot(C3)\n",
    "ElbowPlot(D3)\n",
    "ElbowPlot(A4)\n",
    "ElbowPlot(B4)\n",
    "ElbowPlot(C4)\n",
    "ElbowPlot(D4)\n",
    "ElbowPlot(A5)\n",
    "ElbowPlot(B5)\n",
    "ElbowPlot(C5)\n",
    "ElbowPlot(D5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f0c67123-68ee-4b62-8204-3c64127238b7",
   "metadata": {},
   "outputs": [],
   "source": [
    "#As in PhenoGraph, we first construct a KNN graph based on the euclidean distance in PCA space, and refine the edge weights between any two cells based on the shared overlap in their local neighborhoods (Jaccard similarity). \n",
    "#This step is performed using the FindNeighbors() function, and takes as input the previously defined dimensionality of the dataset (first 20 PCs).\n",
    "B1 <- FindNeighbors(B1, reduction = \"pca\", dims = 1:30)\n",
    "B1 <- FindClusters(B1, verbose = FALSE)\n",
    "B1 <- RunUMAP(B1, reduction= \"pca\", dims = 1:30)\n",
    "\n",
    "\n",
    "A2 <- FindNeighbors(A2, reduction = \"pca\", dims = 1:30)\n",
    "A2 <- FindClusters(A2, verbose = FALSE)\n",
    "A2 <- RunUMAP(A2, reduction = \"pca\", dims = 1:30)\n",
    "B2 <- FindNeighbors(B2, reduction = \"pca\", dims = 1:30)\n",
    "B2 <- FindClusters(B2, verbose = FALSE)\n",
    "B2 <- RunUMAP(B2, reduction = \"pca\", dims = 1:30)\n",
    "C2 <- FindNeighbors(C2, reduction = \"pca\", dims = 1:30)\n",
    "C2 <- FindClusters(C2, verbose = FALSE)\n",
    "C2 <- RunUMAP(C2, reduction = \"pca\", dims = 1:30)\n",
    "D2 <- FindNeighbors(D2, reduction = \"pca\", dims = 1:30)\n",
    "D2 <- FindClusters(D2, verbose = FALSE)\n",
    "D2 <- RunUMAP(D2, reduction = \"pca\", dims = 1:30)\n",
    "\n",
    "\n",
    "A3 <- FindNeighbors(A3, reduction = \"pca\", dims = 1:30)\n",
    "A3 <- FindClusters(A3, verbose = FALSE)\n",
    "A3 <- RunUMAP(A3, reduction = \"pca\", dims = 1:30)\n",
    "B3 <- FindNeighbors(B3, reduction = \"pca\", dims = 1:30)\n",
    "B3 <- FindClusters(B3, verbose = FALSE)\n",
    "B3 <- RunUMAP(B3, reduction = \"pca\", dims = 1:30)\n",
    "C3 <- FindNeighbors(C3, reduction = \"pca\", dims = 1:30)\n",
    "C3 <- FindClusters(C3, verbose = FALSE)\n",
    "C3 <- RunUMAP(C3, reduction = \"pca\", dims = 1:30)\n",
    "D3 <- FindNeighbors(D3, reduction = \"pca\", dims = 1:30)\n",
    "D3 <- FindClusters(D3, verbose = FALSE)\n",
    "D3 <- RunUMAP(D3, reduction = \"pca\", dims = 1:30)\n",
    "\n",
    "\n",
    "A4 <- FindNeighbors(A4, reduction = \"pca\", dims = 1:30)\n",
    "A4 <- FindClusters(A4, verbose = FALSE)\n",
    "A4 <- RunUMAP(A4, reduction = \"pca\", dims = 1:30)\n",
    "B4 <- FindNeighbors(B4, reduction = \"pca\", dims = 1:30)\n",
    "B4 <- FindClusters(B4, verbose = FALSE)\n",
    "B4 <- RunUMAP(B4, reduction = \"pca\", dims = 1:30)\n",
    "C4 <- FindNeighbors(C4, reduction = \"pca\", dims = 1:30)\n",
    "C4 <- FindClusters(C4, verbose = FALSE)\n",
    "C4 <- RunUMAP(C4, reduction = \"pca\", dims = 1:30)\n",
    "D4 <- FindNeighbors(D4, reduction = \"pca\", dims = 1:30)\n",
    "D4 <- FindClusters(D4, verbose = FALSE)\n",
    "D4 <- RunUMAP(D4, reduction = \"pca\", dims = 1:30)\n",
    "\n",
    "A5 <- FindNeighbors(A5, reduction = \"pca\", dims = 1:30)\n",
    "A5 <- FindClusters(A5, verbose = FALSE)\n",
    "A5 <- RunUMAP(A5, reduction = \"pca\", dims = 1:30)\n",
    "B5 <- FindNeighbors(B5, reduction = \"pca\", dims = 1:30)\n",
    "B5 <- FindClusters(B5, verbose = FALSE)\n",
    "B5 <- RunUMAP(B5, reduction = \"pca\", dims = 1:30)\n",
    "C5 <- FindNeighbors(C5, reduction = \"pca\", dims = 1:30)\n",
    "C5 <- FindClusters(C5, verbose = FALSE)\n",
    "C5 <- RunUMAP(C5, reduction = \"pca\", dims = 1:30)\n",
    "D5 <- FindNeighbors(D5, reduction = \"pca\", dims = 1:30)\n",
    "D5 <- FindClusters(D5, verbose = FALSE)\n",
    "D5 <- RunUMAP(D5, reduction = \"pca\", dims = 1:30)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e939295b-b069-4be1-8a98-b072b8677ac4",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Visualizing the PCA analysis results in UMAP space\n",
    "plot00 <- DimPlot(B1, reduction = \"umap\", label = TRUE)\n",
    "plot01 <- DimPlot(A2, reduction = \"umap\", label = TRUE)\n",
    "plot02 <- DimPlot(B2, reduction = \"umap\", label = TRUE)\n",
    "plot03 <- DimPlot(C2, reduction = \"umap\", label = TRUE)\n",
    "plot04 <- DimPlot(D2, reduction = \"umap\", label = TRUE)\n",
    "plot05 <- DimPlot(A3, reduction = \"umap\", label = TRUE)\n",
    "plot06 <- DimPlot(B3, reduction = \"umap\", label = TRUE)\n",
    "plot07 <- DimPlot(C3, reduction = \"umap\", label = TRUE)\n",
    "plot08 <- DimPlot(D3, reduction = \"umap\", label = TRUE)\n",
    "plot09 <- DimPlot(A4, reduction = \"umap\", label = TRUE)\n",
    "plot10 <- DimPlot(B4, reduction = \"umap\", label = TRUE)\n",
    "plot11 <- DimPlot(C4, reduction = \"umap\", label = TRUE)\n",
    "plot12 <- DimPlot(D4, reduction = \"umap\", label = TRUE)\n",
    "plot13 <- DimPlot(A5, reduction = \"umap\", label = TRUE)\n",
    "plot14 <- DimPlot(B5, reduction = \"umap\", label = TRUE)\n",
    "plot15 <- DimPlot(C5, reduction = \"umap\", label = TRUE)\n",
    "plot16 <- DimPlot(D5, reduction = \"umap\", label = TRUE)\n",
    "\n",
    "plot00\n",
    "plot01 + plot02 + plot03 + plot04\n",
    "plot05 + plot06 + plot07 + plot08\n",
    "plot09 + plot10 + plot11 + plot12\n",
    "plot13 + plot14 + plot15 + plot16"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "24552420-f3e9-4948-b160-78cc8b4c3084",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Now that we have performed all quality control steps, we want to save each of these as new objects in .RDS formats for easy recall in integration\n",
    "saveRDS(B1, \"B1.rds\")\n",
    "saveRDS(A2, \"A2.rds\")\n",
    "saveRDS(B2, \"B2.rds\")\n",
    "saveRDS(C2, \"C2.rds\")\n",
    "saveRDS(D2, \"D2.rds\")\n",
    "saveRDS(A3, \"A3.rds\")\n",
    "saveRDS(B3, \"B3.rds\")\n",
    "saveRDS(C3, \"C3.rds\")\n",
    "saveRDS(D3, \"D3.rds\")\n",
    "saveRDS(A4, \"A4.rds\")\n",
    "saveRDS(B4, \"B4.rds\")\n",
    "saveRDS(C4, \"C4.rds\")\n",
    "saveRDS(D4, \"D4.rds\")\n",
    "saveRDS(A5, \"A5.rds\")\n",
    "saveRDS(B5, \"B5.rds\")\n",
    "saveRDS(C5, \"C5.rds\")\n",
    "saveRDS(D5, \"D5.rds\")"
   ]
  }
 ],
 "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
}