Title

Data

acc1_cancer_cells = readRDS("./Data/acc1_cancer_cells_15KnCount_V3.RDS")
all_acc_cancer_cells = readRDS("./Data/acc_cancer_cells_V4.RDS")
acc_all_cells = readRDS("./Data/acc_tpm_nCount_mito_no146_15k_with_ACC1_.RDS")
acc_cancerCells_noACC1 = subset(all_acc_cancer_cells,subset = patient.ident!= "HMSC")

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
luminal_gs = msigdbr(species = "Homo sapiens") %>%as.data.frame() %>% dplyr::filter(gs_name %in% luminal_pathways)%>% dplyr::distinct(gs_name, gene_symbol) %>% as.data.frame()

Parameters

functions

source_from_github(repositoy = "DEG_functions",version = "0.2.24")
source_from_github(repositoy = "HMSC_functions",version = "0.1.13",script_name = "functions.R")
source_from_github(repositoy = "cNMF_functions",version = "0.3.72",script_name = "cnmf_function_Harmony.R")

UMAP

DimPlot(object = acc1_cancer_cells,pt.size = 2)

for (cluster in unique(deg$cluster)) {
  deg_of_cluster = deg[deg$cluster == cluster,]
  enrichment_analysis(deg_of_cluster,background = VariableFeatures(acc1_cancer_cells),fdr_Cutoff = 0.01,ident.1 = paste("Cluster",cluster),ident.2 =  paste("Cluster",cluster),show_by = 1)
}

NA
NA

features UMAP

FeaturePlot(object = acc1_cancer_cells,features = c("TP63","ACTA2","IL12B","CNN1"))

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

Enrichment analysis HMSC vs ACC

all_acc_cancer_cells = SetIdent(all_acc_cancer_cells, value ="patient.ident")
acc_deg <- FindMarkers(all_acc_cancer_cells, ident.1 = "HMSC",logfc.threshold = 1.5,features = VariableFeatures(all_acc_cancer_cells))
enrichment_analysis(acc_deg,background = VariableFeatures(all_acc_cancer_cells),fdr_Cutoff = 0.01,ident.1 = "HMSC",ident.2 = "ACC",show_by = 1)

cell cycle filtering

hallmark_name = "GO_MITOTIC_CELL_CYCLE"
genesets  =getGmt("./Data/h.all.v7.0.symbols.pluscc.gmt")
var_features=acc_cancerCells_noACC1@assays$RNA@var.features
geneIds= genesets[[hallmark_name]]@geneIds
score <- apply(acc_cancerCells_noACC1@assays$RNA@scale.data[intersect(geneIds,var_features),],2,mean)
acc_cancerCells_noACC1=AddMetaData(acc_cancerCells_noACC1,score,hallmark_name)

hallmark_name = "GO_MITOTIC_CELL_CYCLE"
genesets  =getGmt("./Data/h.all.v7.0.symbols.pluscc.gmt")
var_features=acc1_cancer_cells@assays$RNA@var.features
geneIds= genesets[[hallmark_name]]@geneIds
score <- apply(acc1_cancer_cells@assays$RNA@scale.data[intersect(geneIds,var_features),],2,mean)
acc1_cancer_cells=AddMetaData(acc1_cancer_cells,score,hallmark_name)


# hallmark_name = "GO_MITOTIC_CELL_CYCLE"
# genesets  =getGmt("./Data/h.all.v7.0.symbols.pluscc.gmt")
# var_features=acc1_cancer_cells@assays$RNA@var.features
# geneIds= genesets[[hallmark_name]]@geneIds
# score <- apply(acc1_cancer_cells@assays$RNA@scale.data[intersect(geneIds,var_features),],2,mean)
# acc1_cancer_cells=AddMetaData(acc1_cancer_cells,score,hallmark_name)

library(highcharter) 
options(highcharter.theme = hc_theme_smpl(tooltip = list(valueDecimals = 2)))

cc_scores = FetchData(object = acc1_cancer_cells,vars = "GO_MITOTIC_CELL_CYCLE")
 hchart(
  density(cc_scores$GO_MITOTIC_CELL_CYCLE), 
  type = "area", name = "GO_MITOTIC_CELL_CYCLE"
  )

NA

cc_scores = FetchData(object = acc_cancerCells_noACC1,vars = "GO_MITOTIC_CELL_CYCLE")
 hchart(
  density(cc_scores$GO_MITOTIC_CELL_CYCLE), 
  type = "area", name = "GO_MITOTIC_CELL_CYCLE"
  )

NA

Before cc filtering

#filter:
all_acc_cancer_cells_ccFiltered=all_acc_cancer_cells[,all_acc_cancer_cells@meta.data[[hallmark_name]]< 0.16]


min_threshold = min(all_acc_cancer_cells$GO_MITOTIC_CELL_CYCLE)
max_threshold = max(all_acc_cancer_cells$GO_MITOTIC_CELL_CYCLE)

library(viridis)
FeaturePlot(object = all_acc_cancer_cells,features = hallmark_name) + ggtitle("Before cc filtering")  & scale_color_gradientn(colours = plasma(n = 10, direction = -1), limits = c(min_threshold, max_threshold))

Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale.

After cc filtering

FeaturePlot(object = all_acc_cancer_cells_ccFiltered,features = hallmark_name) + ggtitle("After cc filtering") & scale_color_gradientn(colours = plasma(n = 10, direction = -1), limits = c(min_threshold, max_threshold))

Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale.

DimPlot(object = all_acc_cancer_cells_ccFiltered,group.by = "patient.ident")

all_acc_cancer_cells_ccFiltered = SetIdent(object = all_acc_cancer_cells_ccFiltered,value = "orig.ident")
VlnPlot(object = all_acc_cancer_cells_ccFiltered,features = hallmark_name)

