patients_cnmf_analysis

1 Data

lung = readRDS("./Data/lung_cancercells_withTP_onlyPatients.rds")
lung_patients = lung$patient.ident %>% unique() %>% as.character()
lung_patients_filtered = lung_patients[!(lung_patients %in% c("X1055new","X1099"))] # remove patients with less than 100 malignant cells
lung = subset(x = lung,subset = patient.ident %in% lung_patients_filtered)

suffix ="xeno_genes_normalized_0-5sigma_2-7theta"

from cnmf import cNMF
suffix = r.suffix
import pickle
f = open('./Data/cnmf/cnmf_objects/patients_' + suffix + '_cnmf_obj.pckl', 'rb')
cnmf_obj = pickle.load(f)
f.close()

2 Functions

library(stringi)
library(reticulate)
source_from_github(repositoy = "DEG_functions",version = "0.2.24")
source_from_github(repositoy = "cNMF_functions",version = "0.3.9",script_name = "cnmf_function_Harmony.R") 

3 K selection plot

plot_path = paste0("/sci/labs/yotamd/lab_share/avishai.wizel/R_projects/EGFR/Data/cnmf/cNMF_patients_Varnorm_Harmony_",suffix,"/cNMF_patients_Varnorm_Harmony_",suffix,".k_selection.png")
knitr::include_graphics(plot_path)

4 gep scores for all NMF k’s

density_threshold = 0.1
usage_norm, gep_scores3, gep_tpm, topgenes = cnmf_obj.load_results(K=3, density_threshold=density_threshold)
usage_norm, gep_scores4, gep_tpm, topgenes = cnmf_obj.load_results(K=4, density_threshold=density_threshold)
usage_norm, gep_scores5, gep_tpm, topgenes = cnmf_obj.load_results(K=5, density_threshold=density_threshold)
usage_norm, gep_scores6, gep_tpm, topgenes = cnmf_obj.load_results(K=6, density_threshold=density_threshold)
usage_norm, gep_scores7, gep_tpm, topgenes = cnmf_obj.load_results(K=7, density_threshold=density_threshold)
usage_norm, gep_scores8, gep_tpm, topgenes = cnmf_obj.load_results(K=8, density_threshold=density_threshold)
usage_norm, gep_scores9, gep_tpm, topgenes = cnmf_obj.load_results(K=9, density_threshold=density_threshold)

5 Enrichment analysis by top 200 genes of each program

gep_scores3 = py$gep_scores3
gep_scores4 = py$gep_scores4
gep_scores5 = py$gep_scores5
gep_scores6 = py$gep_scores6
gep_scores7 = py$gep_scores7
gep_scores8 = py$gep_scores8
gep_scores9 = py$gep_scores9

all_gep_scores =  list(gep_scores3 = gep_scores3, gep_scores4 = gep_scores4, gep_scores5 = gep_scores5, gep_scores6 = gep_scores6, gep_scores7= gep_scores7, gep_scores8 = gep_scores8, gep_scores9 = gep_scores9)
# canonical_pathways = msigdbr(species = "Homo sapiens", category = "C2") %>% dplyr::filter(gs_subcat != "CGP") %>%  dplyr::distinct(gs_name, gene_symbol) 
for (gep_name in names(all_gep_scores)) {
  gep_scores = all_gep_scores[[gep_name]]
  top_genes_num = 200
  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),top_genes_num) #take top top_genes_num
    res = genes_vec_enrichment(genes = top,background = rownames(gep_scores),homer = T,title = 
                      names(gep_scores)[i],silent = T,return_all = T,custom_pathways = NULL )
     
    plt_list[[i]] = res$plt
  }
  print_tab(plt =   ggarrange(plotlist = plt_list),
            title = gep_name)
}

gep_scores3

gep_scores4

gep_scores5

gep_scores6

gep_scores7

gep_scores8

gep_scores9

NA

6 Chosen K

selected_k = 8
print("selected k = ",selected_k)

selected k =  8

density_threshold = 0.1
cnmf_obj.consensus(k=selected_k, density_threshold=density_threshold,show_clustering=True)
usage_norm, gep_scores, gep_tpm, topgenes = cnmf_obj.load_results(K=selected_k, density_threshold=density_threshold)

