WNN analysis of CITE-seq, RNA + ADT-Part2

1. load libraries

2. Load Seurat Object

 load("../5-SS_ScRNA_Data_Analysis/4-ScSS_MyAnalysis_on_SS/0-Important_R_OBJ/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)

Le chargement a nécessité le package : ggraph

Attachement du package : 'ggraph'

L'objet suivant est masqué depuis '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

#Extract the mean protein expression for each cell line:

# Assuming 'cell_line' is a metadata column in your Seurat object
cell_lines <- unique(All_samples_Merged$cell_line)



# Calculate mean expression for each protein in each cell line
mean_expression <- sapply(cell_lines, function(cl) {
  cells <- WhichCells(All_samples_Merged, expression = cell_line == cl)
  rowMeans(LayerData(All_samples_Merged, assay = "ADT", layer = "data")[, cells])
})


#Calculate distances between cell lines based on protein expression:

# Calculate Euclidean distances
dist_matrix <- dist(t(mean_expression))

Plus d’une classe "dist" est trouvée en cache : Utilisation de la première, depuis l’espace de noms 'spam'
Aussi défini par 'BiocGenerics'

# Convert distance matrix to a phylogenetic tree
tree <- nj(dist_matrix)

#Visualize the tree:

# Plot the tree
plot(tree, main = "Cell Line Similarity Tree Based on Protein Expression")

#Optionally, you can create a heatmap of protein expression:

# Create a heatmap
pheatmap(mean_expression, 
         main = "Protein Expression Heatmap",
         scale = "row",
         cluster_rows = TRUE, 
         cluster_cols = TRUE,
         show_rownames = TRUE,
         show_colnames = TRUE)

NA
NA
NA
NA
NA

DefaultAssay(All_samples_Merged) <- 'ADT'

# ====================================
# Phylogenetic Tree Construction and Visualization
# ====================================

# Load necessary libraries
library(ape)
library(phangorn)
library(plotly)
library(phytools)

# ====================================
# Step 1: Prepare Mean Expression Data
# ====================================

# Assuming 'mean_expression' is a matrix with proteins as rows and cell lines as columns
# Transpose mean_expression to ensure rows are cell lines and columns are proteins
mean_expression_cell_lines <- t(mean_expression)

# Step 2: Compute the Distance Matrix
# ====================================

# Calculate the Euclidean distance matrix between the cell lines
dist_matrix <- dist(mean_expression_cell_lines, method = "euclidean")

Plus d’une classe "dist" est trouvée en cache : Utilisation de la première, depuis l’espace de noms 'spam'
Aussi défini par 'BiocGenerics'

# ====================================
# Step 3: Construct the Phylogenetic Tree
# ====================================

# Create a hierarchical clustering object using UPGMA method
hc <- hclust(dist_matrix, method = "average")

# Convert the hierarchical clustering object into a phylogenetic tree
phylo_tree <- as.phylo(hc)

# ====================================
# Step 4: Visualize the Phylogenetic Tree Using Different Methods
# ====================================

# Method 1: Basic Tree Plot Using `ape`
plot(phylo_tree, main = "Phylogenetic Tree of Sézary Syndrome Cell Lines (ape)", 
     type = "phylo", no.margin = TRUE)
tiplabels(cex = 0.7, adj = 0.5)  # Adjust tip label size and position

library(plotly)
library(phytools)

# Generate the phylogenetic tree (this should be done beforehand)
# phylo_tree <- ...  # your phylogenetic tree construction code here

# Prepare data for plotting
edge_data <- as.data.frame(phylo_tree$edge)
edge_data$length <- phylo_tree$edge.length

# Create a scatter plot for the phylogenetic tree
plot_ly(data = edge_data, 
        x = ~length, 
        y = ~V1,  # Accessing the first column of the edge_data
        type = 'scatter', mode = 'lines+text', 
        text = c(phylo_tree$tip.label, rep("", nrow(edge_data) - length(phylo_tree$tip.label))), 
        textposition = 'top') %>%
  layout(title = "Phylogenetic Tree of Sézary Syndrome Cell Lines (plotly)",
         xaxis = list(title = "Branch Length"),
         yaxis = list(title = "Nodes"))

# Method 4: Using `phytools`
if (requireNamespace("phytools", quietly = TRUE)) {
  library(phytools)
  plot(phylo_tree, main = "Phylogenetic Tree of Sézary Syndrome Cell Lines (phytools)", 
       edge.width = 2, show.tip.label = TRUE, type = "unrooted")
} else {
  message("phytools package is not installed. Please install it for this visualization.")
}

# ====================================
# Step 5: Save the Tree in Newick Format (Optional)
# ====================================

# Save the tree as a Newick file for further analysis or visualization
write.tree(phylo_tree, file = "phylogenetic_tree_cell_lines.newick")

# ====================================
# End of Script
# ====================================

DefaultAssay(All_samples_Merged) <- 'ADT'

# Get cluster information
cluster_info <- All_samples_Merged$SCT_snn_res.0.9
print(table(cluster_info))  # This will show the distribution of cells across clusters

cluster_info
   0    1    2    3    4    5    6    7    8    9   10   11   12   13   14   15   16   17   18 
6303 5408 4930 3875 3742 3330 3291 3230 2696 2192 2022 1975 1779  799  450  400  280  197   77 

# Get unique cluster IDs
clusters <- sort(unique(cluster_info))
print(clusters)  # This will show you what cluster IDs are actually present

 [1] 0  1  2  3  4  5  6  7  8  9  10 11 12 13 14 15 16 17 18
Levels: 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18

# Calculate mean expression for each protein in each cluster
mean_expression_clusters <- sapply(clusters, function(cl) {
  cells <- names(cluster_info)[cluster_info == cl]
  rowMeans(LayerData(All_samples_Merged, assay = "ADT", layer = "data")[, cells, drop = FALSE])
})

# Check the dimensions of the resulting matrix
print(dim(mean_expression_clusters))

[1] 28 19

# Create cluster labels starting from 0
cluster_labels <- paste0("Cluster ", seq(0, length(clusters) - 1))

# Set column names
colnames(mean_expression_clusters) <- cluster_labels

# Create the heatmap
pheatmap(mean_expression_clusters, 
         main = "Protein Expression Heatmap by Cluster",
         scale = "row",
         cluster_rows = TRUE, 
         cluster_cols = TRUE,
         show_rownames = TRUE,
         show_colnames = TRUE,
         fontsize_col = 8,
         angle_col = 45)

# For the tree visualization
dist_matrix_clusters <- dist(t(mean_expression_clusters))

Plus d’une classe "dist" est trouvée en cache : Utilisation de la première, depuis l’espace de noms 'spam'
Aussi défini par 'BiocGenerics'

tree_clusters <- nj(dist_matrix_clusters)

plot(tree_clusters, main = "Cluster Similarity Tree Based on Protein Expression")

NA
NA
NA

4. Data PREPERATION

DefaultAssay(All_samples_Merged) <- 'ADT'


# 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% ~15s          
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=29s  

Calculating kernel bandwidths

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

Avis : 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% ~51s          
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=01m 47s

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

Constructing multimodal KNN graph
Constructing multimodal SNN graph

5. Visualization RNA+ADT

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', group.by = "wsnn_res.0.5",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)

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 = "seurat_clusters",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 = "seurat_clusters",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()

FeaturePlot(All_samples_Merged, features = c("adt_CD45RA","adt_CD45RO","adt_CD5","adt_CD274", "adt_CD95"),
                  reduction = 'wnn.umap', max.cutoff = 2, 
                  cols = c("lightgrey","darkgreen"), ncol = 3)

FeaturePlot(All_samples_Merged, features = c("adt_TCRab", "adt_CD7", "adt_CD3", "adt_CD28"),
                  reduction = 'wnn.umap', max.cutoff = 2, 
                  cols = c("lightgrey","darkgreen"), ncol = 3)

            
 FeaturePlot(All_samples_Merged, features = c("adt_CD26", "adt_CD44", "adt_CD62L","adt_CXCR3", "adt_CD127", "adt_CD45"),
                  reduction = 'wnn.umap', max.cutoff = 2, 
                  cols = c("lightgrey","darkgreen"), ncol = 3)     

 
 FeaturePlot(All_samples_Merged, features = c("adt_CCR6","adt_CCR7","adt_CCR8","adt_CCR10"),
                  reduction = 'wnn.umap', max.cutoff = 2, 
                  cols = c("lightgrey","darkgreen"), ncol = 3) 

 
 FeaturePlot(All_samples_Merged, features = c("adt_CD30","adt_CD40","adt_CCR4","adt_CD4" ,"adt_CD25"),
                  reduction = 'wnn.umap', max.cutoff = 2, 
                  cols = c("lightgrey","darkgreen"), ncol = 3) 

 
 FeaturePlot(All_samples_Merged, features = c("adt_PD1","adt_CD62L","adt_CD95","adt_TCRab", "adt_CXCR4"),
                  reduction = 'wnn.umap', max.cutoff = 2, 
                  cols = c("lightgrey","darkgreen"), ncol = 3)  

 FeaturePlot(All_samples_Merged, features = c("adt_CD2","adt_CD28","adt_CD127","adt_CD19"),
                  reduction = 'wnn.umap', max.cutoff = 2, 
                  cols = c("lightgrey","darkgreen"), ncol = 3)  

 
 
 FeaturePlot(All_samples_Merged, features = c("rna_TP53","rna_CARD11","rna_ARID1A","rna_FAS","rna_CCR4"), 
                  reduction = 'wnn.umap', max.cutoff = 3, ncol = 3)

 
 FeaturePlot(All_samples_Merged, features = c("rna_PLS3","rna_STAT4","rna_GATA3","rna_TRAIL","rna_CD1D","rna_RHOA","rna_TNFRSF1B"),
                  reduction = 'wnn.umap', max.cutoff = 2, 
                  cols = c("lightgrey","darkgreen"), ncol = 3)  

Avis : The following features could not be found rna_TRAILAvis : The following requested variables were not found: rna_TRAIL

 
 FeaturePlot(All_samples_Merged, features = c("rna_KIR3DL2","rna_NKp46","rna_IL2RA","rna_TOX","rna_STAT5A"),
                  reduction = 'wnn.umap', max.cutoff = 2, 
                  cols = c("lightgrey","darkgreen"), ncol = 3)  

Avis : The following features could not be found rna_NKp46Avis : The following requested variables were not found: rna_NKp46

 
 FeaturePlot(All_samples_Merged, features = c("rna_MYC","rna_MNT","rna_EPHA4","rna_DNM3","rna_TWIST1"),
                  reduction = 'wnn.umap', max.cutoff = 2, 
                  cols = c("lightgrey","darkgreen"), ncol = 3)  

 
 FeaturePlot(All_samples_Merged, features = c("rna_TRAF2","rna_SELL","rna_miR21","rna_FCL3","rna_PDCD1","rna_CXCL13"),
                  reduction = 'wnn.umap', max.cutoff = 2, 
                  cols = c("lightgrey","darkgreen"), ncol = 3)  

Avis : The following features could not be found rna_miR21, rna_FCL3Avis : The following requested variables were not found: rna_miR21, rna_FCL3

NA
NA
NA
NA
NA

 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()

 
 table(All_samples_Merged$predicted.celltype.l2, All_samples_Merged$wsnn_res.0.5)

                   
                       0    1   10   11   12   13   14   15   16   17   18   19    2    3    4    5    6    7    8    9
  CD4 CTL              0    0    0    0    0    0    1    0    0    0    0    0    0    0    0    0    5    0    0    0
  CD4 Naive            0    0    0    0    0    0   21    0    0    0    0    4    0    0    0    0    1  679    1    0
  CD4 Proliferating 5025 3002  922 1302 1420  581    0  283   31  144   88    0 5242 3797 2573 2600    0    0   69 1300
  CD4 TCM            959  286 1021   46  179  114  550   84  428   40   17   61  132  477   27   25 2144 1605 1991  602
  CD4 TEM              0    0    0    0    1    0    1    0    0    0    0    0    0    0    0    0   45    0   22    0
  CD8 Naive            6    0    0    1    0    3   41    0    2    0    0    3    0    2    0    0   20  316    0   19
  CD8 TCM              0   10   53    0    2    0    9    1    0    0    0    1    0    0    0    0  232   13   92    0
  CD8 TEM              0    8    0    0    0    0    0    0    0    0    0    0    0    0    0    0  209    0    0    0
  cDC2                 9    0    0    9    0   43    0   25    0    2    0    0    0  147   27    5    0    0    0   69
  dnT                  0    0    0    0    3    0   33    7    0    0    0    1    0    1    0    0   14    2    0    0
  gdT                  0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0   13    0    0    0
  HSPC               158    1    0  399    0   36    2   37    3    5    3    0    9  285  755  653    0    1    0   10
  ILC                  0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    1    0    0    0
  MAIT                 0    0    0    0    0    0    2    0    0    0    0    0    0    0    0    0   54    0    0    0
  NK                   0    0    0    0    0    0    0    0    0    0    0    1    0    0    0    0   91    0    0    0
  NK Proliferating     7 2611    8    0   11    0    0   33    0    0    0    0   11  166    4    5    0    0    0    8
  Treg                 3    0    0    0    1    1   29   30    0    0    0    4    0    1    0    0  102   34    0    0

 # library(clustree)
 # clustree(All_samples_Merged, prefix = "wsnn_res.")

6. Save the Seurat object as an Robj file

save(All_samples_Merged, file = "../5-SS_ScRNA_Data_Analysis/4-ScSS_MyAnalysis_on_SS/0-Important_R_OBJ/All_samples_Merged_WNN_correct_on_HPC.Robj")