DimPlot(object = all_acc_cancer_cells,group.by = "orig.ident")

Enrichment analysis filtered HMSC vs ACC

patient.ident = all_acc_cancer_cells_ccFiltered$patient.ident %>% as.data.frame()
patient.ident[,1] = as.character(patient.ident[,1])
patient.ident[patient.ident[,1] == "ACC1",] = "HMSC"
patient.ident[,1] = as.factor(patient.ident[,1])
all_acc_cancer_cells_ccFiltered = AddMetaData(object = all_acc_cancer_cells_ccFiltered,metadata = patient.ident,col.name = "patient.ident")
all_acc_cancer_cells_ccFiltered = SetIdent(all_acc_cancer_cells_ccFiltered, value ="patient.ident")
acc_deg <- FindMarkers(all_acc_cancer_cells_ccFiltered, ident.1 = "HMSC",logfc.threshold = 1.5,features = VariableFeatures(all_acc_cancer_cells_ccFiltered))
enrichment_analysis(acc_deg,background = VariableFeatures(all_acc_cancer_cells_ccFiltered),fdr_Cutoff = 0.01,ident.1 = "HMSC",ident.2 = "ACC",show_by = 1)

MYB expression

all_acc_cancer_cells = SetIdent(object = all_acc_cancer_cells,value = "patient.ident") #active snn graph
FeaturePlot(object = all_acc_cancer_cells,features = "MYB",label = T)

all_acc_cancer_cells = SetIdent(object = all_acc_cancer_cells,value = "patient.ident") #active snn graph
FeaturePlot(object = all_acc_cancer_cells,features = "MYB",label = T)

CNV

new.cluster.ids <- c("cancer", #0
                     "cancer", #1
                     "CAF", #2
                     "cancer", #3
                     "Endothelial", #4
                     "cancer", #5
                     "cancer", #6
                     "CAF", #7
                     "CAF", #8
                     "CAF", #9
                     "cancer", #10
                     "CAF", #11
                     "cancer", #12
                     "cancer", #13
                     "cancer", #14
                     "cancer", #15
                     "cancer", #16
                     "WBC", #17
                     "CAF" #18
                     )

#rename idents:
acc_all_cells = SetIdent(object = acc_all_cells,value = "RNA_snn_res.1") #active snn graph
names(new.cluster.ids) <- levels(acc_all_cells) #add snn graph levels to new.cluster.ids
acc_all_cells@meta.data[["seurat_clusters"]] = acc_all_cells@meta.data[["RNA_snn_res.1"]]
acc_all_cells = SetIdent(object = acc_all_cells,value = "seurat_clusters")
acc_all_cells <- RenameIdents(acc_all_cells, new.cluster.ids) 

# divide "cancer" into patients:
cell_types = acc_all_cells@active.ident %>% as.data.frame()
cell_types[,1]<- as.character(cell_types[,1])
cell_types = cbind(cell_types,acc_all_cells$patient.ident) %>% setnames(old = names(.), 
         new = c('cell_type','patient'))
cell_types[cell_types$cell_type == "cancer",] = cell_types[cell_types$cell_type == "cancer",2]


# hmsc_rows = (startsWith(x = rownames(cell_types),prefix = "ACC.plate2") | startsWith(x = rownames(cell_types),prefix = "ACC1.")) & cell_types[,1] == "cancer" 
# acc_rows = !(startsWith(x = rownames(cell_types),prefix = "ACC.plate2") | startsWith(x = rownames(cell_types),prefix = "ACC1.")) & cell_types[,1] == "cancer" 
# cell_types[,1][hmsc_rows]  = "HMSC"
# cell_types[,1][acc_rows]  = "ACC"

#add to metadata:
cell_types[,2] = NULL 
cell_types[cell_types$cell_type == "ACC1",] = "HMSC"
acc_all_cells = AddMetaData(object =acc_all_cells ,metadata = cell_types,col.name = "cell.type")

CNV UMAP

library(infercnv)
library(tidyverse)
acc_annotation  = as.data.frame(acc_all_cells@meta.data[,"cell.type",drop = F])
acc_annotation = acc_annotation %>% rownames_to_column("orig.ident") 
acc_annotation = acc_annotation %>% mutate(orig.ident = gsub(x = acc_annotation$orig.ident,pattern = "\\.", replacement = "-") %>% 
  gsub(pattern = "_", replacement = "-"))
  

write.table(acc_annotation, "./Data/inferCNV/acc_annotation.txt", append = FALSE, 
            sep = "\t", dec = ".",row.names = FALSE, col.names = F)

infercnv_obj = CreateInfercnvObject(raw_counts_matrix="./Data/inferCNV/all.4icnv.txt", 
                                    annotations_file="./Data/inferCNV/acc_annotation.txt",
                                    delim="\t",gene_order_file="./Data/inferCNV/gencode_v19_gene_pos.txt"
                                    ,ref_group_names=c("CAF", "Endothelial", "WBC")) #groups of normal cells

infercnv_obj_default = infercnv::run(infercnv_obj, cutoff=1, out_dir='./Data/inferCNV/infercnv_output',
                                     cluster_by_groups=T, plot_steps=FALSE,
                                     denoise=TRUE, HMM=FALSE, no_prelim_plot=TRUE,
                                     png_res=300)
plot_cnv(infercnv_obj_default, output_format = "png",  write_expr_matrix = FALSE,out_dir = "./Data/inferCNV/",png_res   =800,obs_title = "Malignant cells",ref_title = "Normal cells")


cluster.info=FetchData(acc_all_cells,c("ident","orig.ident","UMAP_1","UMAP_2","nCount_RNA","nFeature_RNA","percent.mt","patient.ident","seurat_clusters"))
cluster.info$cell=rownames(cluster.info)

