Myo_lum_score

Parameters:

suffix = "top_25"
threshold = 25
data_to_read = "acc1_cancer_cells_15KnCount_V3.RDS"

• suffix = top_25
• data_to_read = acc1_cancer_cells_15KnCount_V3.RDS
• threshold = 25

functions

source("./.Rprofile")

Welcome to HMSC project.

source_from_github(repositoy = "DEG_functions",version = "0.2.1")

ℹ SHA-1 hash of file is 4aae7e1f13306aa81aa2d8e76dec76035a4dac8f

top_correlated <- function(genes, threshold) {
  markers_expression = FetchData(object = acc1_cancer_cells,vars = genes,slot = "data")
  markers_average = rowMeans(markers_expression) %>% as.data.frame() %>% rename("average" = 1)
  cor_mat = cor(expression %>% t(), markers_average)%>% as.data.frame() 
  cor_mat = cor_mat[complete.cases(cor_mat),,drop=F]  %>% as.data.frame %>%  rename("corr" = 1) #remove rows with NA in at least one col
  if (threshold<1){
    top_genes =   cor_mat %>% as.data.frame %>% select(1) %>% dplyr::filter(.> threshold) %>% rownames()
  }else{
    top_genes =   cor_mat %>%  top_n(threshold,corr) %>% rownames()
  }
  return(top_genes)
}

top_genes_cor_heatmap <- function(top_genes) {
  top_expression = expression %>% dplyr::filter(rownames(expression) %in% top_genes)
colors <- c(seq(-1,1,by=0.01))
my_palette <- c("blue",colorRampPalette(colors = c("blue", "white", "red"))
                                                   (n = length(colors)-3), "red")
pht = pheatmap(mat = cor(top_expression %>% t(), top_expression %>% t()),color = my_palette, breaks = colors)
return(pht)
}

 
enriched_score_umap <- function(enrich_res, genes,col) {
  rownames(enrich_res) = enrich_res$pathway_name
  enriched_genes = enrich_res[col,"geneID"] %>% strsplit(split = "/") %>% .[[1]] %>% c(.,genes) #add original markers
  enriched_genes_score=apply(acc1_cancer_cells@assays$RNA@data[enriched_genes,],2,mean)
  acc1_cancer_cells=AddMetaData(acc1_cancer_cells,enriched_genes_score,"enriched_genes_score")
  print(
    FeaturePlot(object = acc1_cancer_cells,features = "enriched_genes_score")
  )
    return(enriched_genes_score)
}

Data

acc1_cancer_cells = readRDS("./Data/" %>% paste0(data_to_read))
expression = GetAssayData(object = acc1_cancer_cells,assay = "RNA",slot = "data") %>% as.data.frame()

luminal_pathways = c("CHARAFE_BREAST_CANCER_LUMINAL_VS_BASAL_DN","CHARAFE_BREAST_CANCER_LUMINAL_VS_BASAL_UP","CHARAFE_BREAST_CANCER_LUMINAL_VS_MESENCHYMAL_DN","CHARAFE_BREAST_CANCER_LUMINAL_VS_MESENCHYMAL_UP","HUPER_BREAST_BASAL_VS_LUMINAL_DN","LIM_MAMMARY_LUMINAL_PROGENITOR_UP","SMID_BREAST_CANCER_LUMINAL_B_UP" )

# add luminal pathways
msigdb_gene_set = msigdbr(species = "Homo sapiens") 
msigdb_gene_set = msigdb_gene_set %>%as.data.frame() 
luminal_gs = msigdb_gene_set %>%  dplyr::filter(gs_name %in% luminal_pathways)
luminal_gs = luminal_gs %>% dplyr::distinct(gs_name, gene_symbol) %>% as.data.frame()

1.Myo markers

1.a) expression

myo_protein_markers = c("CNN1", "IL12B", "TP63","ACTA2")
FeaturePlot(object = acc1_cancer_cells,features = myo_protein_markers,min.cutoff = 0)

