WNN analysis of CITE-seq, RNA + ADT

1. load libraries

2. Load Seurat Object

 load("/home/bioinfo/Cluster_to_Computer_Transfer_files_folder/All_Normal-PBMC_Abnormal-cellLines_T_cells_Merged_Annotated_UMAP_on_Clusters_to_USE.Robj")
 
  All_samples_Merged

An object of class Seurat 
62625 features across 46976 samples within 6 assays 
Active assay: SCT (25901 features, 3000 variable features)
 3 layers present: counts, data, scale.data
 5 other assays present: RNA, ADT, prediction.score.celltype.l1, prediction.score.celltype.l2, prediction.score.celltype.l3
 4 dimensional reductions calculated: pca, umap, integrated_dr, ref.umap

  
  DimPlot(All_samples_Merged,group.by = "cell_line", 
        reduction = "umap",
        label.size = 3,
        repel = T,
        label = T)

DimPlot(All_samples_Merged,
        group.by = "SCT_snn_res.0.9", 
        reduction = "umap",
        label.size = 3,
        repel = T,
        label = T)

cluster_table <- table(Idents(All_samples_Merged))


DimPlot(All_samples_Merged, group.by = "predicted.celltype.l2", 
        reduction = "umap",
        label.size = 3,
        repel = T,
        label = T)

library(clustree)

Loading required package: ggraph

Attaching package: 'ggraph'

The following object is masked from 'package:sp':

    geometry

clustree(All_samples_Merged, prefix = "SCT_snn_res.")

DimPlot(All_samples_Merged, group.by = "predicted.celltype.l1", 
        reduction = "umap",
        label.size = 3,
        repel = T,
        label = T, label.box = T)

DimPlot(All_samples_Merged, group.by = "predicted.celltype.l1", 
        reduction = "umap",
        label.size = 3,
        repel = T,
        label = F)

DimPlot(All_samples_Merged, group.by = "predicted.celltype.l2", 
        reduction = "umap",
        label.size = 3,
        repel = T,
        label = T, label.box = T)

DimPlot(All_samples_Merged, group.by = "predicted.celltype.l2", 
        reduction = "umap",
        label.size = 3,
        repel = T,
        label = F)

DimPlot(All_samples_Merged, group.by = "predicted.celltype.l2", 
        reduction = "umap",
        label.size = 3,
        repel = T,
        label = T, label.box = T)

table(All_samples_Merged$predicted.celltype.l2, All_samples_Merged$SCT_snn_res.0.9)

                   
                       0    1    2    3    4    5    6    7    8    9   10   11   12   13   14   15   16   17   18
  CD4 CTL              0    0    0    0    0    0    2    0    4    0    0    0    0    0    0    0    0    0    0
  CD4 Naive            0    0    0    0    0    0  595    0  100    0    0    1    0    0    0    0    6    0    4
  CD4 Proliferating 5124 5231 3823 2401 1748 2535    0 2547    0 1263 1329   20 1327  598   41  237    0  154    1
  CD4 TCM            978  158  501 1366   21   24 2170   22 1892  266  585 1859   41  111  406   57  235   35   61
  CD4 TEM              0    0    0    8    0    0   19    0   27    0    0   15    0    0    0    0    0    0    0
  CD8 Naive            6    0    2    0    0    0  304    0   61    0   19    0    1    3    2    0   12    0    3
  CD8 TCM              0    0    0   76    0    0  115    0  137   10    0   72    0    0    0    0    2    0    1
  CD8 TEM              0    0    0    0    0    0   27    0  182    8    0    0    0    0    0    0    0    0    0
  cDC2                 9    0  155    0    0   23    0    5    1    0   67    0    9   43    0   22    0    2    0
  dnT                  1    0    0    2    2    0    4    0   33    0    0    0    0    0    0    7   11    0    1
  gdT                  0    0    0    0    0    0    0    0   13    0    0    0    0    0    0    0    0    0    0
  HSPC               167    8  285    1    0  744    1  652    1    1   12    8  401   43    1   26    0    6    0
  ILC                  0    0    0    0    0    0    0    0    1    0    0    0    0    0    0    0    0    0    0
  MAIT                 0    0    0    0    0    0    4    0   50    0    0    0    0    0    0    0    2    0    0
  NK                   0    0    0    0    0    0    0    0   91    0    0    0    0    0    0    0    0    0    1
  NK Proliferating     7   11  162   20 1971    4    0    4    0  644   10    0    0    0    0   31    0    0    0
  Treg                11    0    2    1    0    0   50    0  103    0    0    0    0    1    0   20   12    0    5