library(limma)
smoothed=apply(infercnv_obj_default@expr.data,2,tricubeMovingAverage, span=0.01)
cnsig=sqrt(apply((smoothed-1)^2,2,mean))
umap=cluster.info
names(cnsig)=umap$cell

acc_all_cells <- AddMetaData(object = acc_all_cells, metadata = cnsig, col.name = "copynumber")
acc_all_cells = SetIdent(object = acc_all_cells,value = "cell.type")
FeaturePlot(acc_all_cells, "copynumber",pt.size = 1, cols = c("lightblue","orange","red","darkred"),label = T,repel = T)

acc_cancer_cells <- AddMetaData(object = acc_cancer_cells, metadata = cnsig, col.name = "copynumber")
acc_cancer_cells = SetIdent(object = acc_cancer_cells,value = "cell.type")
FeaturePlot(acc_cancer_cells, "copynumber",pt.size = 1, cols = c("lightblue","orange","red","darkred"),label = T,repel = T)

CNV plot

## {-}

CNV subtypes

cnv_subtypes = as.data.frame(cutree(infercnv_obj_default@tumor_subclusters[["hc"]][["HMSC"]], k = 2))
names(cnv_subtypes)[1] = "cnv.cluster"
rownames(cnv_subtypes) = rownames(cnv_subtypes) %>% gsub(pattern = "-",replacement = "\\.")
infercnv.observations = data.frame(fread(file = "./Data/inferCNV/infercnv.observations.txt"), row.names=1)

Warning in fread(file = "./Data/inferCNV/infercnv.observations.txt") :
  Detected 1332 column names but the data has 1333 columns (i.e. invalid file). Added 1 extra default column name for the first column which is guessed to be row names or an index. Use setnames() afterwards if this guess is not correct, or fix the file write command that created the file to create a valid file.

names_to_keep = colnames(infercnv.observations) %in% (colnames(acc1_cancer_cells) %>% gsub(pattern = "_",replacement = "\\."))
infercnv.observations = infercnv.observations[,names_to_keep]

rownames(cnv_subtypes) = rownames(cnv_subtypes) %>% gsub(pattern = "2\\.",replacement = "2_")
rownames(cnv_subtypes) = rownames(cnv_subtypes) %>% gsub(pattern = "3\\.",replacement = "3_")
acc1_cancer_cells = AddMetaData(object = acc1_cancer_cells,metadata = cnv_subtypes)

DimPlot(acc1_cancer_cells,group.by = "cnv.cluster",pt.size = 2,cols =colors)

Original score

original_myo_genes = c( "TP63", "TP73", "CAV1", "CDH3", "KRT5", "KRT14", "ACTA2", "TAGLN", "MYLK", "DKK3")
original_lum_genes = c("KIT", "EHF", "ELF5", "KRT7", "CLDN3", "CLDN4", "CD24", "LGALS3", "LCN2", "SLPI" )

calculate_score(dataset = all_acc_cancer_cells,myo_genes = original_myo_genes,lum_genes = original_lum_genes)

correlation of lum score and myo score: -0.45
correlation of lum score and original lum score: 1
correlation of myo score and original myo score: 1

Original score of ACC1

calculate_score(dataset = acc1_cancer_cells,myo_genes = original_myo_genes,lum_genes = original_lum_genes,lum_threshold = 0,myo_threshold = 0)

correlation of lum score and myo score: 0.06
correlation of lum score and original lum score: 1
correlation of myo score and original myo score: 1

0.35 Most correlated score

Myo genes

myo_protein_markers = c("CNN1", "TP63","ACTA2")
top_myo  = top_correlated(dataset = acc1_cancer_cells, genes = myo_protein_markers,threshold = 0.35)
print("Number of genes = " %>% paste(length(top_myo)))

[1] "Number of genes =  20"

message("Names of genes:")

Names of genes:

top_myo %>% head(30)

 [1] "COL16A1"     "RP1-39G22.4" "ACTG2"       "CD200"       "MYLK"        "TP63"        "KCNMB1"     
 [8] "ADAMTS2"     "LOXL2"       "TPM2"        "CLIC3"       "SNCG"        "ACTA2"       "TAGLN"      
[15] "A2M"         "NGFR"        "CNN1"        "PPP1R14A"    "MYL9"        "POM121L9P"  

message("Genes that also apeared in the original score:")

Genes that also apeared in the original score:

base::intersect(top_myo,original_myo_genes) 

[1] "MYLK"  "TP63"  "ACTA2" "TAGLN"

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

myo_enrich_res

Lum genes

lum_protein_markers = c("KIT")
top_lum  = top_correlated(dataset = acc1_cancer_cells, genes = lum_protein_markers,threshold = 0.35,n_vargenes = 5000)

Calculating gene variances
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Calculating feature variances of standardized and clipped values
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

print("Number of genes = " %>% paste(length(top_lum)))

[1] "Number of genes =  5"

message("Names of genes:")

Names of genes:

top_lum %>% head(30)

[1] "B3GNT2"  "GLRB"    "EFNA5"   "ALDH3B2" "CCND1"  

message("Genes that also apeared in the original score:")

Genes that also apeared in the original score:

base::intersect(top_lum,original_lum_genes) 

character(0)

lum_enrich_res = genes_vec_enrichment(genes = top_lum,background = rownames(acc1_cancer_cells),homer = T,title = "lum top enrichment",custom_pathways = luminal_gs)

lum_enrich_res

top correlated score

calculate_score(dataset = acc1_cancer_cells,myo_genes = top_myo,lum_genes = top_lum,lum_threshold = 0,myo_threshold = 0)

correlation of lum score and myo score: -0.05
correlation of lum score and original lum score: 0.58
correlation of myo score and original myo score: 0.77