patients_geps = py$gep_scores

6.1 programs enrichment with canonical

gep_scores = py$gep_scores

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

top_genes_num = 200
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),top_genes_num) #take top top_genes_num
  res = genes_vec_enrichment(genes = top,background = rownames(gep_scores),homer = T,title = 
                    names(gep_scores)[i],silent = T,return_all = T,custom_pathways = canonical_pathways )
   
  plt_list[[i]] = res$plt
}
gridExtra::grid.arrange(grobs = plt_list)

7 Correlation of programs

cor_res = cor(patients_geps)
breaks <- c(seq(-1,1,by=0.01))
colors_for_plot <- colorRampPalette(colors = c("blue", "white", "red"))(n = length(breaks))

pht = pheatmap(cor_res,color = colors_for_plot,breaks = breaks)

print_tab(pht,title = "correlation")

##   correlation {.unnumbered }  
 

#correlation based on top 200
all_top= c()
for (i in 1:ncol(patients_geps)) {
  top_genes = patients_geps  %>%  arrange(desc(patients_geps[i])) #sort by score a
  top = head(rownames(top_genes),200) #take top top_genes_num
all_top = c(all_top,top)
}
gep_scores_top = patients_geps[all_top,]
cor_res = cor(gep_scores_top)

breaks <- c(seq(-1,1,by=0.01))
colors_for_plot <- colorRampPalette(colors = c("blue", "white", "red"))(n = length(breaks))

# correlation by top 150 combined
pht = pheatmap(cor_res,color = colors_for_plot,breaks = breaks)

print_tab(pht,title = "top 200 genes correlation")

##   top 200 genes correlation {.unnumbered }  
 

8 Combine similar programs

groups_list = c(4,3,6)
patients_geps = union_programs(groups_list = groups_list,all_metagenes = patients_geps)

plt_list = list()

for (program  in names (patients_geps)) {
 p = ggplot(patients_geps, aes(x=!!ensym(program))) +
  geom_density()+xlab(program)+
   geom_vline(
    aes(xintercept=sort(patients_geps[,program],TRUE)[200]  ,color="top200"),
          linetype="dashed", size=1)+
   geom_vline(
    aes(xintercept=sort(patients_geps[,program],TRUE)[100]  ,color="top100"),
          linetype="dashed", size=1)+
      geom_vline(
    aes(xintercept=sort(patients_geps[,program],TRUE)[50]  ,color="top50"),
          linetype="dashed", size=1)+
         geom_vline(
    aes(xintercept=sort(patients_geps[,program],TRUE)[150]  ,color="top150"),
          linetype="dashed", size=1)+
   scale_color_manual(name = "top n genes", values = c(top200 = "blue",top100 = "red",top150 = "yellow",top50 = "green"))
   plt_list[[program]] <- p

}
 
ggarrange(plotlist = plt_list)

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

9 Calculate usage

# get expression with genes in cnmf input
lung = FindVariableFeatures(object = lung,nfeatures = 2000)
genes = rownames(lung)[rownames(lung) %in% VariableFeatures(object = xeno)[1:2000]]

lung_expression = t(as.matrix(GetAssayData(lung,slot='data'))) 
lung_expression = 2**lung_expression #convert from log2(tpm+1) to tpm
lung_expression = lung_expression-1
lung_expression = lung_expression[,genes] %>% as.data.frame()

all_0_genes = colnames(lung_expression)[colSums(lung_expression==0, na.rm=TRUE)==nrow(lung_expression)] #delete rows that have all 0
genes = genes[!genes %in% all_0_genes]
lung_expression = lung_expression[,!colnames(lung_expression) %in% all_0_genes]
gc(verbose = F)

lung_expression = r.lung_expression
genes = r.genes
gep_scores = r.patients_geps
usage_by_calc = get_usage_from_score(counts=lung_expression,tpm=lung_expression,genes=genes,cnmf_obj=cnmf_obj,k=selected_k,sumTo1=False)

