search

stuart-lab / signac

R toolkit for the analysis of single-cell chromatin data
https://stuartlab.org/signac/
Other
331 stars 88 forks source link

Integrating multiple multiome samples: unable to combine integrated ATAC and RNA into one object. #1574

Closed hannahvanm closed 10 months ago

hannahvanm commented 11 months ago

First, thank you so much for all of the help thus far on using Signac packages.  I'm working with two multiome samples that I want to integrate for analysis. Looking through other github responses, it seems people have had success integrating the RNA and ATAC seperately, then using FindMultiModalNeighbors to carry on with WNN analysis https://github.com/stuart-lab/signac/discussions/438 I integrated our RNA using IntegrateLayers for CCA Integration, then integrated our ATAC using the Signac scATAC-seq data integration vignette.

I'm struggling to successfully integrate the combined RNA and ATAC into one object so that I can go do FindMultiModalNeighbors, using this code: T20_T24_integrated[['ATAC']] <- T20_T24_ATAC_integrated[['ATAC']] T20_T24_integrated[['integrated_lsi']] <- T20_T24_ATAC_integrated[['integrated_lsi']]

Here's the error I get: Error in[[<-`: ! Cannot add new cells with [[<- Backtrace:

  1. methods (local) [[<-(*tmp*, "ATAC", value = <ChrmtnAs[,161566]>)
  2. SeuratObject (local) [[<-(*tmp*, "ATAC", value = <ChrmtnAs[,161566]>)`

I think part of it the issue is that the way cells are called in the ATAC vignette is not matching the actual number of cells in our sample- it far exceeds our cell count. Is there a way to have the cells from the RNA assay be used as cells in the ATAC (for example when creating the feature matrix) or any other fix someone could suggest? I've tried a few different ways already, but so far I'm still unable to combine the integrated RNA and ATAC objects into one.

Dexiao211 commented 11 months ago

Hi @hannahvanm, I am currently working on a project similar to yours, focusing on the integration of four multiome samples. I noticed that you might have encountered a similar issue to the one I am facing. Would it be possible for you to share your codes? Thank you for considering my request.

hannahvanm commented 11 months ago

Hi @Dexiao211 Sure, I'm happy to share (chaotic as it might be) 

T20rna_counts <- T20$`Gene Expression`
T24rna_counts <- T24$`Gene Expression`

T20_CNC <- CreateSeuratObject(counts = T20rna_counts, project = "T20_CNC", min.cells = 3, min.features = 100)
T24_TNC <- CreateSeuratObject(counts = T24rna_counts, project = "T24_CNC", min.cells = 3, min.features = 100)

T20_CNC <- RenameCells(T20_CNC, add.cell.id = "T20")
T24_TNC <- RenameCells(T24_TNC, add.cell.id = "T24")
head(T20_CNC)

T20_T24_RNA_merged <- merge(T20_CNC, y = T24_TNC, project = "T20_T24_RNA_merged")
Cells(T20_T24_RNA_merged)

T20_T24_RNA_merged[["percent.mt"]] <- PercentageFeatureSet(T20_T24_RNA_merged, pattern = "^mt-")
VlnPlot(T20_T24_RNA_merged, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
T20_T24_RNA_merged <- subset(T20_T24_RNA_merged, subset = nFeature_RNA < 3100 & nCount_RNA < 10000 & percent.mt < 3)

T20_T24_RNA_merged <- NormalizeData(T20_T24_RNA_merged)
T20_T24_RNA_merged <- FindVariableFeatures(T20_T24_RNA_merged)
T20_T24_RNA_merged <- ScaleData(T20_T24_RNA_merged)
T20_T24_RNA_merged <- RunPCA(T20_T24_RNA_merged)
ElbowPlot(T20_T24_RNA_merged)
T20_T24_RNA_merged <- FindNeighbors(T20_T24_RNA_merged, dims = 1:15, reduction = "pca")
T20_T24_RNA_merged <- FindClusters(T20_T24_RNA_merged, resolution = 1, cluster.name = "unintegrated_clusters")

T20_T24_RNA_merged <- RunUMAP(T20_T24_RNA_merged, dims = 1:15, reduction = "pca", reduction.name = "umap")
DimPlot(T20_T24_RNA_merged, reduction = "umap", group.by = c("orig.ident", "seurat_clusters"))

T20_T24_integrated <- IntegrateLayers(object = T20_T24_RNA_merged, method = CCAIntegration, orig.reduction = "pca", new.reduction = "integrated.cca",
    verbose = FALSE)

T20_T24_integrated[["RNA"]] <- JoinLayers(T20_T24_integrated[["RNA"]])
T20_T24_integrated <- FindNeighbors(T20_T24_integrated, reduction = "integrated.cca", dims = 1:20)
T20_T24_integrated <- FindClusters(T20_T24_integrated, resolution = 1)

DefaultAssay(T20_T24_integrated) <- "RNA"
T20_T24_integrated <- SCTransform(T20_T24_integrated, verbose = FALSE) %>% RunPCA() %>% RunUMAP(dims = 1:50, reduction.name = 'umap.rna', reduction.key = 'rnaUMAP_')

T20_peaks <- read.table(
  file = "/atac_peaks.bed",
  col.names = c("chr", "start", "end")
)
T24_peaks <- read.table(
  file = "/atac_peaks.bed",
  col.names = c("chr", "start", "end")
)
tail(T20_peaks)

T20_peaks <- makeGRangesFromDataFrame(T20_peaks)
T24_peaks <- makeGRangesFromDataFrame(T24_peaks)
T20_T24_peaks <- reduce(x = c(T20_peaks, T24_peaks))
peakwidths <- width(T20_T24_peaks)
T20_T24_peaks <- T20_T24_peaks[peakwidths  < 10000 & peakwidths > 20]
T20_T24_peaks

T20fragpath <- "/atac_fragments.tsv.gz"
T24fragpath <- "/atac_fragments.tsv.gz"
T20_fragcounts <- CountFragments(fragments = T20fragpath)
T24_fragcounts <- CountFragments(fragments = T24fragpath)

T20_atac_cells <- T20_fragcounts[T20_fragcounts$frequency_count > 10, "CB"]
T24_atac_cells <- T24_fragcounts[T24_fragcounts$frequency_count > 10, "CB"]
T20_frags <- CreateFragmentObject(path = T20fragpath, cells = T20_atac_cells)
T24_frags <- CreateFragmentObject(path = T24fragpath, cells = T24_atac_cells)

t20.counts <- FeatureMatrix(
  fragments = T20_frags,
  features = T20_T24_peaks,
  cells = T20_atac_cells
)
t24.counts <- FeatureMatrix(
  fragments = T24_frags,
  features = T20_T24_peaks,
  cells = T24_atac_cells
)
annotations <- rtracklayer::import('/Volumes/BigHannah/2023_Hannah_Multiome/new.edited_SL_allfilled_eGFP_forSignac_v6.gtf')

t20_assay <- CreateChromatinAssay(t20.counts, fragments = T20_frags, genome = 'GCA_010993605.1', annotation = annotations)
t20_object <- CreateSeuratObject(t20_assay, assay = "ATAC", project = "T20_CNC")
t24_assay <- CreateChromatinAssay(t24.counts, fragments = T24_frags, genome = 'GCA_010993605.1', annotation = annotations)
t24_object <- CreateSeuratObject(t24_assay, assay = "ATAC", project = "T24_TNC")

t20_object$dataset <- 'T20_CNC'
t24_object$dataset <- 'T24_TNC'

t20_object <- RenameCells(t20_object, add.cell.id = "T20")
t24_object <- RenameCells(t24_object, add.cell.id = "T24")

T20_T24_ATAC_merged <- merge(
  x = t20_object,
  y = list(t24_object)
)

# added to run integration anchors 
t20_object <- FindTopFeatures(t20_object, min.cutoff = 10)
t20_object <- RunTFIDF(t20_object)
t20_object <- RunSVD(t20_object)
t20_object <- RunUMAP(t20_object, reduction = "lsi", dims = 1:30)
t24_object <- FindTopFeatures(t24_object, min.cutoff = 10)
t24_object <- RunTFIDF(t24_object)
t24_object <- RunSVD(t24_object)
t24_object <- RunUMAP(t24_object, reduction = "lsi", dims = 1:30)

DefaultAssay(T20_T24_ATAC_merged) <- "ATAC"
T20_T24_ATAC_merged <- RunTFIDF(T20_T24_ATAC_merged)
T20_T24_ATAC_merged <- FindTopFeatures(T20_T24_ATAC_merged, min.cutoff = 'q0')
T20_T24_ATAC_merged <- RunSVD(T20_T24_ATAC_merged)
T20_T24_ATAC_merged <- RunUMAP(T20_T24_ATAC_merged, reduction = 'lsi', dims = 1:50, reduction.name = "umap.atac", reduction.key = "atacUMAP_")

integration.anchors <- FindIntegrationAnchors(
  object.list = list(t20_object, t24_object),
  anchor.features = rownames(T20_T24_ATAC_merged_processed),
  reduction = "rlsi",
  dims = 1:30
)
# integrate LSI embeddings 
T20_T24_ATAC_integrated <- IntegrateEmbeddings(
  anchorset = integration.anchors,
  reductions = T20_T24_ATAC_merged_processed[["lsi"]],
  new.reduction.name = "integrated_lsi",
  dims.to.integrate = 1:30
)

T20_T24_integrated[['ATAC']] <- T20_T24_ATAC_integrated[['ATAC']]
T20_T24_integrated[['integrated_lsi']] <- T20_T24_ATAC_integrated[['integrated_lsi']]
T20_T24_integrated

I also tried these alternatives, both coming from other folk's github inquiries, but neither were quite successful for me: 

DefaultAssay(T20_T24_RNA_merged) <- "RNA"
T20_T24_integrated <- SCTransform(T20_T24_RNA_merged, verbose = FALSE) %>% RunPCA() %>% RunUMAP(dims = 1:50, reduction.name = 'umap.rna', reduction.key = 'rnaUMAP_')

features <- SelectIntegrationFeatures(T20_T24_integrated, nfeatures = 3000)
T20_T24_integrated <- PrepSCTIntegration(
  T20_T24_integrated,
  anchor.features = features
)
T20_T24_integrated <- lapply(
  X = T20_T24_integrated,
  FUN = RunPCA,
  features = features,
  verbose = FALSE
)
integration_anchors <- FindIntegrationAnchors(
  object.list = T20_T24_integrated,
  normalization.method = "SCT",
  anchor.features = features,
  dims = 1:30,
  reduction = "rpca",
  k.anchor = 20,
)
T20_T24_integrated_test <- IntegrateData(
  anchorset = integration_anchors,
  normalization.method = "SCT",
  new.assay.name = "integratedRNA",
  dims = 1:30
)
#run LSI on new seurat object with integrated RNA assay
DefaultAssay(T20_T24_integrated_test) <- "ATAC"
T20_T24_integrated_test <- RunTFIDF(T20_T24_integrated)
T20_T24_integrated_test <- FindTopFeatures(so, min.cutoff = "q25")
T20_T24_integrated_test <- RunSVD(T20_T24_integrated)
#integrate embeddings and output new object to prevent overwriting integrated RNA
#I can't quite remember why I did this the first time, but it was the work around that I needed
T20_T24_integrated_test_ATAC <- IntegrateEmbeddings(
  anchorset = integration_anchors,
  new.reduction.name = "integratedLSI",
  reductions = T20_T24_integrated_test@reductions$lsi
)
#copy integrated LSI from duplicate seurat object to original object
T20_T24_integrated_test@reductions$integratedLSI <- T20_T24_integrated_test_ATAC@reductions$integratedLSI
T20_T24_RNA_merged.list <- SplitObject(T20_T24_RNA_merged, split.by = "orig.ident")
T20_T24_RNA_merged.list <- lapply(X = T20_T24_RNA_merged.list, FUN = function(x) {
    x <- NormalizeData(x)
    x <- FindVariableFeatures(x, selection.method = "vst", nfeatures = 2000)
})

features <- SelectIntegrationFeatures(object.list = T20_T24_RNA_merged.list)
T20_T24_RNA_merged.anchors <- FindIntegrationAnchors(object.list = T20_T24_RNA_merged.list, anchor.features = features)
T20_T24_RNA_integrated <- IntegrateData(anchorset = T20_T24_RNA_merged.anchors)
T20_T24_RNA_integrated <- RunPCA(T20_T24_RNA_integrated)

# Then do the ATAC processing, but integrate using the same anchors as from the RNA. 
T20_peaks <- read.table(
  file = "/atac_peaks.bed",
  col.names = c("chr", "start", "end")
)
T24_peaks <- read.table(
  file = "/atac_peaks.bed",
  col.names = c("chr", "start", "end")
)

T20_peaks <- makeGRangesFromDataFrame(T20_peaks)
T24_peaks <- makeGRangesFromDataFrame(T24_peaks)

T20_T24_peaks <- reduce(x = c(T20_peaks, T24_peaks))
peakwidths <- width(T20_T24_peaks)
T20_T24_peaks <- T20_T24_peaks[peakwidths  < 10000 & peakwidths > 20]

T20fragpath <- "/atac_fragments.tsv.gz"
T24fragpath <- "/atac_fragments.tsv.gz"
T20_fragcounts <- CountFragments(fragments = T20fragpath)
T24_fragcounts <- CountFragments(fragments = T24fragpath)
T20_frags <- CreateFragmentObject(path = T20fragpath)
T24_frags <- CreateFragmentObject(path = T24fragpath)

t20.counts <- FeatureMatrix(
  fragments = T20_frags,
  features = T20_T24_peaks,
  cells = T20_RNA_cells
)
t24.counts <- FeatureMatrix(
  fragments = T24_frags,
  features = T20_T24_peaks,
  cells = T24_RNA_cells
)

annotations <- rtracklayer::import('/file.gtf')

t20_assay <- CreateChromatinAssay(t20.counts, fragments = T20_frags, genome = 'GCA_010993605.1', annotation = annotations)
t20_object <- CreateSeuratObject(t20_assay, assay = "ATAC", project = "T20_CNC")
t24_assay <- CreateChromatinAssay(t24.counts, fragments = T24_frags, genome = 'GCA_010993605.1', annotation = annotations)
t24_object <- CreateSeuratObject(t24_assay, assay = "ATAC", project = "T24_TNC")

t20_object$dataset <- 'T20_CNC'
t24_object$dataset <- 'T24_TNC'
t20_object <- RenameCells(t20_object, add.cell.id = "T20")
t24_object <- RenameCells(t24_object, add.cell.id = "T24")

T20_T24_ATAC_merged <- merge(
  x = t20_object,
  y = list(t24_object)
)
DefaultAssay(T20_T24_ATAC_merged) <- "ATAC"
T20_T24_ATAC_merged <- RunTFIDF(T20_T24_ATAC_merged)
T20_T24_ATAC_merged <- FindTopFeatures(T20_T24_ATAC_merged, min.cutoff = 'q0')
T20_T24_ATAC_merged <- RunSVD(T20_T24_ATAC_merged)
T20_T24_ATAC_merged <- RunUMAP(T20_T24_ATAC_merged, reduction = 'lsi', dims = 1:50, reduction.name = "umap.atac", reduction.key = "atacUMAP_")

T20_T24_RNA_integrated <- IntegrateEmbeddings(
  anchorset = T20_T24_RNA_merged.anchors,
  query = T20_T24_ATAC_merged,
  reductions = T20_T24_ATAC_merged[["lsi"]],
  new.reduction.name = "integrated_lsi",
  dims.to.integrate = 1:30
)
T20_T24_multi <- FindMultiModalNeighbors(T24_multi, reduction.list = list("pca", "lsi"), dims.list = list(1:50, 2:50))
T24_multi <- RunUMAP(T24_multi, nn.name = "weighted.nn", reduction.name = "wnn.umap", reduction.key = "wnnUMAP_")
T24_multi <- FindClusters(T24_multi, graph.name = "wsnn", algorithm = 3, verbose = FALSE)
Dexiao211 commented 11 months ago

Hi @hannahvanm Thanks for sharing your codes # create fragment objects frags.500 <-CreateFragmentObject( path = "../vignette_data/pbmc500/atac_pbmc_500_nextgem_fragments.tsv.gz", cells = rownames(md.500) ) You can import barcodes (or cells) from the RNA assay into the 'cells' parameter in this line. I believe this can resolve your issue.

hannahvanm commented 11 months ago

@Dexiao211 Thank you so much for that suggestion! I had tried adding cells like that elsewhere, but it didn't work. The rownames didn't work for me, but what did was object_cells = Cells(rna_obj) then use object_cells in createfragment object. So in total, here's the complete code that worked for me, thanks to your help! 

T20rna_counts <- T20$`Gene Expression`
T24rna_counts <- T24$`Gene Expression`

T20_CNC <- CreateSeuratObject(counts = T20rna_counts, project = "T20_CNC", min.cells = 3, min.features = 100)
T24_TNC <- CreateSeuratObject(counts = T24rna_counts, project = "T24_CNC", min.cells = 3, min.features = 100)

T20_CNC[["percent.mt"]] <- PercentageFeatureSet(T20_CNC, pattern = "^mt-")
T24_TNC[["percent.mt"]] <- PercentageFeatureSet(T24_TNC, pattern = "^mt-")
VlnPlot(T20_CNC, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
VlnPlot(T24_TNC, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
T20_CNC <- subset(T20_CNC, subset = nFeature_RNA < 3100 & percent.mt < 3)
T24_TNC <- subset(T24_TNC, subset = nFeature_RNA < 3100 & percent.mt < 3)

T20_CNC <- NormalizeData(T20_CNC)
T20_CNC <- FindVariableFeatures(T20_CNC)
T20_CNC <- ScaleData(T20_CNC)
T20_CNC <- RunPCA(T20_CNC)
T24_TNC <- NormalizeData(T24_TNC)
T24_TNC <- FindVariableFeatures(T24_TNC)
T24_TNC <- ScaleData(T24_TNC)
T24_TNC <- RunPCA(T24_TNC)

# because we're going to use the cells from here to put into the ATAC, should we use the raw or filtered cells. Additionally, should we analyze and scale everything before merging? Or not because we're integrating later? 
head(T20_CNC)
rownames(T20_CNC)
T20_cells <- Cells(T20_CNC)
T24_cells <- Cells(T24_TNC)
head(T20_cells)

T20_CNC <- RenameCells(T20_CNC, add.cell.id = "T20")
T24_TNC <- RenameCells(T24_TNC, add.cell.id = "T24")
head(T20_CNC)

T20_T24_RNA_merged <- merge(T20_CNC, y = T24_TNC, project = "T20_T24_RNA_merged")

T20_T24_RNA_merged <- NormalizeData(T20_T24_RNA_merged)
T20_T24_RNA_merged <- FindVariableFeatures(T20_T24_RNA_merged)
T20_T24_RNA_merged <- ScaleData(T20_T24_RNA_merged)
T20_T24_RNA_merged <- RunPCA(T20_T24_RNA_merged)
ElbowPlot(T20_T24_RNA_merged)
T20_T24_RNA_merged <- FindNeighbors(T20_T24_RNA_merged, dims = 1:20, reduction = "pca")
T20_T24_RNA_merged <- FindClusters(T20_T24_RNA_merged, resolution = 1, cluster.name = "unintegrated_clusters")
T20_T24_RNA_merged <- RunUMAP(T20_T24_RNA_merged, dims = 1:20, reduction = "pca", reduction.name = "umap")

T20_T24_RNA_merged.list <- SplitObject(T20_T24_RNA_merged, split.by = "orig.ident")
T20_T24_RNA_merged.list <- lapply(X = T20_T24_RNA_merged.list, FUN = function(x) {
    x <- NormalizeData(x)
    x <- FindVariableFeatures(x, selection.method = "vst", nfeatures = 2000)
})

features <- SelectIntegrationFeatures(object.list = T20_T24_RNA_merged.list)
T20_T24_RNA_merged.anchors <- FindIntegrationAnchors(object.list = T20_T24_RNA_merged.list, anchor.features = features, reduction = "rpca")
T20_T24_RNA_integrated <- IntegrateData(anchorset = T20_T24_RNA_merged.anchors)

T20_T24_RNA_integrated <- ScaleData(T20_T24_RNA_integrated)
T20_T24_RNA_integrated <- RunPCA(T20_T24_RNA_integrated)

# Then do the ATAC processing, but integrate using the same anchors as from the RNA. 
T20_peaks <- read.table(
  file = "/Volumes/BigHannah/2023_Hannah_Multiome/T20_Hannah_multiome/outs/atac_peaks.bed",
  col.names = c("chr", "start", "end")
)
T24_peaks <- read.table(
  file = "/Volumes/BigHannah/2023_Hannah_Multiome/T24_Hannah_multiome/outs/atac_peaks.bed",
  col.names = c("chr", "start", "end")
)
T20_peaks <- makeGRangesFromDataFrame(T20_peaks)
T24_peaks <- makeGRangesFromDataFrame(T24_peaks)

T20_T24_peaks <- reduce(x = c(T20_peaks, T24_peaks))

peakwidths <- width(T20_T24_peaks)
T20_T24_peaks <- T20_T24_peaks[peakwidths  < 10000 & peakwidths > 20]

T20fragpath <- "/Volumes/BigHannah/2023_Hannah_Multiome/T20_Hannah_multiome/outs/atac_fragments.tsv.gz"
T24fragpath <- "/Volumes/BigHannah/2023_Hannah_Multiome/T24_Hannah_multiome/outs/atac_fragments.tsv.gz"
T20_fragcounts <- CountFragments(fragments = T20fragpath)
T24_fragcounts <- CountFragments(fragments = T24fragpath)
head(T20_fragcounts)

T20_frags <- CreateFragmentObject(path = T20fragpath, cells = T20_cells)
#782 CB
T24_frags <- CreateFragmentObject(path = T24fragpath, cells = T24_cells)
rownames(T20_CNC)
#3280 CB

t20.counts <- FeatureMatrix(
  fragments = T20_frags,
  features = T20_T24_peaks,
  cells = T20_cells
)
t24.counts <- FeatureMatrix(
  fragments = T24_frags,
  features = T20_T24_peaks,
  cells = T24_cells
)
annotations <- rtracklayer::import('/Volumes/BigHannah/2023_Hannah_Multiome/new.edited_SL_allfilled_eGFP_forSignac_v6.gtf')

t20_assay <- CreateChromatinAssay(t20.counts, fragments = T20_frags, genome = 'GCA_010993605.1', annotation = annotations)
t20_object <- CreateSeuratObject(t20_assay, assay = "ATAC", project = "T20_CNC")
t24_assay <- CreateChromatinAssay(t24.counts, fragments = T24_frags, genome = 'GCA_010993605.1', annotation = annotations)
t24_object <- CreateSeuratObject(t24_assay, assay = "ATAC", project = "T24_TNC")

t20_object <- RenameCells(t20_object, add.cell.id = "T20")
t24_object <- RenameCells(t24_object, add.cell.id = "T24")
head(t20_object)

T20_T24_ATAC_merged <- merge(
  x = t20_object,
  y = list(t24_object)
)

DefaultAssay(T20_T24_ATAC_merged) <- "ATAC"
T20_T24_ATAC_merged <- RunTFIDF(T20_T24_ATAC_merged)
T20_T24_ATAC_merged <- FindTopFeatures(T20_T24_ATAC_merged, min.cutoff = 'q0')
T20_T24_ATAC_merged <- RunSVD(T20_T24_ATAC_merged)
T20_T24_ATAC_merged <- RunUMAP(T20_T24_ATAC_merged, reduction = 'lsi', dims = 1:50, reduction.name = "umap.atac", reduction.key = "atacUMAP_")

T20_T24_ATAC_merged.list <- SplitObject(T20_T24_ATAC_merged, split.by = "orig.ident")

integration.anchors <- FindIntegrationAnchors(
  object.list = T20_T24_ATAC_merged.list,
  reduction = "rlsi",
  dims = 1:30
)
#anchor.features = rownames(T20_T24_ATAC_merged)

T20_T24_ATAC_integrated <- IntegrateEmbeddings(
  anchorset = integration.anchors,
  reductions = T20_T24_ATAC_merged[["lsi"]],
  new.reduction.name = "integrated_lsi",
  dims.to.integrate = 1:30
)

T20_T24_RNA_integrated[['ATAC']] <- T20_T24_ATAC_integrated[['ATAC']]
T20_T24_RNA_integrated[['integrated_lsi']] <- T20_T24_ATAC_integrated[['integrated_lsi']]

T20_T24_RNA_integrated <- FindMultiModalNeighbors(T20_T24_RNA_integrated, reduction.list = list("pca", "integrated_lsi"), dims.list = list(1:50, 2:50))
T20_T24_RNA_integrated <- RunUMAP(T20_T24_RNA_integrated, nn.name = "weighted.nn", reduction.name = "wnn.umap", reduction.key = "wnnUMAP_")
T20_T24_RNA_integrated <- FindClusters(T20_T24_RNA_integrated, graph.name = "wsnn", algorithm = 3, verbose = FALSE)
T20_T24_RNA_integrated

p1 <- DimPlot(T20_T24_RNA_integrated, reduction = "pca", label = TRUE, label.size = 2.5, repel = TRUE) + ggtitle("RNA")
p2 <- DimPlot(T20_T24_RNA_integrated, reduction = "integrated_lsi",  label = TRUE, label.size = 2.5, repel = TRUE) + ggtitle("ATAC")
p3 <- DimPlot(T20_T24_RNA_integrated, reduction = "wnn.umap", label = TRUE, label.size = 2.5, repel = TRUE) + ggtitle("WNN")
p3
p1 + p2 + p3 & NoLegend() & theme(plot.title = element_text(hjust = 0.5))
p3 & NoLegend() & theme(plot.title = element_text(hjust = 0.5))

p4 <- DimPlot(T20_T24_RNA_integrated, reduction = "wnn.umap", label = TRUE, label.size = 2.5, repel = TRUE, group.by = "orig.ident") + ggtitle("WNN")