enriched genes score

rownames(lum_enrich_res) = lum_enrich_res$pathway_name
lum_enriched_genes = lum_enrich_res[1,"geneID"] %>% strsplit(split = "/") %>% .[[1]] %>% c(.,lum_protein_markers) #add original markers

rownames(myo_enrich_res) = myo_enrich_res$pathway_name
myo_enriched_genes = myo_enrich_res[1,"geneID"] %>% strsplit(split = "/") %>% .[[1]] %>% c(.,myo_protein_markers) #add original markers

calculate_score(dataset = acc1_cancer_cells,myo_genes = myo_enriched_genes,lum_genes = lum_enriched_genes,lum_threshold = -2.5,myo_threshold = -2.5)

correlation of lum score and myo score: -0.16
correlation of lum score and original lum score: 0.43
correlation of myo score and original myo score: 0.79

0.2 Most correlated score

myo Genes

n_vargenes = 2000
myo_protein_markers = c("CNN1", "TP63","ACTA2")
top_myo  = top_correlated(dataset = acc1_cancer_cells, genes = myo_protein_markers,threshold = 0.2,n_vargenes = n_vargenes)
print("Number of genes = " %>% paste(length(top_myo)))

[1] "Number of genes =  48"

message("Names of genes:")

Names of genes:

top_myo %>% head(30)

 [1] "PLOD1"         "RP1-39G22.4"   "PDZK1"         "CSRP1"         "CHI3L1"        "HIST3H3"      
 [7] "AC104699.1"    "ACTG2"         "LYG1"          "DALRD3"        "CD200"         "RP11-627C21.1"
[13] "FGFBP2"        "IGFBP7"        "HSPB3"         "RP11-168A11.4" "SPARC"         "CD83"         
[19] "HLA-DOB"       "TBCC"          "NFKBIE"        "MIR3662"       "SMOC2"         "FBXL6"        
[25] "TPM2"          "SNCG"          "ACTA2"         "DKK3"          "MIR7113"       "CTA-797E19.3" 

message("Genes that also apeared in the original score:")

Genes that also apeared in the original score:

base::intersect(top_myo,original_myo_genes) 

[1] "ACTA2" "DKK3"  "TAGLN"

myo_enrich_res = genes_vec_enrichment(genes = top_myo,background = VariableFeatures(acc1_cancer_cells) %>% head(n_vargenes),homer = T,title = "myo top enrichment",custom_pathways = luminal_gs)

myo_enrich_res

Lum Genes

acc1_cancer_cells = FindVariableFeatures(object = acc1_cancer_cells,nfeatures = 15000)

Calculating gene variances
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Calculating feature variances of standardized and clipped values
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

lum_protein_markers = c("KIT")
n_vargenes = 2000
top_lum  = top_correlated(dataset = acc1_cancer_cells, genes = lum_protein_markers,threshold = 0.3,n_vargenes = n_vargenes)
print("Number of genes = " %>% paste(length(top_lum)))

[1] "Number of genes =  13"

message("Names of genes:")

Names of genes:

top_lum %>% head(30)

 [1] "CSDE1"   "RRP9"    "TKT"     "TFG"     "ATP1B3"  "EFNA5"   "GABRP"   "CALML5"  "MMP7"    "SNORD68"
[11] "APP"     "SDF2L1"  "SAT1"   

message("Genes that also apeared in the original score:")

Genes that also apeared in the original score:

base::intersect(top_lum,original_lum_genes) 

character(0)

lum_enrich_res = genes_vec_enrichment(genes = top_lum,background = VariableFeatures(acc1_cancer_cells) %>% head(n_vargenes),homer = T,title = "lum top enrichment",custom_pathways = luminal_gs)

lum_enrich_res

calculate_score(dataset = acc1_cancer_cells,myo_genes = top_myo,lum_genes = top_lum,lum_threshold = 2,myo_threshold = 1)

correlation of lum score and myo score: -0.02
correlation of lum score and original lum score: 0.65
correlation of myo score and original myo score: 0.77

myo_intersected = intersect(top_myo,original_myo_genes) 
lum_intersected = intersect(top_lum,original_lum_genes) 
message("genes in myo score:")

genes in myo score:

myo_intersected

[1] "ACTA2" "DKK3"  "TAGLN"

message("genes in lum score:")

genes in lum score:

lum_intersected

[1] "LGALS3"

calculate_score(dataset = acc1_cancer_cells,myo_genes = myo_intersected,lum_genes = lum_intersected,lum_threshold = 2,myo_threshold = 1)

correlation of lum score and myo score: 0.03
correlation of lum score and original lum score: 0.69
correlation of myo score and original myo score: 0.81

enriched genes

rownames(lum_enrich_res) = lum_enrich_res$pathway_name
lum_enriched_genes = lum_enrich_res[3,"geneID"] %>% strsplit(split = "/") %>% .[[1]] %>% c(.,lum_protein_markers) #add original markers

rownames(myo_enrich_res) = myo_enrich_res$pathway_name
myo_enriched_genes = myo_enrich_res[3,"geneID"] %>% strsplit(split = "/") %>% .[[1]] %>% c(.,myo_protein_markers) #add original markers

message("genes in myo score:")

genes in myo score:

myo_enriched_genes

 [1] "FBLIM1"  "NEXN"    "NMNAT2"  "MYLK"    "CCDC50"  "IGFBP7"  "RAI14"   "ARAP3"   "SPARC"   "CALD1"   "LOXL2"  
[12] "COL5A1"  "ACTA2"   "DKK3"    "MSRB3"   "COL4A1"  "ACTN1"   "TPM1"    "TGFB1I1" "ADORA2B" "MXRA7"   "CNN1"   
[23] "TP63"    "ACTA2"  