usage_by_calc = py$usage_by_calc
groups_list = c(4,3,6)
usage_by_calc = union_programs(groups_list = groups_list,all_metagenes = usage_by_calc)
usage_by_calc = apply(usage_by_calc, MARGIN = 1, sum_2_one) %>% t() %>% as.data.frame()
usage_by_calc =usage_by_calc %>% rename(cell_cycle = gep4.3.6, hypoxia_like = gep2, interferon_like = gep1, TNFa =  gep5, INF2 = gep7)

10 Usgage UMAP

all_metagenes= usage_by_calc

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

Note: Using an external vector in selections is ambiguous.
ℹ Use `all_of(i)` instead of `i` to silence this message.
ℹ See <https://tidyselect.r-lib.org/reference/faq-external-vector.html>.
This message is displayed once per session.

FeaturePlot(object = lung,features = colnames(all_metagenes),ncol = 3)

Warning: The following variables were found in both object metadata and the default assay: INF2
Returning metadata; if you want the feature, please use the assay's key (eg. rna_INF2)

DimPlot(object = lung,group.by = "time.point",pt.size = 0.5)
pre_treatment_cells = FetchData(object = lung,vars = "time.point") %>% filter(time.point == "pre-treatment") %>% rownames()
on_treatment_cells = FetchData(object = lung,vars = "time.point") %>% filter(time.point == "on-treatment") %>% rownames()

cells_to_highlight =  list(pre_treatment_cells = pre_treatment_cells)
DimPlot(object = lung, cells.highlight = cells_to_highlight, cols.highlight = c("red"), cols = "gray", order = TRUE, label = T, repel = T)

cells_to_highlight =  list(on_treatment_cells = on_treatment_cells)
DimPlot(object = lung, cells.highlight = cells_to_highlight, cols.highlight = c("cyan3"), cols = "gray", order = TRUE,  label = T, repel = T)

DotPlot(object = lung, features = c("hypoxia_like","interferon_like","cell_cycle","TNFa", "INF2"),scale = T,group.by  = 'time.point')

Warning: The following variables were found in both object metadata and the default assay: INF2
Returning metadata; if you want the feature, please use the assay's key (eg. rna_INF2)

11 Assignment

larger_by = 1.25
lung = program_assignment(dataset = lung,larger_by = larger_by,program_names = colnames(all_metagenes))

p = cell_percentage(dataset = lung,time.point_var = "time.point",by_program = T)
print_tab(plt = p,title = "by program")
p = cell_percentage(dataset = lung,time.point_var = "time.point",by_tp = T,x_order = NULL)
print_tab(plt = p,title = "by timepoint")

# colors =  rainbow(ncol(all_metagenes))
# fc <- colorRampPalette(c("lightgreen", "darkgreen"))
# greens = fc(4)
# colors[1] = "blue"
# colors[2:3] = greens[1:2]
# colors[4] = "red"
# colors[5:6] = greens[3:4]
# colors = c(colors,"grey")
# DimPlot(lung,group.by = "program.assignment",pt.size = 0.5,cols =colors)


# colors =  rainbow(ncol(all_metagenes))
# colors = c(colors,"grey")
# DimPlot(lung,group.by = "program.assignment",pt.size = 0.5,cols =colors)

DimPlot(lung,group.by = "program.assignment",pt.size = 0.5)
DimPlot(lung,group.by = "patient.ident",pt.size = 0.5)
DimPlot(lung,group.by = "time.point",pt.size = 0.5)

12 Score regulation

metagenes_mean_compare(dataset = lung, time.point_var = "time.point",prefix = "patient",patient.ident_var = "patient.ident",pre_on = c("pre-treatment","on-treatment"),axis.text.x = 8,programs = c("hypoxia_like","interferon_like","cell_cycle"))

hypoxia_like per patient

hypoxia_like

interferon_like per patient

interferon_like

cell_cycle per patient

cell_cycle

NA

13 CC signature

# get genes and plot cor heatmap
hallmark_name = "HALLMARK_G2M_CHECKPOINT"
genesets  =getGmt("./Data/h.all.v7.0.symbols.pluscc.gmt")
geneIds= genesets[[hallmark_name]]@geneIds
hallmars_exp = FetchData(object = lung,vars = c(geneIds))

Warning in FetchData.Seurat(object = lung, vars = c(geneIds)) : The following requested variables were not found: AC027237.1, JPT1, KNL1, TENT4A, AC091021.1, PTTG3P

