CD45_analysis_sepcificAssign_primary

1 Functions

scAssign_db <- function(tissue) {
  lapply(c("dplyr","Seurat","HGNChelper"), library, character.only = T)
  source("https://raw.githubusercontent.com/IanevskiAleksandr/sc-type/master/R/gene_sets_prepare.R"); source("https://raw.githubusercontent.com/IanevskiAleksandr/sc-type/master/R/sctype_score_.R")
  
  db_ = "https://raw.githubusercontent.com/IanevskiAleksandr/sc-type/master/ScTypeDB_full.xlsx";
  tissue = tissue # e.g. Immune system,Pancreas,Liver,Eye,Kidney,Brain,Lung,Adrenal,Heart,Intestine,Muscle,Placenta,Spleen,Stomach,Thymus 
  
  # prepare gene sets
  gs_list = gene_sets_prepare(db_, tissue)
  return(gs_list)
}

scAssign <- function(seuratObj,gs_list) {
  es.max = sctype_score(scRNAseqData = seuratObj[["RNA"]]@scale.data, scaled = TRUE, 
                      gs = gs_list$gs_positive, gs2 = gs_list$gs_negative) 
  # merge by cluster
  cL_resutls = do.call("rbind", lapply(unique(seuratObj@meta.data$seurat_clusters), function(cl){
      es.max.cl = sort(rowSums(es.max[ ,rownames(seuratObj@meta.data[seuratObj@meta.data$seurat_clusters==cl, ])]), decreasing = !0)
      head(data.frame(cluster = cl, type = names(es.max.cl), scores = es.max.cl, ncells = sum(seuratObj@meta.data$seurat_clusters==cl)), 10)
  }))
  sctype_scores = cL_resutls %>% group_by(cluster) %>% top_n(n = 1, wt = scores)  
  
  # set low-confident (low ScType score) clusters to "unknown"
  sctype_scores$type[as.numeric(as.character(sctype_scores$scores)) < sctype_scores$ncells/4] = "Unknown"
  print(sctype_scores[,1:3])
  
  seuratObj@meta.data$cell_identity = ""
  for(j in unique(sctype_scores$cluster)){
    cl_type = sctype_scores[sctype_scores$cluster==j,]; 
    seuratObj@meta.data$cell_identity[seuratObj@meta.data$seurat_clusters == j] = as.character(cl_type$type[1])
  }
  return(seuratObj)
}

2 Data

acc_immune = LoadH5Seurat(file = "./Data/acc_immune_5KvarGenes.h5seurat")

3 Clustering

acc_immune <- FindClusters(acc_immune, resolution = 5,verbose = F, algorithm = 1)
DimPlot(acc_immune,label = T)

4 Cell type assignment

4.1 cluster scores

gs_list = scAssign_db(tissue = "Immune system")