message("genes in lum score:")

genes in lum score:

lum_enriched_genes

 [1] "PADI2"      "PATJ"       "PEX11B"     "APH1A"      "C1orf43"    "EFNA4"      "NECTIN4"    "SCYL3"      "ELF3"      
[10] "SOX13"      "IRF6"       "MED28"      "EPB41L4A"   "RGL2"       "C6orf132"   "TPD52L1"    "ICA1"       "MACC1"     
[19] "TRPS1"      "FAM83H"     "RASEF"      "ARRDC1"     "COMMD3"     "ANK3"       "GSTO2"      "PDCD4"      "EHF"       
[28] "ALDH3B2"    "SHANK2"     "SORL1"      "FKBP4"      "PTPN6"      "DUSP16"     "RERG"       "ADCY6"      "ERBB3"     
[37] "ERP29"      "SUSD6"      "RPS6KA5"    "SPINT1"     "FEM1B"      "TLE3"       "SCAMP2"     "CLN3"       "ADGRG1"    
[46] "ATP2C2"     "GGT6"       "MYO1D"      "ST6GALNAC2" "CYB5A"      "BLVRB"      "VRK3"       "SYCP2"      "TMPRSS2"   
[55] "LIMK2"      "KIT"       

calculate_score(dataset = acc1_cancer_cells,myo_genes = myo_enriched_genes,lum_genes = lum_enriched_genes,lum_threshold = 0,myo_threshold = -1)

correlation of lum score and myo score: 0.1
correlation of lum score and original lum score: 0.71
correlation of myo score and original myo score: 0.76

enriched genes and in original score

myo_enriched_genes = myo_enriched_genes[myo_enriched_genes %in% original_myo_genes]
lum_enriched_genes = lum_enriched_genes[lum_enriched_genes %in% original_lum_genes]

message("genes in myo score:")

genes in myo score:

myo_enriched_genes

[1] "MYLK"  "ACTA2" "DKK3"  "TP63"  "ACTA2"

message("genes in lum score:")

genes in lum score:

lum_enriched_genes

[1] "EHF" "KIT"

calculate_score(dataset = acc1_cancer_cells,myo_genes = myo_enriched_genes,lum_genes = lum_enriched_genes,lum_threshold = 2,myo_threshold = 2)

correlation of lum score and myo score: -0.07
correlation of lum score and original lum score: 0.62
correlation of myo score and original myo score: 0.77

HPV

Only HMSC cancer cells:

HPV33_P3 = fread("./Data/HPV33_P3.txt",col.names = c("plate","reads")) %>% as.data.frame()
HPV33_P3.df = HPV33_P3 %>% mutate(
  plate = gsub(x =HPV33_P3$plate, replacement = "",pattern = "_.*$") 
  %>% gsub(pattern = "-P",replacement = ".P") 
  %>% gsub(pattern = "-",replacement = "_",)
  )
HPV33_P3.df = HPV33_P3.df %>% dplyr::filter(HPV33_P3.df$plate %in% colnames(acc1_cancer_cells))
rownames(HPV33_P3.df)  <- HPV33_P3.df$plate
HPV33_P3.df$plate = NULL


HPV33_P2 = fread("./Data/HPV33_P2.txt",col.names = c("plate","reads")) %>% as.data.frame()
HPV33_P2.df = HPV33_P2 %>% mutate(
  plate = gsub(x =HPV33_P2$plate, replacement = "",pattern = "_.*$") 
  %>% gsub(pattern = "plate2-",replacement = "plate2_",)
  %>% gsub(pattern = "-",replacement = "\\.",)
  )
HPV33_P2.df = HPV33_P2.df %>% dplyr::filter(HPV33_P2.df$plate %in% colnames(acc1_cancer_cells))
rownames(HPV33_P2.df)  <- HPV33_P2.df$plate
HPV33_P2.df$plate = NULL

HPV33 = rbind(HPV33_P3.df,HPV33_P2.df)
acc1_cancer_cells = AddMetaData(object = acc1_cancer_cells,metadata = HPV33,col.name = "HPV33.reads")
FeaturePlot(acc1_cancer_cells,features = "HPV33.reads",max.cutoff = 10)

library(plotly)
data = FetchData(object = acc1_cancer_cells,vars = "HPV33.reads")
data = data$HPV33
data = data[data<100]
plot_ly(x = data, type = "histogram", nbinsx = 10)

data = FetchData(object = acc1_cancer_cells,vars = "HPV33.reads")
print(
  data %>% 
  ggplot(aes( x=HPV33.reads)) + 
  geom_histogram() + scale_x_continuous(breaks=c(1,3,10,100))
)

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

hpv33_positive = HPV33 %>% dplyr::mutate(hpv33_positive = case_when(reads >= 10 ~ "positive",
                                                                    reads < 10 ~ "negative")
)



hpv33_positive$reads = NULL
acc1_cancer_cells = AddMetaData(object = acc1_cancer_cells,metadata = hpv33_positive)

DimPlot(object = acc1_cancer_cells,group.by  = c("hpv33_positive"),pt.size = 2)

cNMF

library(reticulate)

acc1_cancer_cells = FindVariableFeatures(object = acc1_cancer_cells,nfeatures = nfeatures)

Calculating gene variances
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Calculating feature variances of standardized and clipped values
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

from cnmf import cNMF
import numpy as np
nfeatures_name = r.nfeatures_name
name = 'HMSC_cNMF_'+nfeatures_name+'Kvargenes'
outdir = './Data/cNMF'
K_range = np.arange(3,10)
cnmf_obj = cNMF(output_dir=outdir, name=name)
counts_fn='./Data/cNMF/hmsc_expressionData_'+nfeatures_name+'Kvargenes.txt'
tpm_fn = counts_fn ## This is a weird case where because this dataset is not 3' end umi sequencing, we opted to use the TPM matrix as the input matrix rather than the count matrix