hallmars_exp = hallmars_exp[,colSums(hallmars_exp[])>0] #remove no expression genes
hallmark_cor = cor(hallmars_exp)
pht1 = pheatmap(mat = hallmark_cor,silent = T)

# make annotations
num_of_clusters = 4
clustering_distance = "euclidean"
myannotation = as.data.frame(cutree(pht1[["tree_row"]], k = num_of_clusters)) #split into k clusters
names(myannotation)[1] = "cluster"
myannotation$cluster = as.factor(myannotation$cluster)
palette1 <-brewer.pal(num_of_clusters, "Paired")
names(palette1) = unique(myannotation$cluster)
ann_colors = list (cluster = palette1)
annotation = list(ann_colors = ann_colors, myannotation = myannotation)

#choose colors
colors <- c(seq(-1,1,by=0.01))
my_palette <- c("blue",colorRampPalette(colors = c("blue", "white", "red"))
                                                 (n = length(colors)-3), "red")


print_tab(plt = 
              pheatmap(mat = hallmark_cor,annotation_col =  annotation[["myannotation"]], annotation_colors = annotation[["ann_colors"]], clustering_distance_rows = clustering_distance,clustering_distance_cols = clustering_distance,color = my_palette,breaks = colors,show_rownames = F,show_colnames = F)
            ,title = "genes expression heatmap")

genes expression heatmap

NA

#choose clusters
chosen_clusters = c(1,2,3,4)
print ("chosen_clusters= ", chosen_clusters)

[1] “chosen_clusters=”

#UMAP expression of signature
chosen_genes = annotation[["myannotation"]] %>% dplyr::filter(cluster %in% chosen_clusters) %>% rownames() #take relevant genes
score <- apply(lung@assays$RNA@data[chosen_genes,],2,mean)
lung=AddMetaData(lung,score,hallmark_name)

print_tab(FeaturePlot(object = lung, features = hallmark_name),title = "Expression")

Expression

#plot signature distribution
cc_scores = FetchData(object = lung,vars = "HALLMARK_G2M_CHECKPOINT")
plt  =  ggplot(cc_scores, aes(x=HALLMARK_G2M_CHECKPOINT)) +
  geom_density()+
   geom_vline(
    aes(xintercept=mean(cc_scores$HALLMARK_G2M_CHECKPOINT) + sd(cc_scores$HALLMARK_G2M_CHECKPOINT) ,color="1 SD"),
          linetype="dashed", size=1)+
    geom_vline(
    aes(xintercept=mean(cc_scores$HALLMARK_G2M_CHECKPOINT) + 2*sd(cc_scores$HALLMARK_G2M_CHECKPOINT) ,color="2 SD"),
          linetype="dashed", size=1)

print_tab(plt = plt,title = "dist")

dist

Warning: Use of cc_scores$HALLMARK_G2M_CHECKPOINT is discouraged. Use HALLMARK_G2M_CHECKPOINT instead. Warning: Use of cc_scores$HALLMARK_G2M_CHECKPOINT is discouraged. Use HALLMARK_G2M_CHECKPOINT instead. Warning: Use of cc_scores$HALLMARK_G2M_CHECKPOINT is discouraged. Use HALLMARK_G2M_CHECKPOINT instead. Warning: Use of cc_scores$HALLMARK_G2M_CHECKPOINT is discouraged. Use HALLMARK_G2M_CHECKPOINT instead.

# Plot assignment 

cc_scores = cc_scores %>% mutate(is_cycling = if_else(condition = 
                                          HALLMARK_G2M_CHECKPOINT > mean(HALLMARK_G2M_CHECKPOINT) + sd(HALLMARK_G2M_CHECKPOINT),
                                        true = "cycling",
                                        false = "non_cycling"))
lung = AddMetaData(object = lung,metadata = cc_scores$is_cycling,col.name = "is_cycling")
print_tab(plt = DimPlot(object = lung,group.by = "is_cycling") , title = "assignment 1 sd")

assignment 1 sd

