Thank you for the nice words on our paper. From what I understood your have question regarding sample integration. We actually used the Seurat v3 pipeline which relies on anchor cells to compute batch-correction vectors between datasets. https://satijalab.org/seurat/archive/v3.0/integration.html
You can find the Seurat code below. I have also sent you this script via email since bioprotocol doesn't allow to send R scripts.
library(Seurat)
library(gplots)
library(paletteer)
library(tidyverse)
# First we generate one seurat object for each patient and put them in a list
eleven_cd3.list <- c(P34_Tumor,P35_Tumor,P42_Tumor,P43_Tumor,P46_Tumor,P47_Tumor,P55_Tumor, P57_Tumor, P58_Tumor, P60_Tumor, P61_Tumor)
for (i in 1:length(x = eleven_cd3.list)) {
eleven_cd3.list[[i]] <- NormalizeData(object = eleven_cd3.list[[i]], verbose = T)
eleven_cd3.list[[i]] <- FindVariableFeatures(object = eleven_cd3.list[[i]],
selection.method = "vst", nfeatures = 2000, verbose = T)
}
eleven_cd3.anchors <- FindIntegrationAnchors(object.list = eleven_cd3.list, dims = 1:30)
eleven.tils.cd3.integrated <- IntegrateData(anchorset = eleven_cd3.anchors, dims = 1:30)
# switch to integrated assay. The variable features of this assay are automatically set during IntegrateData
DefaultAssay(object = eleven.tils.cd3.integrated) <- "integrated"
# Run the standard workflow for visualization and clustering
eleven.tils.cd3.integrated <- ScaleData(object = eleven.tils.cd3.integrated, verbose = FALSE)
eleven.tils.cd3.integrated <- RunPCA(object = eleven.tils.cd3.integrated, npcs = 30, verbose = FALSE)
eleven.tils.cd3.integrated <- FindNeighbors(object = eleven.tils.cd3.integrated, dims = 1:30)
eleven.tils.cd3.integrated <- FindClusters(object = eleven.tils.cd3.integrated, resolution = 1.2)
eleven.tils.cd3.integrated <- RunUMAP(object = eleven.tils.cd3.integrated, reduction = "pca",
dims = 1:30)
DimPlot(object = eleven.tils.cd3.integrated, reduction = "umap", group.by = "orig.ident", pt.size = 1)
DimPlot(object = eleven.tils.cd3.integrated, reduction = "umap", pt.size = 1, label = T, repel = T)
FeaturePlot(eleven.tils.cd3.integrated, "CD3D", pt.size = 1, min.cutoff = "q3",order = T, max.cutoff = "q97")