Warning in FeaturePlot(object = acc1_cancer_cells, features = myo_protein_markers,  :
  All cells have the same value (0) of IL12B.

IL12B is 0 in all cells

myo_protein_markers = c("CNN1", "TP63","ACTA2")

1.b) top 30 correlated genes with TP63+ACTA2+CNN1

top_myo  = top_correlated(genes = myo_protein_markers,threshold = 25)

Warning in cor(expression %>% t(), markers_average) :
  the standard deviation is zero

Top correlated genes

print("#genes = " %>% paste(length(top_myo)))

[1] "#genes =  25"

top_myo %>% head(30)

 [1] "COL16A1"     "RP1-39G22.4" "CHI3L1"      "ACTG2"       "CD200"       "MYLK"        "TP63"        "IGFBP7"      "KCNMB1"     
[10] "ADAMTS2"     "CALD1"       "LOXL2"       "TPM2"        "CORO2A"      "CLIC3"       "SNCG"        "ACTA2"       "TAGLN"      
[19] "A2M"         "NGFR"        "MIR7974"     "CNN1"        "PPP1R14A"    "MYL9"        "POM121L9P"  

Correlation

top_genes_cor_heatmap(top_genes = top_myo)

Enrichment analysis

enrich_res = genes_vec_enrichment(genes = top_myo,background = rownames(acc1_cancer_cells),homer = T,title = "myo top enrichment")

enrich_res

UMAP score of enriched genes

lst = enriched_score_umap(enrich_res = enrich_res,genes = myo_protein_markers,col = 1)

myo_score = lst[["score"]]
enriched_genes = lst[["genes"]]

genes in score:

lst[["genes"]]

 [1] "COL16A1" "MYLK"    "CALD1"   "LOXL2"   "TPM2"    "ACTA2"   "TAGLN"   "MYL9"    "CNN1"    "TP63"    "ACTA2"  

CALD1 Lysyl oxidase gene (LOXL2) expression in the stromal reaction to in situ and invasive ductal breast carcinoma TPM2 = tropomyosin ## {.unnumbered}

2. Luminal markers

2.a) expression

lum_protein_markers = c("KIT")
FeaturePlot(object = acc1_cancer_cells,features = lum_protein_markers,min.cutoff = 0)

2 top correlated genes with KIT

top_kit = top_correlated(genes = lum_protein_markers,threshold = 25)

Warning in cor(expression %>% t(), markers_average) :
  the standard deviation is zero

Top correlated genes

print("#genes = " %>% paste(length(top_kit)))

[1] "#genes =  25"

top_kit

 [1] "FBXO44"        "USP48"         "MPC2"          "SLC19A2"       "ATL2"          "B3GNT2"        "AC093162.5"    "MAP2"         
 [9] "GK5"           "KIT"           "GLRB"          "EFNA5"         "PCDHGB7"       "SLC29A1"       "SASH1"         "ALDH3B2"      
[17] "CCND1"         "RP11-254B13.1" "VSIG10"        "NDFIP2"        "SUSD6"         "SPPL2A"        "DYNC1LI2"      "DSC2"         
[25] "INSR"         

Correlation

top_genes_cor_heatmap(top_genes = top_kit)

Enrichment analysis

enrich_res = genes_vec_enrichment(genes = top_kit,background = rownames(expression),homer = T,title = "lum top enrichment",custom_pathways = luminal_gs)

enrich_res

UMAP score of enriched genes

lst = enriched_score_umap(enrich_res = enrich_res,genes = lum_protein_markers,col = 1)

lum_score = lst[["score"]]
enriched_genes = lst[["genes"]]

genes in score:

lst[["genes"]]

[1] "SLC19A2" "GLRB"    "CCND1"   "KIT"    

CCND1 Copy Number Variation in Circulating Tumor DNA from Luminal B Breast Cancer Patients

Correlation between luminal and myo scores

cor(luminal_score,myo_score)

[1] -0.1540176

# res = FindAllMarkers(object = acc1_cancer_cells)

# markers_cluster3 = res[res$cluster == 3 & res$p_val<0.00001,"gene"]

# genes_vec_enrichment(genes = markers_cluster3,background = rownames(acc1_cancer_cells),homer = T,title = "myo top enrichment",custom_pathways = luminal_gs)