# 
# cc_scores = cc_scores %>% mutate(is_cycling = if_else(condition = 
#                                           HALLMARK_G2M_CHECKPOINT > mean(HALLMARK_G2M_CHECKPOINT) + 2*sd(cc_scores$HALLMARK_G2M_CHECKPOINT),
#                                         true = "cycling",
#                                         false = "non_cycling"))
# lung = AddMetaData(object = lung,metadata = cc_scores$is_cycling,col.name = "is_cycling")
# print_tab(plt = DimPlot(object = lung,group.by = "is_cycling") , title = "assignment 2 sd")
# 

 df  = FetchData(object = lung,vars = c("is_cycling","time.point")) %>% 
    filter (time.point %in% c("pre-treatment","on-treatment")) %>% 
    droplevels() 
  test = fisher.test(table(df))
    
  library(ggstatsplot)

    plt = ggbarstats(
    df, is_cycling, time.point,
    results.subtitle = FALSE,
    subtitle = paste0(
      "Fisher's exact test", ", p-value = ",
       round(test$p.value,13))
    )
  
print_tab(plt = plt,title = "fisher")

fisher

NA

#correlation to NMF

cor_res = cor(lung$cell_cycle,lung$HALLMARK_G2M_CHECKPOINT)
print(paste("correlation of cc prgoram to HALLMARK_G2M_CHECKPOINT:", cor_res))

[1] "correlation of cc prgoram to HALLMARK_G2M_CHECKPOINT: 0.77656941037506"

14 Hypoxia signature

# get genes and create heatmap
hallmark_name = "HALLMARK_HYPOXIA"
genesets  =getGmt("./Data/h.all.v7.0.symbols.pluscc.gmt")
geneIds= genesets[[hallmark_name]]@geneIds
hallmars_exp = FetchData(object = lung,vars = c(geneIds))

Warning in FetchData.Seurat(object = lung, vars = c(geneIds)) : The following requested variables were not found: CAVIN3, CCN5, LARGE1, CAVIN1, CCN1, CCN2

hallmars_exp = hallmars_exp[,colSums(hallmars_exp[])>0] #remove no expression genes
hallmark_cor = cor(hallmars_exp)
pht1 = pheatmap(mat = hallmark_cor,silent = T)

num_of_clusters = 4
clustering_distance = "euclidean"
myannotation = as.data.frame(cutree(pht1[["tree_row"]], k = num_of_clusters)) #split into k clusters
 
names(myannotation)[1] = "cluster"
  myannotation$cluster = as.factor(myannotation$cluster)
  
  palette1 <-brewer.pal(num_of_clusters, "Paired")

  names(palette1) = unique(myannotation$cluster)
  ann_colors = list (cluster = palette1)
  annotation = list(ann_colors = ann_colors, myannotation = myannotation)
  
  colors <- c(seq(-1,1,by=0.01))
  my_palette <- c("blue",colorRampPalette(colors = c("blue", "white", "red"))
                                                   (n = length(colors)-3), "red")


  print_tab(plt = 
                pheatmap(mat = hallmark_cor,annotation_col =  annotation[["myannotation"]], annotation_colors = annotation[["ann_colors"]], clustering_distance_rows = clustering_distance,clustering_distance_cols = clustering_distance,color = my_palette,breaks = colors,show_rownames = F,show_colnames = F)
            ,title = "genes expression heatmap")

genes expression heatmap

NA

chosen_genes = annotation[["myannotation"]] %>% dplyr::filter(cluster == 1 | cluster == 2) %>% rownames() #take relevant genes
score <- apply(lung@assays$RNA@data[chosen_genes,],2,mean)
lung=AddMetaData(lung,score,hallmark_name)

print_tab(FeaturePlot(object = lung, features = hallmark_name),title = "Expression")

Expression

NA

cc_scores = FetchData(object = lung,vars = hallmark_name)

plt  =  ggplot(cc_scores, aes(x=HALLMARK_HYPOXIA)) +
  geom_density()+
   geom_vline(
    aes(xintercept=mean(cc_scores$HALLMARK_HYPOXIA) + sd(cc_scores$HALLMARK_HYPOXIA) ,color="1 SD"),
          linetype="dashed", size=1)+
    geom_vline(
    aes(xintercept=mean(cc_scores$HALLMARK_HYPOXIA) + 2*sd(cc_scores$HALLMARK_HYPOXIA) ,color="2 SD"),
          linetype="dashed", size=1)

