microenv_analysis

Parameters

suffix = ""
data_to_read = "./Data/acc_cancer_no146_primaryonly15k_cancercells.rds"

functions

Data

acc = read_rds(file = data_to_read)

DimPlot(acc,group.by = "patient.ident")

Immune checkpoints

FeaturePlot(object = acc,features = c("CD247","PDCD1LG2","CTLA4"))

cc_receptor = rownames(acc)[startsWith(x =rownames(acc), prefix = "CCR")]
cc_ligand = rownames(acc)[startsWith(x =rownames(acc), prefix = "CCL")]

CC receptors

i=1
exit = F
while (!exit) {
  
  cat("### ",cc_receptor[i:(i+3)]," \n")
    print(
    FeaturePlot(object = acc,features = cc_receptor[i:(i+3)])
  )
  i= i+4
  if (i > length(cc_receptor)) { exit = T}
   cat(' \n\n')

  plot.new()
  dev.off()

}

CCR4 CCR9 CCR3 CCR1

Warning in FeaturePlot(object = acc, features = cc_receptor[i:(i + 3)]) : All cells have the same value (0) of CCR4.

CCR2 CCR5 CCRL2 CCR6

Warning in FeaturePlot(object = acc, features = cc_receptor[i:(i + 3)]) : All cells have the same value (0) of CCR2. Warning in FeaturePlot(object = acc, features = cc_receptor[i:(i + 3)]) : All cells have the same value (0) of CCR5. Warning in FeaturePlot(object = acc, features = cc_receptor[i:(i + 3)]) : All cells have the same value (0) of CCR6.

CCR7 CCR10 NA NA