3. Data PREPERATION

DefaultAssay(All_samples_Merged) <- 'ADT'

VariableFeatures(All_samples_Merged) <- rownames(All_samples_Merged[["ADT"]])

# we will use all ADT features for dimensional reduction
# we set a dimensional reduction name to avoid overwriting the 
All_samples_Merged <- NormalizeData(All_samples_Merged, normalization.method = 'CLR', margin = 2) %>% 
  ScaleData() %>% RunPCA(reduction.name = 'apca', npcs =28, maxit = 5000)

Normalizing across cells

  |                                                  | 0 % ~calculating  
  |+                                                 | 1 % ~00s          
  |+                                                 | 2 % ~00s          
  |++                                                | 3 % ~00s          
  |++                                                | 4 % ~00s          
  |+++                                               | 5 % ~00s          
  |+++                                               | 6 % ~00s          
  |++++                                              | 7 % ~00s          
  |++++                                              | 8 % ~00s          
  |+++++                                             | 9 % ~00s          
  |+++++                                             | 10% ~00s          
  |++++++                                            | 11% ~00s          
  |++++++                                            | 12% ~00s          
  |+++++++                                           | 13% ~00s          
  |+++++++                                           | 14% ~00s          
  |++++++++                                          | 15% ~00s          
  |++++++++                                          | 16% ~00s          
  |+++++++++                                         | 17% ~00s          
  |+++++++++                                         | 18% ~00s          
  |++++++++++                                        | 19% ~00s          
  |++++++++++                                        | 20% ~00s          
  |+++++++++++                                       | 21% ~00s          
  |+++++++++++                                       | 22% ~00s          
  |++++++++++++                                      | 23% ~00s          
  |++++++++++++                                      | 24% ~00s          
  |+++++++++++++                                     | 25% ~00s          
  |+++++++++++++                                     | 26% ~00s          
  |++++++++++++++                                    | 27% ~00s          
  |++++++++++++++                                    | 28% ~00s          
  |+++++++++++++++                                   | 29% ~00s          
  |+++++++++++++++                                   | 30% ~00s          
  |++++++++++++++++                                  | 31% ~00s          
  |++++++++++++++++                                  | 32% ~00s          
  |+++++++++++++++++                                 | 33% ~00s          
  |+++++++++++++++++                                 | 34% ~00s          
  |++++++++++++++++++                                | 35% ~00s          
  |++++++++++++++++++                                | 36% ~00s          
  |+++++++++++++++++++                               | 37% ~00s          
  |+++++++++++++++++++                               | 38% ~00s          
  |++++++++++++++++++++                              | 39% ~00s          
  |++++++++++++++++++++                              | 40% ~00s          
  |+++++++++++++++++++++                             | 41% ~00s          
  |+++++++++++++++++++++                             | 42% ~00s          
  |++++++++++++++++++++++                            | 43% ~00s          
  |++++++++++++++++++++++                            | 44% ~00s          
  |+++++++++++++++++++++++                           | 45% ~00s          
  |+++++++++++++++++++++++                           | 46% ~00s          
  |++++++++++++++++++++++++                          | 47% ~00s          
  |++++++++++++++++++++++++                          | 48% ~00s          
  |+++++++++++++++++++++++++                         | 49% ~00s          
  |+++++++++++++++++++++++++                         | 50% ~00s          
  |++++++++++++++++++++++++++                        | 51% ~00s          
  |++++++++++++++++++++++++++                        | 52% ~00s          
  |+++++++++++++++++++++++++++                       | 53% ~00s          
  |+++++++++++++++++++++++++++                       | 54% ~00s          
  |++++++++++++++++++++++++++++                      | 55% ~00s          
  |++++++++++++++++++++++++++++                      | 56% ~00s          
  |+++++++++++++++++++++++++++++                     | 57% ~00s          
  |+++++++++++++++++++++++++++++                     | 58% ~00s          
  |++++++++++++++++++++++++++++++                    | 59% ~00s          
  |++++++++++++++++++++++++++++++                    | 60% ~00s          
  |+++++++++++++++++++++++++++++++                   | 61% ~00s          
  |+++++++++++++++++++++++++++++++                   | 62% ~00s          
  |++++++++++++++++++++++++++++++++                  | 63% ~00s          
  |++++++++++++++++++++++++++++++++                  | 64% ~00s          
  |+++++++++++++++++++++++++++++++++                 | 65% ~00s          
  |+++++++++++++++++++++++++++++++++                 | 66% ~00s          
  |++++++++++++++++++++++++++++++++++                | 67% ~00s          
  |++++++++++++++++++++++++++++++++++                | 68% ~00s          
  |+++++++++++++++++++++++++++++++++++               | 69% ~00s          
  |+++++++++++++++++++++++++++++++++++               | 70% ~00s          
  |++++++++++++++++++++++++++++++++++++              | 71% ~00s          
  |++++++++++++++++++++++++++++++++++++              | 72% ~00s          
  |+++++++++++++++++++++++++++++++++++++             | 73% ~00s          
  |+++++++++++++++++++++++++++++++++++++             | 74% ~00s          
  |++++++++++++++++++++++++++++++++++++++            | 75% ~00s          
  |++++++++++++++++++++++++++++++++++++++            | 76% ~00s          
  |+++++++++++++++++++++++++++++++++++++++           | 77% ~00s          
  |+++++++++++++++++++++++++++++++++++++++           | 78% ~00s          
  |++++++++++++++++++++++++++++++++++++++++          | 79% ~00s          
  |++++++++++++++++++++++++++++++++++++++++          | 80% ~00s          
  |+++++++++++++++++++++++++++++++++++++++++         | 81% ~00s          
  |+++++++++++++++++++++++++++++++++++++++++         | 82% ~00s          
  |++++++++++++++++++++++++++++++++++++++++++        | 83% ~00s          
  |++++++++++++++++++++++++++++++++++++++++++        | 84% ~00s          
  |+++++++++++++++++++++++++++++++++++++++++++       | 85% ~00s          
  |+++++++++++++++++++++++++++++++++++++++++++       | 86% ~00s          
  |++++++++++++++++++++++++++++++++++++++++++++      | 87% ~00s          
  |++++++++++++++++++++++++++++++++++++++++++++      | 88% ~00s          
  |+++++++++++++++++++++++++++++++++++++++++++++     | 89% ~00s          
  |+++++++++++++++++++++++++++++++++++++++++++++     | 90% ~00s          
  |++++++++++++++++++++++++++++++++++++++++++++++    | 91% ~00s          
  |++++++++++++++++++++++++++++++++++++++++++++++    | 92% ~00s          
  |+++++++++++++++++++++++++++++++++++++++++++++++   | 93% ~00s          
  |+++++++++++++++++++++++++++++++++++++++++++++++   | 94% ~00s          
  |++++++++++++++++++++++++++++++++++++++++++++++++  | 95% ~00s          
  |++++++++++++++++++++++++++++++++++++++++++++++++  | 96% ~00s          
  |+++++++++++++++++++++++++++++++++++++++++++++++++ | 97% ~00s          
  |+++++++++++++++++++++++++++++++++++++++++++++++++ | 98% ~00s          
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 99% ~00s          
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=00s  