print_tab(plt = plt,title = "dist")

dist

Warning: Use of cc_scores$HALLMARK_HYPOXIA is discouraged. Use HALLMARK_HYPOXIA instead. Warning: Use of cc_scores$HALLMARK_HYPOXIA is discouraged. Use HALLMARK_HYPOXIA instead. Warning: Use of cc_scores$HALLMARK_HYPOXIA is discouraged. Use HALLMARK_HYPOXIA instead. Warning: Use of cc_scores$HALLMARK_HYPOXIA is discouraged. Use HALLMARK_HYPOXIA instead.

NA

14.1 correlation to nmf

cor_res = cor(lung$hypoxia_like,lung$HALLMARK_HYPOXIA)
print(paste("correlation of hypoxia program to HALLMARK_HYPOXIA:", cor_res))

[1] "correlation of hypoxia program to HALLMARK_HYPOXIA: 0.487874254360783"

15 glycolysis signature

hallmark_name = "HALLMARK_GLYCOLYSIS"
genesets  =getGmt("./Data/h.all.v7.0.symbols.pluscc.gmt")
geneIds= genesets[[hallmark_name]]@geneIds
hallmars_exp = FetchData(object = lung,vars = c(geneIds))

Warning in FetchData.Seurat(object = lung, vars = c(geneIds)) : The following requested variables were not found: AC010618.1, AC074143.1

hallmars_exp = hallmars_exp[,colSums(hallmars_exp[])>0] #remove no expression genes
hallmark_cor = cor(hallmars_exp)
pht1 = pheatmap(mat = hallmark_cor,silent = T)

num_of_clusters = 5
clustering_distance = "euclidean"
myannotation = as.data.frame(cutree(pht1[["tree_row"]], k = num_of_clusters)) #split into k clusters
 
names(myannotation)[1] = "cluster"
  myannotation$cluster = as.factor(myannotation$cluster)
  
  palette1 <-brewer.pal(num_of_clusters, "Paired")

  names(palette1) = unique(myannotation$cluster)
  ann_colors = list (cluster = palette1)
  annotation = list(ann_colors = ann_colors, myannotation = myannotation)
  
  colors <- c(seq(-1,1,by=0.01))
  my_palette <- c("blue",colorRampPalette(colors = c("blue", "white", "red"))
                                                   (n = length(colors)-3), "red")


  print_tab(plt = 
                pheatmap(mat = hallmark_cor,annotation_col =  annotation[["myannotation"]], annotation_colors = annotation[["ann_colors"]], clustering_distance_rows = clustering_distance,clustering_distance_cols = clustering_distance,color = my_palette,breaks = colors,show_rownames = F,show_colnames = F)
            ,title = "genes expression heatmap")

genes expression heatmap

NA

chosen_genes = annotation[["myannotation"]] %>% dplyr::filter(cluster == 1 | cluster == 2 | cluster == 5) %>% rownames() #take relevant genes
score <- apply(lung@assays$RNA@data[chosen_genes,],2,mean)
lung=AddMetaData(lung,score,hallmark_name)

print_tab(FeaturePlot(object = lung, features = hallmark_name),title = "Expression")

Expression

NA

cc_scores = FetchData(object = lung,vars = hallmark_name)

plt  =  ggplot(cc_scores, aes(x=HALLMARK_GLYCOLYSIS)) +
  geom_density()+
   geom_vline(
    aes(xintercept=mean(HALLMARK_GLYCOLYSIS) + sd(HALLMARK_GLYCOLYSIS) ,color="1 SD"),
          linetype="dashed", size=1)+
    geom_vline(
    aes(xintercept=mean(HALLMARK_GLYCOLYSIS) + 2*sd(HALLMARK_GLYCOLYSIS) ,color="2 SD"),
          linetype="dashed", size=1)+
    geom_vline(
    aes(xintercept=mean(HALLMARK_GLYCOLYSIS) ,color="mean"),
          linetype="dashed", size=1)

print_tab(plt = plt,title = "dist")

dist

NA