Warning in FetchData.Seurat(object = object, vars = c(dims, “ident”, features), : The following requested variables were not found: NA

CC ligands

i=1
exit = F
while (!exit) {
  
  cat("### ",cc_ligand[i:(i+3)]," \n")
    print(
    FeaturePlot(object = acc,features = cc_ligand[i:(i+3)])
  )
  i= i+4
  if (i > length(cc_ligand)) { exit = T}
   cat(' \n\n')

  plot.new()
  dev.off()

}

CCL20 CCL28 CCL26 CCL19

Warning in FeaturePlot(object = acc, features = cc_ligand[i:(i + 3)]) : All cells have the same value (0) of CCL26.

CCL21 CCL22 CCL17 CCL2

CCL7 CCL11 CCL8 CCL13

Warning in FeaturePlot(object = acc, features = cc_ligand[i:(i + 3)]) : All cells have the same value (0) of CCL11. Warning in FeaturePlot(object = acc, features = cc_ligand[i:(i + 3)]) : All cells have the same value (0) of CCL13.

CCL5 CCL16 CCL14 CCL15

Warning in FeaturePlot(object = acc, features = cc_ligand[i:(i + 3)]) : All cells have the same value (0) of CCL16. Warning in FeaturePlot(object = acc, features = cc_ligand[i:(i + 3)]) : All cells have the same value (0) of CCL15.

CCL23 CCL18 CCL3 CCL4

Warning in FeaturePlot(object = acc, features = cc_ligand[i:(i + 3)]) : All cells have the same value (0) of CCL23.

CCL3L3 CCL4L2 CCL25 NA

Warning in FetchData.Seurat(object = object, vars = c(dims, “ident”, features), : The following requested variables were not found: NA Warning in FeaturePlot(object = acc, features = cc_ligand[i:(i + 3)]) : All cells have the same value (0) of CCL3L3.

lumScore = FetchData(object = acc,vars = c("luminal_over_myo"))
lumScore = lumScore %>% mutate(lum_or_myo = case_when( 
  luminal_over_myo > 1  ~ "luminal",
  luminal_over_myo < (-1)  ~ "myo",
  TRUE ~ "NA"))

acc = AddMetaData(object = acc,metadata = lumScore)
DimPlot(acc,group.by = "lum_or_myo")

lumScore_vs_program = FetchData(object = acc,vars = c("CCL28","lum_or_myo"))
lumScore_vs_program = lumScore_vs_program %>% dplyr::filter(lum_or_myo != "NA")
ggboxplot(lumScore_vs_program, x = "lum_or_myo", y = "CCL28",
          palette = "jco",
          add = "jitter")+ stat_compare_means(method = "wilcox.test",comparisons = list(c("luminal","myo")))+
  stat_summary(fun.data = function(x) data.frame(y=11, label = paste("Mean=",mean(x))), geom="text") +
     theme(legend.position="none")

lumScore_vs_program = FetchData(object = acc,vars = c("CCL22","lum_or_myo"))
lumScore_vs_program = lumScore_vs_program %>% dplyr::filter(lum_or_myo != "NA")
ggboxplot(lumScore_vs_program, x = "lum_or_myo", y = "CCL22",
          palette = "jco",
          add = "jitter")+ stat_compare_means(method = "wilcox.test",comparisons = list(c("luminal","myo")))+ 
     stat_summary(fun.data = function(x) data.frame(y=11, label = paste("Mean=",mean(x))), geom="text") +
     theme(legend.position="none")

cxc_receptor = rownames(acc)[startsWith(x =rownames(acc), prefix = "CXCR")]
cxc_ligand = rownames(acc)[startsWith(x =rownames(acc), prefix = "CXCL")]

CXC receptors

i=1
exit = F
while (!exit) {
  
  cat("### ",cxc_receptor[i:(i+3)]," \n")
    print(
    FeaturePlot(object = acc,features = cxc_receptor[i:(i+3)])
  )
  i= i+4
  if (i > length(cxc_receptor)) { exit = T}
   cat(' \n\n')

  plot.new()
  dev.off()

}

CXCR4 CXCR2 CXCR6 CXCR5

CXCR3 NA NA NA

Warning in FetchData.Seurat(object = object, vars = c(dims, “ident”, features), : The following requested variables were not found: NA

CXC ligands

i=1
exit = F
while (!exit) {
  
  cat("### ",cxc_ligand[i:(i+3)]," \n")
    print(
    FeaturePlot(object = acc,features = cxc_ligand[i:(i+3)])
  )
  i= i+4
  if (i > length(cxc_ligand)) { exit = T}
   cat(' \n\n')

  plot.new()
  dev.off()

}

CXCL8 CXCL6 CXCL1 CXCL5

CXCL3 CXCL2 CXCL9 CXCL10

CXCL11 CXCL13 CXCL14 CXCL12

Warning in FeaturePlot(object = acc, features = cxc_ligand[i:(i + 3)]) : All cells have the same value (0) of CXCL13.

CXCL16 CXCL17 NA NA

Warning in FetchData.Seurat(object = object, vars = c(dims, “ident”, features), : The following requested variables were not found: NA

lumScore_vs_program = FetchData(object = acc,vars = c("CXCL14","lum_or_myo"))
lumScore_vs_program = lumScore_vs_program %>% dplyr::filter(lum_or_myo != "NA")
ggboxplot(lumScore_vs_program, x = "lum_or_myo", y = "CXCL14",
          palette = "jco",
          add = "jitter")+ stat_compare_means(method = "wilcox.test",comparisons = list(c("luminal","myo")))+
  stat_summary(fun.data = function(x) data.frame(y=16, label = paste("Mean=",mean(x))), geom="text") +
     theme(legend.position="none")

NA
NA
NA

IL receptors genes

il_receptors_genes = rownames(acc)[grepl("^IL\\d*R.*", rownames(acc))]

i=1
exit = F
while (!exit) {
  
  cat("### ",il_receptors_genes[i:(i+3)]," \n")
    print(
    FeaturePlot(object = acc,features = il_receptors_genes[i:(i+3)])
  )
  i= i+4
  if (i > length(il_receptors_genes)) { exit = T}
   cat(' \n\n')

  plot.new()
  dev.off()

}

IL22RA1 IL23R IL12RB2 IL6R

IL1R2 IL1R1 IL1RL2 IL1RL1

IL18R1 IL18RAP IL36RN IL1RN

IL8RBP IL5RA IL17RE IL17RC

IL17RB IL17RD IL20RB IL1RAP

IL7R IL31RA IL20RA IL22RA2

Warning in FeaturePlot(object = acc, features = il_receptors_genes[i:(i + : All cells have the same value (0) of IL22RA2.

IL11RA IL9RP2 IL15RA IL2RA

IL10RA IL4R IL21R IL9RP4

Warning in FeaturePlot(object = acc, features = il_receptors_genes[i:(i + : All cells have the same value (0) of IL21R.

IL27RA IL12RB1 IL10RB IL17RA

IL2RB IL17REL IL3RA IL1RAPL1

IL2RG IL1RAPL2 IL13RA2 IL13RA1

Warning in FeaturePlot(object = acc, features = il_receptors_genes[i:(i + : All cells have the same value (0) of IL13RA2.

IL9R NA NA NA

Warning in FetchData.Seurat(object = object, vars = c(dims, “ident”, features), : The following requested variables were not found: NA

Classical HLA 1 genes

hla_genes = rownames(acc)[startsWith(x =rownames(acc), prefix = "HLA")]
hla_class_2 = hla_genes[startsWith(x =hla_genes, prefix = "HLA-D")]
hla_class_1 = hla_genes[!startsWith(x =hla_genes, prefix = "HLA-D")]

FeaturePlot(object = acc,features = c("HLA-A","HLA-B","HLA-C"))

Non-classical HLA 1 genes

i=1
exit = F
non_classical  = hla_class_1 [!hla_class_1 %in% c("HLA-A","HLA-B","HLA-C")] %>% sort()
while (!exit) {
  
  cat("### ",non_classical[i:(i+3)]," \n")
    print(
    FeaturePlot(object = acc,features = non_classical[i:(i+3)])
  )
  i= i+4
  if (i > length(non_classical)) { exit = T}
   cat(' \n\n')

  plot.new()
  dev.off()

}

HLA-E HLA-F HLA-G HLA-H

HLA-J HLA-K HLA-L HLA-N

HLA-P HLA-S HLA-T HLA-U

HLA-V HLA-W HLA-Z NA

Warning in FetchData.Seurat(object = object, vars = c(dims, “ident”, features), : The following requested variables were not found: NA

acc= AddModuleScore(object = acc,features = list(c("HLA-T","HLA-G","HLA-H","HLA-K","HLA-W","HLA-L","HLA-F","HLA-J")),name = "hla_score")

FeaturePlot(object = acc,features = "hla_score1")

HLA 2 genes

i=1

exit = F
while (!exit) {
    cat("### ",hla_class_2[i:(i+3)]," \n")

    print(
    FeaturePlot(object = acc,features = hla_class_2[i:(i+3)])
  )
 i= i+4
  if (i > length(hla_class_2)) { exit = T}
   cat(' \n\n')

  plot.new()
  dev.off()
}

IL17RD

    FeaturePlot(object = acc,features = "IL17RD")

acc = SetIdent(object = acc,value = "lum_or_myo")
deg = FindMarkers(object = acc,ident.1 = "luminal",ident.2 = "myo",densify = T,logfc.threshold = 0)

deg  =deg %>% mutate(fdr = p.adjust(p_val,method = "fdr"))%>% #add fdr
    filter((avg_log2FC>1 & fdr<0.1) | (avg_log2FC< (-1) & fdr<0.1))  #filter significant

deg

str_detect(rownames(deg), "HLA") %>% table()

.
FALSE  TRUE 
 2681     7 

 rownames(deg)[which(str_detect(rownames(deg), "HLA"))]

[1] "HLA-B"    "HLA-L"    "HLA-DRB5" "HLA-DRB1" "HLA-W"    "HLA-DRB6" "HLA-F"   

deg_hla = rownames(deg)[which(str_detect(rownames(deg), "HLA"))]

deg %>% dplyr::filter(row.names(deg) %in% deg_hla)

lumScore_vs_program = FetchData(object = acc,vars = c("HLA-B","lum_or_myo"))
lumScore_vs_program = lumScore_vs_program %>% dplyr::filter(lum_or_myo != "NA")
ggboxplot(lumScore_vs_program, x = "lum_or_myo", y = "HLA-B",
          palette = "jco",
          add = "jitter")+ stat_compare_means(method = "wilcox.test",comparisons = list(c("luminal","myo")))+
  stat_summary(fun.data = function(x) data.frame(y=16, label = paste("Mean=",mean(x))), geom="text") +
     theme(legend.position="none")

NA
NA
NA

intersect with GO “immune response GO:0006955”

# ensembl = useMart("ensembl",dataset="hsapiens_gene_ensembl") #uses human ensembl annotations
# 
#  immune_genes <- getBM(attributes=c('hgnc_symbol'),
#        filters = 'go', values = 'GO:0002376', mart = ensembl)

immune_genes = fread(input = "./Data/GO_term_summary_20230218_165937.txt",sep = "\t", select = 2) %>% pull(1)

Warning in fread(input = "./Data/GO_term_summary_20230218_165937.txt", sep = "\t",  :
  Detected 11 column names but the data has 12 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.

intersect(rownames(deg),immune_genes)

[1] "C3"

FeaturePlot(object = acc_cancer,features = "C3")

FeaturePlot(object = acc_cancer,features = "RARA")

library(msigdbr)

genes = msigdbr(species = "Homo sapiens", category = "C2") %>% filter(gs_subcat != "CGP" & grepl('RETINOIC', gs_name))

retinoic_acid_pathways = unique(genes$gs_name)
retinoic_acid_pathways = retinoic_acid_pathways[c(1,2)] #take relevant pathways


for (pathway in retinoic_acid_pathways) {
  pathway_genes = genes %>% filter(gs_name == pathway) %>% pull("gene_symbol")
  print(
    deg %>% dplyr::filter(row.names(deg) %in% intersect(rownames(deg),pathway_genes))
  )
  sig_genes = deg %>% dplyr::filter(row.names(deg) %in% intersect(rownames(deg),pathway_genes)) %>% rownames()

#   for (gene  in sig_genes) {
#     lumScore_vs_gene = FetchData(object = acc,vars = c(gene,"lum_or_myo"))
# lumScore_vs_gene = lumScore_vs_gene %>% dplyr::filter(lum_or_myo != "NA")
# p = ggboxplot(lumScore_vs_gene, x = "lum_or_myo", y = gene,
#           palette = "jco",
#           add = "jitter")+ stat_compare_means(method = "wilcox.test",comparisons = list(c("luminal","myo")))+
#   stat_summary(fun.data = function(x) data.frame(y=16, label = paste("Mean=",mean(x))), geom="text") +
#      theme(legend.position="none")
# 
# print(p)
#   }
}

NA
NA
NA
NA

https://reactome.org/PathwayBrowser/#/R-HSA-5362517

for (pathway in retinoic_acid_pathways) {
  pathway_genes = genes %>% filter(gs_name == pathway) %>% pull("gene_symbol")
  geneIds = intersect(pathway_genes, VariableFeatures(acc)[1:10000])
  score <- apply(acc@assays$RNA@data[geneIds,],2,mean)
  acc=AddMetaData(acc,score,col.name = pathway)
}