Centering and scaling data matrix

  |                                                                                                                                
  |                                                                                                                          |   0%
  |                                                                                                                                
  |==========================================================================================================================| 100%
Warning: You're computing too large a percentage of total singular values, use a standard svd instead.Warning: Requested number is larger than the number of available items (28). Setting to 28.Warning: Requested number is larger than the number of available items (28). Setting to 28.Warning: Requested number is larger than the number of available items (28). Setting to 28.Warning: Requested number is larger than the number of available items (28). Setting to 28.Warning: Requested number is larger than the number of available items (28). Setting to 28.PC_ 1 
Positive:  CD5, CD274, TCRab, CD7, CD3, CD28, CD45RO, CD26, CD44, CD62L 
       CXCR3, CD127, CD45, CCR8 
Negative:  CD25, CD30, CD2, CD95, CD45RA, CD4, CCR6, CCR4, CXCR4, CD40 
       CD19, CCR7, PD1, CCR10 
PC_ 2 
Positive:  CD45, CD44, CD45RO, CD4, CD3, CD40, CCR4, CD45RA, CD26, CD95 
       CD2, CD5, CD62L, CD28 
Negative:  CD127, CD19, PD1, CCR7, CCR10, CXCR3, CCR8, CCR6, CD7, CD274 
       CD30, CXCR4, CD25, TCRab 