16 INF signature

hallmark_name = "HALLMARK_INTERFERON_GAMMA_RESPONSE"
genesets  =getGmt("./Data/h.all.v7.0.symbols.pluscc.gmt")
geneIds= genesets[[hallmark_name]]@geneIds
hallmars_exp = FetchData(object = lung,vars = c(geneIds))

Warning in FetchData.Seurat(object = lung, vars = c(geneIds)) : The following requested variables were not found: AC124319.1

hallmars_exp = hallmars_exp[,colSums(hallmars_exp[])>0] #remove no expression genes
hallmark_cor = cor(hallmars_exp)
pht1 = pheatmap(mat = hallmark_cor,silent = T)

num_of_clusters = 4
clustering_distance = "euclidean"
myannotation = as.data.frame(cutree(pht1[["tree_row"]], k = num_of_clusters)) #split into k clusters
 
names(myannotation)[1] = "cluster"
  myannotation$cluster = as.factor(myannotation$cluster)
  
  palette1 <-brewer.pal(num_of_clusters, "Paired")

  names(palette1) = unique(myannotation$cluster)
  ann_colors = list (cluster = palette1)
  annotation = list(ann_colors = ann_colors, myannotation = myannotation)
  
  colors <- c(seq(-1,1,by=0.01))
  my_palette <- c("blue",colorRampPalette(colors = c("blue", "white", "red"))
                                                   (n = length(colors)-3), "red")


  print_tab(plt = 
                pheatmap(mat = hallmark_cor,annotation_col =  annotation[["myannotation"]], annotation_colors = annotation[["ann_colors"]], clustering_distance_rows = clustering_distance,clustering_distance_cols = clustering_distance,color = my_palette,breaks = colors,show_rownames = F,show_colnames = F)
            ,title = "genes expression heatmap")

genes expression heatmap

NA

chosen_clusters = c(2,3)
chosen_genes = annotation[["myannotation"]] %>% dplyr::filter(cluster %in% chosen_clusters ) %>% rownames() #take relevant genes
score <- apply(lung@assays$RNA@data[chosen_genes,],2,mean)
lung=AddMetaData(lung,score,hallmark_name)

print_tab(FeaturePlot(object = lung, features = hallmark_name),title = "Expression")

Expression

NA

cc_scores = FetchData(object = lung,vars = hallmark_name)

plt  =  ggplot(cc_scores, aes(x=HALLMARK_INTERFERON_GAMMA_RESPONSE)) +
  geom_density()+
   geom_vline(
    aes(xintercept=mean(HALLMARK_INTERFERON_GAMMA_RESPONSE) + sd(HALLMARK_INTERFERON_GAMMA_RESPONSE) ,color="1 SD"),
          linetype="dashed", size=1)+
    geom_vline(
    aes(xintercept=mean(HALLMARK_INTERFERON_GAMMA_RESPONSE) + 2*sd(HALLMARK_INTERFERON_GAMMA_RESPONSE) ,color="2 SD"),
          linetype="dashed", size=1)+
    geom_vline(
    aes(xintercept=mean(HALLMARK_INTERFERON_GAMMA_RESPONSE) ,color="mean"),
          linetype="dashed", size=1)

print_tab(plt = plt,title = "dist")

dist

NA

16.1 correlatino to nmf

cor_res = cor(lung$interferon_like,lung$HALLMARK_INTERFERON_GAMMA_RESPONSE)
print(paste("correlation of ifn program to HALLMARK_INTERFERON_GAMMA_RESPONSE:", cor_res))

[1] "correlation of ifn program to HALLMARK_INTERFERON_GAMMA_RESPONSE: 0.056908124816368"

17 APC signature

hallmark_name = "KEGG_ANTIGEN_PROCESSING_AND_PRESENTATION"
# genesets  =getGmt("./Data/h.all.v7.0.symbols.pluscc.gmt")
# geneIds= genesets[[hallmark_name]]@geneIds
geneIds = canonical_pathways %>% filter(gs_name == hallmark_name) %>% pull(gene_symbol) %>% intersect(rownames(lung))