Warning in checkGeneSymbols(markers_all) : Human gene symbols should be all upper-case except for the ‘orf’ in open reading frames. The case of some letters was corrected. Warning in checkGeneSymbols(markers_all) : x contains non-approved gene symbols Warning in checkGeneSymbols(markers_all) : x contains non-approved gene symbols Warning in checkGeneSymbols(markers_all) : x contains non-approved gene symbols Warning in checkGeneSymbols(markers_all) : x contains non-approved gene symbols Warning in checkGeneSymbols(markers_all) : x contains non-approved gene symbols Warning in checkGeneSymbols(markers_all) : x contains non-approved gene symbols Warning in checkGeneSymbols(markers_all) : x contains non-approved gene symbols Warning in checkGeneSymbols(markers_all) : x contains non-approved gene symbols Warning in checkGeneSymbols(markers_all) : x contains non-approved gene symbols Warning in checkGeneSymbols(markers_all) : x contains non-approved gene symbols Warning in checkGeneSymbols(markers_all) : x contains non-approved gene symbols Warning in checkGeneSymbols(markers_all) : x contains non-approved gene symbols Warning in checkGeneSymbols(markers_all) : x contains non-approved gene symbols Warning in checkGeneSymbols(markers_all) : x contains non-approved gene symbols Warning in checkGeneSymbols(markers_all) : x contains non-approved gene symbols Warning in checkGeneSymbols(markers_all) : x contains non-approved gene symbols Warning in checkGeneSymbols(markers_all) : x contains non-approved gene symbols Warning in checkGeneSymbols(markers_all) : x contains non-approved gene symbols Warning in checkGeneSymbols(markers_all) : x contains non-approved gene symbols Warning in checkGeneSymbols(markers_all) : x contains non-approved gene symbols Warning in checkGeneSymbols(markers_all) : x contains non-approved gene symbols Warning in checkGeneSymbols(markers_all) : x contains non-approved gene symbols Warning in checkGeneSymbols(markers_all) : x contains non-approved gene symbols Warning in checkGeneSymbols(markers_all) : x contains non-approved gene symbols Warning in checkGeneSymbols(markers_all) : x contains non-approved gene symbols Warning in checkGeneSymbols(markers_all) : x contains non-approved gene symbols Warning in checkGeneSymbols(markers_all) : x contains non-approved gene symbols Warning in checkGeneSymbols(markers_all) : x contains non-approved gene symbols Warning in checkGeneSymbols(markers_all) : x contains non-approved gene symbols Warning in checkGeneSymbols(markers_all) : x contains non-approved gene symbols Warning in checkGeneSymbols(markers_all) : x contains non-approved gene symbols Warning in checkGeneSymbols(markers_all) : x contains non-approved gene symbols Warning in checkGeneSymbols(markers_all) : x contains non-approved gene symbols Warning in checkGeneSymbols(markers_all) : x contains non-approved gene symbols Warning in checkGeneSymbols(markers_all) : x contains non-approved gene symbols Warning in checkGeneSymbols(markers_all) : x contains non-approved gene symbols Warning in checkGeneSymbols(markers_all) : x contains non-approved gene symbols Warning in checkGeneSymbols(markers_all) : x contains non-approved gene symbols Warning in checkGeneSymbols(markers_all) : x contains non-approved gene symbols Warning in checkGeneSymbols(markers_all) : x contains non-approved gene symbols

acc_immune %<>% scAssign(gs_list = gs_list)

acc_immune %<>%  SetIdent(value = "cell_identity")
print_tab(DimPlot(acc_immune,group.by = "cell_identity"),title = "UMAP")

UMAP

print_tab(plt = acc_immune$cell_identity %>% table() %>% as.data.frame(), title = "cell count")

cell count

NA

5 Markers

FeaturePlot(acc_immune, features = c("CD8A","MS4A1", "CD4","CD14","MS4A2"))

6 Markers heatmaps

print_tab(plt = 
            DoHeatmap(acc_immune, features =c("CD8A","MS4A1", "CD4","CD14","MS4A2"), size = 3,slot = "data")
          ,title = "markers")

markers

print_tab(plt = 
            DoHeatmap(subset(acc_immune,subset = cell_identity == "Macrophages"), features
                      =gs_list$gs_positive$Macrophages, size = 3,slot = "data")
          ,title = "Macrophages markers")

Macrophages markers