cnmf_obj.prepare(counts_fn=counts_fn, components=K_range, seed=14,tpm_fn=tpm_fn)

cnmf_obj.factorize(worker_i=0, total_workers=1)

cnmf_obj.combine()
cnmf_obj.k_selection_plot()

Save object

# import pickle
# f = open('./Data/cNMF/HMSC_cNMF_'+nfeatures_name+'Kvargenes/cnmf_obj.pckl', 'wb')
# pickle.dump(cnmf_obj, f)
# f.close()

Load object

from cnmf import cNMF
import pickle
nfeatures = "2K"
f = open('./Data/cNMF/HMSC_cNMF_' + nfeatures+ 'vargenes/cnmf_obj.pckl', 'rb')
cnmf_obj = pickle.load(f)
f.close()

selected_k = 4
density_threshold = 0.1
cnmf_obj.consensus(k=selected_k, density_threshold=density_threshold,show_clustering=True)

/sci/labs/yotamd/lab_share/avishai.wizel/python_envs/miniconda/envs/cnmf_env_6/lib/python3.7/site-packages/scanpy/preprocessing/_simple.py:843: UserWarning: Received a view of an AnnData. Making a copy.
  view_to_actual(adata)

usage_norm, gep_scores, gep_tpm, topgenes = cnmf_obj.load_results(K=selected_k, density_threshold=density_threshold)

gep_scores = py$gep_scores
gep_tpm = py$gep_tpm
all_metagenes= py$usage_norm

Enrichment analysis by top 200 genes of each program

plt_list = list()
for (i in 1:ncol(gep_scores)) {
  top_genes = gep_scores  %>%  arrange(desc(gep_scores[i])) #sort by score a
  top = head(rownames(top_genes),200) #take top top_genes_num
  res = genes_vec_enrichment(genes = top,background = rownames(gep_scores),homer = T,title = 
                    i,silent = T,return_all = T)
   
  plt_list[[i]] = res$plt
}
gridExtra::grid.arrange(grobs = plt_list)

canonical_pathways = msigdbr(species = "Homo sapiens", category = "C2") %>% dplyr::filter(gs_subcat != "CGP") %>%  dplyr::distinct(gs_name, gene_symbol) 

plt_list = list()
for (i in 1:ncol(gep_scores)) {
  top_genes = gep_scores  %>%  arrange(desc(gep_scores[i])) #sort by score a
  top = head(rownames(top_genes),200) #take top top_genes_num
  res = genes_vec_enrichment(genes = top,background = rownames(gep_scores),homer = T,title = 
                    i,silent = T,return_all = T,custom_pathways = canonical_pathways)
   
  plt_list[[i]] = res$plt
}
gridExtra::grid.arrange(grobs = plt_list)

# lum genes in metagenes
message("lum genes in metagenes")

lum genes in metagenes

for (i in 1:ncol(gep_scores)) {
  top_genes = gep_scores  %>%  arrange(desc(gep_scores[i])) #sort by score a
  top = head(rownames(top_genes),200) #take top top_genes_num
  cat(paste0("metagene ",i,": "))
  print(original_lum_genes[original_lum_genes %in% top])

}

metagene 1: [1] "KRT7" "SLPI"
metagene 2: [1] "CLDN4"
metagene 3: character(0)
metagene 4: [1] "KRT7"   "LGALS3" "LCN2"   "SLPI"  

cat("\n")

# myo genes in metagenes
message("myo genes in metagenes")

myo genes in metagenes

for (i in 1:ncol(gep_scores)) {
  top_genes = gep_scores  %>%  arrange(desc(gep_scores[i])) #sort by score a
  top = head(rownames(top_genes),200) #take top top_genes_num
  cat(paste0("metagene ",i,": "))
  print(original_myo_genes[original_myo_genes %in% top])

}

metagene 1: [1] "KRT14" "ACTA2" "DKK3" 
metagene 2: character(0)
metagene 3: character(0)
metagene 4: character(0)

cat("\n")

notch_genes = c("JAG1","JAG2","NOTCH3","NOTCH2","NOTCH1","DLL1","MYB","HES4","HEY1","HEY2","NRARP")
# notch genes in metagenes
message("myo genes in metagenes")

myo genes in metagenes

for (i in 1:ncol(gep_scores)) {
  top_genes = gep_scores  %>%  arrange(desc(gep_scores[i])) #sort by score a
  top = head(rownames(top_genes),200) #take top top_genes_num
  cat(paste0("metagene ",i,": "))
  print(notch_genes[notch_genes %in% top])

}

metagene 1: character(0)
metagene 2: character(0)
metagene 3: character(0)
metagene 4: character(0)

cat("\n")

plt_list = list()
for (i in 1:ncol(gep_scores)) {
  top_genes = gep_scores  %>%  arrange(desc(gep_scores[i])) #sort by score a
  top = head(rownames(top_genes),10) #take top top_genes_num
  message(paste("program ",i,"top genes:"))
  print(top)

}

program  1 top genes:

 [1] "IGKV5-2" "NDUFB7"  "LGALS1"  "UBL5"    "S100A6"  "S100A11" "CUTA"    "IFITM3"  "JTB"     "IL17RB" 

program  2 top genes:

 [1] "EGR1"    "C6orf62" "JUNB"    "DNAJA1"  "ATF3"    "IER2"    "ERRFI1"  "KLF6"    "CDKN1A"  "MTHFD2" 

program  3 top genes:

 [1] "RP1-128M12.3"  "RP11-374F3.3"  "RP11-403A21.3" "RP11-454L9.2"  "AC097500.1"    "RP11-463H24.4" "RP11-515O17.2"
 [8] "RP11-536L3.4"  "PALD1"         "AL049758.2"   