hallmars_exp = FetchData(object = lung,vars = c(geneIds))
hallmars_exp = hallmars_exp[,colSums(hallmars_exp[])>0] #remove no expression genes
hallmark_cor = cor(hallmars_exp)
pht1 = pheatmap(mat = hallmark_cor,silent = T)

num_of_clusters = 4
clustering_distance = "euclidean"
myannotation = as.data.frame(cutree(pht1[["tree_row"]], k = num_of_clusters)) #split into k clusters
 
names(myannotation)[1] = "cluster"
  myannotation$cluster = as.factor(myannotation$cluster)
  
  palette1 <-brewer.pal(num_of_clusters, "Paired")

  names(palette1) = unique(myannotation$cluster)
  ann_colors = list (cluster = palette1)
  annotation = list(ann_colors = ann_colors, myannotation = myannotation)
  
  colors <- c(seq(-1,1,by=0.01))
  my_palette <- c("blue",colorRampPalette(colors = c("blue", "white", "red"))
                                                   (n = length(colors)-3), "red")


  print_tab(plt = 
                pheatmap(mat = hallmark_cor,annotation_col =  annotation[["myannotation"]], annotation_colors = annotation[["ann_colors"]], clustering_distance_rows = clustering_distance,clustering_distance_cols = clustering_distance,color = my_palette,breaks = colors,show_rownames = F,show_colnames = F)
            ,title = "genes expression heatmap")

genes expression heatmap

NA

NA

chosen_clusters = c(4)
chosen_genes = annotation[["myannotation"]] %>% dplyr::filter(cluster %in% chosen_clusters ) %>% rownames() #take relevant genes
score <- apply(lung@assays$RNA@data[chosen_genes,],2,mean)
lung=AddMetaData(lung,score,hallmark_name)

print_tab(FeaturePlot(object = lung, features = hallmark_name),title = "Expression")

Expression

NA

17.1 correlatino to nmf

cor_res = cor(lung$interferon_like,lung$KEGG_ANTIGEN_PROCESSING_AND_PRESENTATION)
print(paste("correlation of ifn program to KEGG_ANTIGEN_PROCESSING_AND_PRESENTATION:", cor_res))

[1] "correlation of ifn program to KEGG_ANTIGEN_PROCESSING_AND_PRESENTATION: 0.421455421412037"

18 INF like signature

hallmark_name = "IFN_like_genes"
i = 1
top_genes = patients_geps  %>%  arrange(desc(patients_geps[i])) #sort by score a
top = head(rownames(top_genes),100) #take top top_genes_num
chosen_genes = top %>% intersect(rownames(lung))
score <- apply(lung@assays$RNA@data[chosen_genes,],2,mean)
lung=AddMetaData(lung,score,hallmark_name)

print_tab(FeaturePlot(object = lung, features = hallmark_name),title = "Expression")

Expression

NA

18.1 correlatino to nmf

cor_res = cor(lung$interferon_like,lung$IFN_like_genes)
print(paste("correlation of ifn program to KEGG_ANTIGEN_PROCESSING_AND_PRESENTATION:", cor_res))

[1] "correlation of ifn program to KEGG_ANTIGEN_PROCESSING_AND_PRESENTATION: 0.417529710113569"

19 Signature regulation

metagenes_mean_compare(dataset = lung, time.point_var = "time.point",prefix = "patient",patient.ident_var = "patient.ident",pre_on = c("pre-treatment","on-treatment"),axis.text.x = 8,programs = c("HALLMARK_HYPOXIA","HALLMARK_INTERFERON_GAMMA_RESPONSE","KEGG_ANTIGEN_PROCESSING_AND_PRESENTATION","IFN_like_genes","HALLMARK_G2M_CHECKPOINT"))

HALLMARK_HYPOXIA per patient

HALLMARK_HYPOXIA

HALLMARK_INTERFERON_GAMMA_RESPONSE per patient

HALLMARK_INTERFERON_GAMMA_RESPONSE

KEGG_ANTIGEN_PROCESSING_AND_PRESENTATION per patient

KEGG_ANTIGEN_PROCESSING_AND_PRESENTATION

IFN_like_genes per patient

IFN_like_genes

HALLMARK_G2M_CHECKPOINT per patient

HALLMARK_G2M_CHECKPOINT

NA