PC_ 3 
Positive:  CD62L, TCRab, CD28, CD3, CD4, CXCR4, CCR7, CD45RA, CD127, CD5 
       PD1, CCR4, CD95, CD2 
Negative:  CD26, CXCR3, CD45RO, CD274, CD40, CD44, CCR10, CCR6, CD7, CD30 
       CD19, CD25, CD45, CCR8 
PC_ 4 
Positive:  CCR4, CD5, CD274, CD4, CD30, CCR6, CD44, CD28, CD25, CD45RA 
       CD19, CCR10, CXCR3, PD1 
Negative:  CD40, CD45, CD2, CD62L, CD7, CD26, CXCR4, CCR7, CD127, CD95 
       CD45RO, CCR8, CD3, TCRab 
PC_ 5 
Positive:  CD45RO, CD28, CCR8, CD26, CCR4, PD1, CD4, CCR10, CD40, CD127 
       CD19, TCRab, CD2, CCR7 
Negative:  CD7, CD45RA, CD5, CD44, CCR6, CXCR3, CXCR4, CD274, CD62L, CD30 
       CD3, CD95, CD45, CD25 
Warning: Key 'PC_' taken, using 'apca_' instead

# Identify multimodal neighbors. These will be stored in the neighbors slot, 
# and can be accessed using bm[['weighted.nn']]
# The WNN graph can be accessed at bm[["wknn"]], 
# and the SNN graph used for clustering at bm[["wsnn"]]
# Cell-specific modality weights can be accessed at bm$RNA.weight

All_samples_Merged <- FindMultiModalNeighbors(
  All_samples_Merged, reduction.list = list("pca", "apca"), 
  dims.list = list(1:20, 1:18), modality.weight.name = "RNA.weight"
)

Calculating cell-specific modality weights
Finding 20 nearest neighbors for each modality.

  |                                                  | 0 % ~calculating  
  |+++++++++++++++++++++++++                         | 50% ~08s          
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=17s  

Calculating kernel bandwidths

  |                                                  | 0 % ~calculating  
  |+++++++++++++++++++++++++                         | 50% ~01s          
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=02s  

Warning: The number of provided modality.weight.name is not equal to the number of modalities. SCT.weight ADT.weight are used to store the modality weightsFinding multimodal neighbors

  |                                                  | 0 % ~calculating  
  |+++++++++++++++++++++++++                         | 50% ~29s          
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=01m 01s

  |                                                  | 0 % ~calculating  
  |+++++++++++++++++++++++++                         | 50% ~01s          
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=02s  

Constructing multimodal KNN graph
Constructing multimodal SNN graph

4. Visualization RNA+ADT

All_samples_Merged <- RunUMAP(All_samples_Merged, nn.name = "weighted.nn", reduction.name = "wnn.umap", reduction.key = "wnnUMAP_")

Warning: The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
This message will be shown once per session17:46:04 UMAP embedding parameters a = 0.9922 b = 1.112
17:46:04 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 20
17:46:06 Initializing from normalized Laplacian + noise (using RSpectra)
17:46:09 Commencing optimization for 200 epochs, with 1470228 positive edges
Using method 'umap'
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
17:46:26 Optimization finished

All_samples_Merged <- FindClusters(All_samples_Merged, graph.name = "wsnn", algorithm = 3, resolution = 0.9, verbose = FALSE)

p1 <- DimPlot(All_samples_Merged, reduction = 'wnn.umap', label = TRUE, repel = TRUE, label.size = 2.5) + NoLegend()
p2 <- DimPlot(All_samples_Merged, reduction = 'wnn.umap', group.by = 'predicted.celltype.l2', label = TRUE, repel = TRUE, label.size = 2.5) + NoLegend()
p1 + p2

DimPlot(All_samples_Merged, reduction = 'wnn.umap', group.by = "cell_line",label = TRUE, repel = TRUE, label.size = 2.5)

DimPlot(All_samples_Merged, reduction = 'wnn.umap', label = TRUE, repel = TRUE, label.size = 2.5)