program  4 top genes:

 [1] "LTF"           "PIGR"          "FMO2"          "HSD17B2"       "MLPH"          "PRR15L"        "RF00019.219"  
 [8] "RP11-817J15.2" "CD14"          "CLU"          

meta3 = FetchData(object = acc1_cancer_cells,vars = c("metagene.3"))
ggplot(meta3, aes(x=metagene.3)) + 
  geom_density()

meta3[,1] = as.numeric(meta3[,1])
sum(meta3[,1]>0,na.rm = T )

[1] 82

sum(meta3[,1]==0,na.rm = T )

[1] 53

assignment UMAP

larger_by = 2
message(paste("larger_by = ", larger_by))

larger_by =  2

acc1_cancer_cells = program_assignment(dataset = acc1_cancer_cells,larger_by = larger_by,program_names = colnames(all_metagenes))
selected_k =4
colors =  rainbow(selected_k)
colors = c(colors,"grey")
DimPlot(acc1_cancer_cells,group.by = "program.assignment",pt.size = 2,cols =colors)

show cell cycle program:

hallmark_name = "GO_MITOTIC_CELL_CYCLE"
genesets  =GSEABase::getGmt("./Data/h.all.v7.0.symbols.pluscc.gmt")
var_features=acc1_cancer_cells@assays$RNA@var.features
geneIds= genesets[[hallmark_name]]@geneIds
score <- apply(acc1_cancer_cells@assays$RNA@data[intersect(geneIds,var_features),],2,mean)
acc1_cancer_cells=AddMetaData(acc1_cancer_cells,score,hallmark_name)

FeaturePlot(object = acc1_cancer_cells,features = hallmark_name)

cc_vs_program2 = FetchData(object = acc1_cancer_cells,vars = c("metagene.2",hallmark_name))
cor(cc_vs_program2[1],cc_vs_program2[2])

           GO_MITOTIC_CELL_CYCLE
metagene.2             0.6127281

Comparisions

cnv_vs_hpv = FetchData(object = acc1_cancer_cells,vars = c("cnv.cluster","hpv33_positive"))
test <- fisher.test(table(cnv_vs_hpv))
ggbarstats(
  cnv_vs_hpv, cnv.cluster, hpv33_positive,
  results.subtitle = FALSE,
  subtitle = paste0(
    "Fisher's exact test", ", p-value = ",
    ifelse(test$p.value < 0.001, "< 0.001", round(test$p.value, 3))
  )
)

cnv_vs_cnmf = FetchData(object = acc1_cancer_cells,vars = c("program.assignment","cnv.cluster"))
cnv_vs_cnmf = cnv_vs_cnmf %>% dplyr::filter(program.assignment == "1" |program.assignment == "2" )
test <- fisher.test(table(cnv_vs_cnmf))
ggbarstats(
  cnv_vs_cnmf, program.assignment, cnv.cluster,
  results.subtitle = FALSE,
  subtitle = paste0(
    "Fisher's exact test", ", p-value = ",
    ifelse(test$p.value < 0.001, "< 0.001", round(test$p.value, 3))
  )
)

cnmf_vs_hpv = FetchData(object = acc1_cancer_cells,vars = c("program.assignment","hpv33_positive"))
cnmf_vs_hpv = cnmf_vs_hpv %>% dplyr::filter(program.assignment == "1" |program.assignment == "2" )
test <- fisher.test(table(cnmf_vs_hpv))
ggbarstats(
  cnmf_vs_hpv, program.assignment, hpv33_positive,
  results.subtitle = FALSE,
  subtitle = paste0(
    "Fisher's exact test", ", p-value = ",
    ifelse(test$p.value < 0.001, "< 0.001", round(test$p.value, 3))
  )
)

myb_vs_cnv = FetchData(object = acc1_cancer_cells,vars = c("cnv.cluster","MYB"))
myb_vs_cnv $cnv.cluster = as.character(myb_vs_cnv $cnv.cluster )

ggboxplot(myb_vs_cnv, x = "cnv.cluster", y = "MYB",
          palette = "jco",
          add = "jitter")+ stat_compare_means(method = "wilcox.test",comparisons = list(c("1","2")))

myb_vs_hpv = FetchData(object = acc1_cancer_cells,vars = c("hpv33_positive","MYB"))
myb_vs_hpv $hpv33_positive = as.character(myb_vs_hpv $hpv33_positive )

ggboxplot(myb_vs_hpv, x = "hpv33_positive", y = "MYB",
          palette = "jco",
          add = "jitter")+ stat_compare_means(method = "wilcox.test",comparisons = list(c("positive","negative")))+ stat_summary(fun.data = function(x) data.frame(y=15, label = paste("Mean=",round(mean(x),digits = 2))), geom="text") +ylab("log2(MYB)")

hpvReads_vs_myb = FetchData(object = acc1_cancer_cells,vars = c("HPV33.reads","MYB"))
corr = cor(hpvReads_vs_myb$HPV33.reads,hpvReads_vs_myb$MYB)
print("correlation of MYB abd hpv33_reads:" %>% paste(corr %>% round(digits = 2)))

acc1_cancer_cells_data = acc1_cancer_cells@assays[["RNA"]]@data %>% as.data.frame()
acc1_cancer_cells_data = cor(acc1_cancer_cells_data)
annotation = FetchData(object = acc1_cancer_cells,vars = c("orig.ident"))

colors <- c(seq(-1,1,by=0.01))
  my_palette <- c("blue",colorRampPalette(colors = c("blue", "white", "red"))
                                                   (n = length(colors)-3), "red")
  