Warning in DoHeatmap(subset(acc_immune, subset = cell_identity == “Macrophages”), : The following features were omitted as they were not found in the data slot for the RNA assay: PF4

NA

7 Cell type correlation

acc_immune = SetIdent(object = acc_immune,value = "cell_identity")
correlation = cor(AverageExpression(object = acc_immune) %>% as.data.frame())
colnames(correlation) <- gsub(pattern = "RNA.",replacement = "",x = colnames(correlation))
rownames(correlation) <- gsub(pattern = "RNA.",replacement = "",x = rownames(correlation))

pheatmap(mat = correlation,main = "Expression correlation")

8 Antigen presenting machinery

apm_genes = c("HLA-A","HLA-B","HLA-C","B2M","TAP1","TAP2", "TAPBP")
apm_score = FetchData(acc_immune,vars = apm_genes,slot = "data") %>% rowMeans()
acc_immune = AddMetaData(object = acc_immune,metadata = apm_score,col.name = "APM_score")

print_tab(plt = FeaturePlot(acc_immune,features = apm_genes),title = "genes")
print_tab(plt = FeaturePlot(acc_immune,features = "APM_score"),title = "score")

9 Exhaustion markers

exhausted_genes = c("PDCD1","CD244","CD160","CTLA4","HAVCR2")
FeaturePlot(acc_immune,features = exhausted_genes)

10 Immune receptors

receptors = c("CCR3", "CCR4", "CCR10","CXCR2", "CXCR3", "CXCR4", "IL17A")
FeaturePlot(acc_immune,features = receptors)

11 CellphoneDB

acc_cancer_cells_pri = readRDS("/sci/labs/yotamd/lab_share/avishai.wizel/R_projects/ACC_microenv/Data/acc_cancer_no146_primaryonly15k_cancercells.rds")
acc_caf = readRDS("/sci/labs/yotamd/lab_share/ACC/ACC_sc/analysis/acc_tpm_nCount_mito_no146_cafs.rds")
acc_immune_pri = subset(acc_immune, subset = origin == "Primary")

# merge cancer and immune
common_genes = intersect(rownames(acc_cancer_cells_pri),rownames(acc_immune_pri)) %>% intersect(rownames(acc_caf_cells))
acc_cancer_and_cd45 = merge(acc_cancer_cells_pri[common_genes,],acc_immune_pri[common_genes,])
overlapping_cells = colnames(acc_cancer_cells_pri) %>% intersect(colnames(acc_caf)) 
acc_cancer_cd45_caf = merge(acc_cancer_and_cd45[common_genes,],acc_caf_cells[common_genes,!colnames(acc_caf_cells) %in% overlapping_cells] )

  #write metadata

#create lum or myo
lum_over_myo = FetchData(object = acc_cancer_cells_pri,vars = "luminal_over_myo")
lum_over_myo$lum_or_myo = "Unknown"
lum_over_myo$lum_or_myo [lum_over_myo$luminal_over_myo>1]  = "Luminal"
lum_over_myo$lum_or_myo [lum_over_myo$luminal_over_myo<(-1)]  = "Myo"
lum_or_myo = lum_over_myo[,"lum_or_myo",drop = F]
names(lum_or_myo)[1] = "cell_identity"
lum_or_myo = lum_or_myo %>% filter(cell_identity != "Unknown") #remove unknown cells
  
# combine
immune_identity =FetchData(object = acc_immune_pri,vars = "cell_identity")
caf_identity =FetchData(object = acc_caf_cells,vars = "cell.type")
names(caf_identity)[1] = "cell_identity"
all_identity = do.call("rbind", list(lum_or_myo, immune_identity, caf_identity))

#rename and sort columns
all_identity$barcode_sample = rownames(all_identity)
all_identity = all_identity %>% rename(cell_type = cell_identity)
all_identity = all_identity[,c(2,1)]

write.table(x = all_identity,file = "./Data/CellphoneDB/metadata_primary.tsv",row.names =F,sep = "\t")

#write normalized counts
acc_cancer_cd45_caf = acc_cancer_cd45_caf[,rownames(all_identity)]
count_matrix = as.data.frame(acc_cancer_cd45_caf@assays[["RNA"]]@data)
fwrite(count_matrix, file = "./Data/CellphoneDB/counts_primary.txt",sep = "\t",row.names = T)

library(ktplots)
library(reticulate)
acc_cancer_cd45_caf$cell_type = all_identity[,2,drop = F] # add cells identities to Seurat

#read data:
pvals =  py$pvalues
means = py$means

#or:
pvals = read.delim("./Data/CellphoneDB/output_primary/statistical_analysis_pvalues_08_01_2023_16:13:32.txt",check.names = F)
means = read.delim("./Data/CellphoneDB/output_primary/statistical_analysis_means_08_01_2023_16:13:32.txt",check.names = F)

12 significant interactions heatmap

plot_cpdb_heatmap(scdata = acc_cancer_cd45_caf, idents = 'cell_type',pvals =  pvals,main = "Number of significant interactions",alpha = 0.05,treeheight_row = 50)

trace(plot_cpdb,edit = T)
untrace(plot_cpdb)
plot_cpdb_with_col_clustering = plot_cpdb
to_insert =quote(
   if (ncol(means_mat) > 2) {
              d <- dist(as.data.frame(t(means_mat)))
              h <- hclust(d)
              means_mat <- means_mat[, h$order, drop = FALSE]
              pvals_mat <- pvals_mat[, h$order, drop = FALSE]
              plot(h, main = cell_type2)
  }
)


body(plot_cpdb_with_col_clustering) <- body(plot_cpdb_with_col_clustering) %>% as.list %>% append(to_insert,after =  41) %>% as.call
trace(plot_cpdb_with_col_clustering,edit = T)

uniq_interactions <- function(cell_type1, cell_type2_a ,cell_type2_b, gene.family = NULL,genes = NULL) {
  require(purrr)
  a = plot_cpdb(cell_type1 = cell_type1, cell_type2 =cell_type2_a, scdata = acc_cancer_cd45_caf,
                      idents = 'cell_type', means = means, pvals = pvals,
                      gene.family = gene.family,return_table = T, p.adjust.method = "fdr",keep_significant_only = T, genes = genes)
  
  b = plot_cpdb(cell_type1 = cell_type1, cell_type2 = cell_type2_b, scdata = acc_cancer_cd45_caf,
                      idents = 'cell_type', means = means, pvals = pvals,
                      gene.family = gene.family,return_table = T,p.adjust.method = "fdr",keep_significant_only = T, genes = genes)
  
  
  all = list()
  i=0
  nelements = a$Var2 %>% unique() %>% length()
  a_couples = a$Var2 %>% unique()%>% as.vector()  
  b_couples = b$Var2 %>% unique()%>% as.vector()  
  for (i in 1:nelements) {
    sig_couples_a = a %>% filter(Var2 == a_couples[i]) %>% filter(pvals_adj <= 0.05) %>%  pull(Var1) %>% as.vector()  
    sig_couples_b = b %>% filter(Var2 == b_couples[i]) %>% filter(pvals_adj <= 0.05) %>%  pull(Var1) %>% as.vector()  
    only_in_a = sig_couples_a[! sig_couples_a %in% sig_couples_b]
    all[[a_couples[i]]] = only_in_a 
    i = i+1
  }
  all = t(map_dfr(all, ~as_tibble(t(.)))) %>% as.data.frame() %>%  set_names(a_couples)
  return(all)
}

13 All interactions

acc_immune <- StashIdent(object = acc_immune, save.name = "cell_identity")

With Seurat 3.X, stashing identity classes can be accomplished with the following:
acc_immune[["cell_identity"]] <- Idents(object = acc_immune)

SeuratDisk::SaveH5Seurat(object = acc_immune,filename = "./Data/acc_immune_5KvarGenes_V2_cellIdentity.h5seurat")

Creating h5Seurat file for version 3.1.5.9900
Adding counts for RNA
Adding data for RNA
Adding scale.data for RNA
Adding variable features for RNA
Adding feature-level metadata for RNA
Adding cell embeddings for pca
Adding loadings for pca
No projected loadings for pca
Adding standard deviations for pca
No JackStraw data for pca
Adding cell embeddings for umap
No loadings for umap
No projected loadings for umap
No standard deviations for umap
No JackStraw data for umap

14 Coinhibitory interactions

print_tab(plot_cpdb(cell_type1 = 'CD|B|FC', cell_type2 = 'Luminal', scdata = acc_cancer_cd45_caf,
          idents = 'cell_type', means = means, pvals = pvals,
          gene.family = 'coinhibitory',return_table = F,max_size = 4,p.adjust.method = "fdr",keep_significant_only = F,cluster_rows = F)+
  ggtitle("coinhibitory Luminal"),title = "Luminal")

print_tab(plot_cpdb(cell_type1 = 'CD|B|FC', cell_type2 = 'Myo', scdata = acc_cancer_cd45_caf,
          idents = 'cell_type', means = means, pvals = pvals,
          gene.family = 'coinhibitory',return_table = F,max_size = 4,p.adjust.method = "fdr",keep_significant_only = F,cluster_rows = F)+
  ggtitle("coinhibitory Myo"),title = "Myo")

print_tab(plt = 
            uniq_interactions(cell_type1 = 'CD|B|FC',cell_type2_a  = 'Luminal',cell_type2_b = "Myo",gene.family = "coinhibitory")
,title = "unique in luminal")


print_tab(plt = 
            uniq_interactions(cell_type1 = 'CD|B|FC',cell_type2_a  = 'Myo',cell_type2_b = "Luminal",gene.family = "coinhibitory")
,title = "unique in myo")

print_tab(plt = 
            uniq_interactions(cell_type1 = 'CD|B|FC',cell_type2_a  = 'CAF',cell_type2_b = "Myo|Luminal",gene.family = "coinhibitory")
,title = "unique in CAF")

15 Chemokines interactions

print_tab(
  plot_cpdb(cell_type1 = 'CD|B|FC', cell_type2 = 'Luminal', scdata = acc_cancer_cd45_caf,
          idents = 'cell_type', means = means, pvals = pvals,
          gene.family = 'chemokines',return_table = F,max_size = 4,p.adjust.method = "fdr",keep_significant_only = F,cluster_rows = F)+
  ggtitle("chemokines Luminal"),title = "Luminal")

print_tab(
  plot_cpdb(cell_type1 = 'CD|B|FC', cell_type2 = 'Myo', scdata = acc_cancer_cd45_caf,
          idents = 'cell_type', means = means, pvals = pvals,
          gene.family = 'chemokines',return_table = F,max_size = 4,p.adjust.method = "fdr",keep_significant_only = F,cluster_rows = F)+
  ggtitle("chemokines Myo"),title = "Myo")


print_tab(plt = 
            uniq_interactions(cell_type1 = 'CD|B|FC',cell_type2_a  = 'Luminal',cell_type2_b = "Myo",gene.family = "chemokines")
,title = "unique in luminal")


print_tab(plt = 
            uniq_interactions(cell_type1 = 'CD|B|FC',cell_type2_a  = 'Myo',cell_type2_b = "Luminal",gene.family = "chemokines")
,title = "unique in myo")

print_tab(plt = 
            uniq_interactions(cell_type1 = 'CD|B|FC',cell_type2_a  = 'CAF',cell_type2_b = "Myo|Luminal",gene.family = "chemokines")
,title = "unique in CAF")

16 Chemokine ligands

genes = c("CXCL1\\D", "CXCL2\\D","CXCL3\\D","CXCL17","C3","CXCL14")
print_tab(plot_cpdb(cell_type1 = 'CD|B|FC', cell_type2 = 'Myo', scdata = acc_cancer_cd45_caf,
    idents = 'cell_type', means = means, pvals = pvals,
 genes = genes,return_table = F,max_size = 4,p.adjust.method = "fdr" ,keep_significant_only = F) 
 ,title = "Myo")

print_tab(plot_cpdb(cell_type1 = 'CD|B|FC', cell_type2 = 'Luminal', scdata = acc_cancer_cd45_caf,
    idents = 'cell_type', means = means, pvals = pvals,
  genes = genes,return_table = F,max_size = 4,p.adjust.method = "fdr" ,keep_significant_only = F) 
,title = "Luminal")


print_tab(plt = 
            uniq_interactions(cell_type1 = 'CD|B|FC',cell_type2_a  = 'Luminal',cell_type2_b = "Myo", genes = genes)
,title = "unique in luminal")


print_tab(plt = 
            uniq_interactions(cell_type1 = 'CD|B|FC',cell_type2_a  = 'Myo',cell_type2_b = "Luminal", genes = genes)
,title = "unique in myo")

print_tab(plt = 
            uniq_interactions(cell_type1 = 'CD|B|FC',cell_type2_a  = 'CAF',cell_type2_b = "Myo|Luminal", genes = genes)
,title = "unique in CAF")

17 CCL22 and CCL28

print_tab(plot_cpdb(cell_type1 = 'CD|B|FC', cell_type2 = 'Myo', scdata = acc_cancer_cd45_caf,
    idents = 'cell_type', means = means, pvals = pvals,
 genes = c("CCL22", "CCL28" ),return_table = F,max_size = 4,p.adjust.method = "fdr" ,keep_significant_only = F),title = "Myo")

print_tab(plot_cpdb(cell_type1 = 'CD|B|FC', cell_type2 = 'Luminal', scdata = acc_cancer_cd45_caf,
    idents = 'cell_type', means = means, pvals = pvals,
  genes = c("CCL22", "CCL28" ),return_table = F,max_size = 6,p.adjust.method = "fdr" ,keep_significant_only = F),title = "Luminal")

plot_cpdb(cell_type1 = 'CD|B|FC', cell_type2 = 'Myo', scdata = acc_cancer_cd45_caf,
    idents = 'cell_type', means = means, pvals = pvals,
 genes = c("JAG", "MYB" ),return_table = F,max_size = 4,p.adjust.method = "fdr" ,keep_significant_only = F) 

plot_cpdb(cell_type1 = 'CD|B|FC', cell_type2 = 'Luminal', scdata = acc_cancer_cd45_caf,
    idents = 'cell_type', means = means, pvals = pvals,
  genes = c("JAG", "MYB" , "NOTCH","HES1","HEY"),return_table = F,max_size = 4,p.adjust.method = "fdr" ,keep_significant_only = F)