DimPlot(All_samples_Merged, reduction = 'wnn.umap', group.by = 'predicted.celltype.l2', label = TRUE, repel = TRUE, label.size = 2.5)

All_samples_Merged <- RunUMAP(All_samples_Merged, reduction = 'pca', dims = 1:30, assay = 'RNA', 
              reduction.name = 'rna.umap', reduction.key = 'rnaUMAP_')

17:46:59 UMAP embedding parameters a = 0.9922 b = 1.112
17:46:59 Read 46976 rows and found 30 numeric columns
17:46:59 Using Annoy for neighbor search, n_neighbors = 30
17:46:59 Building Annoy index with metric = cosine, n_trees = 50
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
17:47:03 Writing NN index file to temp file /tmp/Rtmp5U7myn/filecd834167b5c1
17:47:03 Searching Annoy index using 1 thread, search_k = 3000
17:47:13 Annoy recall = 100%
17:47:14 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
17:47:16 Initializing from normalized Laplacian + noise (using RSpectra)
17:47:23 Commencing optimization for 200 epochs, with 2014172 positive edges
Using method 'umap'
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
17:47:39 Optimization finished

All_samples_Merged <- RunUMAP(All_samples_Merged, reduction = 'apca', dims = 1:18, assay = 'ADT', 
              reduction.name = 'adt.umap', reduction.key = 'adtUMAP_')

17:47:40 UMAP embedding parameters a = 0.9922 b = 1.112
17:47:40 Read 46976 rows and found 18 numeric columns
17:47:40 Using Annoy for neighbor search, n_neighbors = 30
17:47:40 Building Annoy index with metric = cosine, n_trees = 50
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
17:47:43 Writing NN index file to temp file /tmp/Rtmp5U7myn/filecd831fbd05d3
17:47:43 Searching Annoy index using 1 thread, search_k = 3000
17:47:55 Annoy recall = 100%
17:47:56 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
17:47:58 Initializing from normalized Laplacian + noise (using RSpectra)
17:47:59 Commencing optimization for 200 epochs, with 2051920 positive edges
Using method 'umap'
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
17:48:15 Optimization finished

p3 <- DimPlot(All_samples_Merged, reduction = 'rna.umap', group.by = 'predicted.celltype.l2', label = TRUE, 
              repel = TRUE, label.size = 2.5) + NoLegend()
p4 <- DimPlot(All_samples_Merged, reduction = 'adt.umap', group.by = 'predicted.celltype.l2', label = TRUE, 
              repel = TRUE, label.size = 2.5) + NoLegend()
p3 + p4

DimPlot(All_samples_Merged, reduction = 'rna.umap', group.by = "cell_line",label = TRUE, repel = TRUE, label.size = 2.5)

DimPlot(All_samples_Merged, reduction = 'rna.umap', group.by = 'predicted.celltype.l2', label = TRUE, 
              repel = TRUE, label.size = 2.5)

DimPlot(All_samples_Merged, reduction = 'adt.umap', group.by = "cell_line",label = TRUE, repel = TRUE, label.size = 2.5)

DimPlot(All_samples_Merged, reduction = 'adt.umap', group.by = 'predicted.celltype.l2', label = TRUE, 
              repel = TRUE, label.size = 2.5) + NoLegend()

p5 <- FeaturePlot(All_samples_Merged, features = c("adt_CD45RA","adt_CD45RO","adt_CD90"),
                  reduction = 'wnn.umap', max.cutoff = 2, 
                  cols = c("lightgrey","darkgreen"), ncol = 3)

Warning: The following features could not be found adt_CD90Warning: The following requested variables were not found: adt_CD90

p6 <- FeaturePlot(All_samples_Merged, features = c("rna_TRDC","rna_MPO","rna_AVP"), 
                  reduction = 'wnn.umap', max.cutoff = 3, ncol = 3)
p5 / p6

 VlnPlot(All_samples_Merged, features = "SCT.weight", group.by = 'predicted.celltype.l2', sort = TRUE, pt.size = 0.1) +
  NoLegend()

 VlnPlot(All_samples_Merged, features = "ADT.weight", group.by = 'predicted.celltype.l2', sort = TRUE, pt.size = 0.1) +
  NoLegend()

5. Save the Seurat object as an Robj file

save(All_samples_Merged, file = "../0-imp_OBJ_SS/All_samples_Merged_WNN_correct.Robj")