pht1 = pheatmap(acc1_cor,annotation_col  = annotation,fontsize = 6,breaks = colors, color = my_palette,show_colnames = F,show_rownames = F)

acc1_cancer_cells_data = acc1_cancer_cells@assays[["RNA"]]@data %>% as.data.frame()
acc1_cancer_cells_data = cor(acc1_cancer_cells_data)
annotation = FetchData(object = acc1_cancer_cells,vars = c("cnv.cluster"))

colors <- c(seq(-1,1,by=0.01))
  my_palette <- c("blue",colorRampPalette(colors = c("blue", "white", "red"))
                                                   (n = length(colors)-3), "red")
  
pht1 = pheatmap(acc1_cor,annotation_col  = annotation,fontsize = 6,breaks = colors, color = my_palette,show_colnames = F,show_rownames = F)

cnv_vs_plate = FetchData(object = acc1_cancer_cells,vars = c("cnv.cluster","orig.ident"))
test <- fisher.test(table(cnv_vs_plate))
ggbarstats(
  cnv_vs_plate, cnv.cluster, orig.ident,
  results.subtitle = FALSE,
  subtitle = paste0(
    "Fisher's exact test", ", p-value = ",
    ifelse(test$p.value < 0.001, "< 0.001", round(test$p.value, 3))
  )
)

# creat colours for each group
cnv_vs_plate = FetchData(object = acc1_cancer_cells,vars = c("cnv.cluster","orig.ident"))

rownames(cnv_vs_plate) = rownames(cnv_vs_plate) %>% gsub(pattern = "_",replacement = "\\.")
cnv_vs_plate$cnv.cluster = as.character(cnv_vs_plate$cnv.cluster)
cnv_vs_plate = cnv_vs_plate %>% dplyr::rename(plate = "orig.ident")
annoCol <- list(plate = c(ACC1.P3 = "red",ACC.plate2 = "green"),cnv.cluster =c("1"="green","2" = "red"))

pheatmap(infercnv.observations2,cluster_cols = F,cluster_rows = F, show_rownames = F,show_colnames = F, breaks = breaks,color = colorRampPalette(rev(c("darkred", "white", "darkblue")))(15),annotation_row = cnv_vs_plate,annotation_colors = annoCol)

cnmf_vs_plate = FetchData(object = acc1_cancer_cells,vars = c("program.assignment","orig.ident"))
cnmf_vs_plate= cnmf_vs_plate %>% dplyr::filter(program.assignment == "1" | program.assignment == "2")
test <- fisher.test(table(cnmf_vs_plate))
ggbarstats(
  cnmf_vs_plate, program.assignment, orig.ident,
  results.subtitle = FALSE,
  subtitle = paste0(
    "Fisher's exact test", ", p-value = ",
    ifelse(test$p.value < 0.001, "< 0.001", round(test$p.value, 3))
  )
)

plate_3 = cnmf_vs_plate %>% dplyr::filter((program.assignment == 1 & orig.ident == "ACC1.P3") ) %>% rownames() 
plate_2 = cnmf_vs_plate %>% dplyr::filter((program.assignment == 2 & orig.ident == "ACC.plate2")) %>% rownames()
cells = list(ACC1.P3 = plate_3,ACC.plate2  = plate_2)

Results

exceptions

exceptions_plt = DimPlot(object = acc1_cancer_cells, cells.highlight = cells, cols.highlight = c("green","red"), cols = "gray", order = TRUE,pt.size = 2,sizes.highlight = 2,combine = F) 

Warning messages: 1: No Python documentation found for ‘cnmf_obj.prepare’. 2: No Python documentation found for ‘cnmf_obj.prepare’.

colors = c("green","red","cyan","purple","grey")
program.assignment_plt = DimPlot(acc1_cancer_cells,group.by = "program.assignment",pt.size = 2,cols = colors,combine = F)
metagene.1_plt <- FeaturePlot(object = acc1_cancer_cells,features = c("metagene.1"),combine = F)
metagene.2_plt = FeaturePlot(object = acc1_cancer_cells,features = c("metagene.2"),combine = F)

lst = list(exceptions = exceptions_plt, program.assignment = program.assignment_plt,metagene.1 = metagene.1_plt,metagene.2 = metagene.2_plt)

for (i in 1:length(lst)) {
  cat("### ",names(lst)[i]," \n")
  print(
    lst[[i]]
    )
  plot.new()
  dev.off()
  cat(' \n\n')
}

exceptions

[[1]] Warning in grSoftVersion() : unable to load shared object ‘/usr/local/lib/R/modules//R_X11.so’: libXt.so.6: cannot open shared object file: No such file or directory

program.assignment

[[1]]

metagene.1

[[1]]

metagene.2

[[1]]

NA

no_neg <- function(x) {
  x = x + abs(min(x))
  x
}
# gep_scores_norm = apply(gep_scores, MARGIN = 2, FUN = min_max_normalize)%>% as.data.frame()
# gep_scores_norm = gep_scores
gep_scores_norm = apply(gep_scores, MARGIN = 2, FUN = no_neg)%>% as.data.frame()
gep_scores_norm = sum2one(gep_scores_norm)
all_metagenes = expression_mult(gep_scores = gep_scores_norm,dataset = all_acc_cancer_cells,top_genes = T,z_score = F)

# Make metagene names
for (i in 1:ncol(all_metagenes)) {
  colnames(all_metagenes)[i] = "metagene." %>% paste0(i)
}


#add each metagene to metadata
for (i in 1:ncol(all_metagenes)) {
  metage_metadata = all_metagenes %>% select(i)
  all_acc_cancer_cells = AddMetaData(object = all_acc_cancer_cells,metadata = metage_metadata)
}

FeaturePlot(object = all_acc_cancer_cells,features = colnames(all_metagenes),max.cutoff = 100)