microenv_analysis_clean

1 Parameters

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

2 functions

3 Data

acc = read_rds(file = data_to_read)

4 ACC primary cancer cells

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

by patient

print_tab(plt = DimPlot(acc,group.by = "lum_or_myo"),title = "cell type")

cell type

print_tab(plt = FeaturePlot(acc,features  = "luminal_over_myo"),title = "luminal_over_myo score")

luminal_over_myo score

NA

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

acc_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

5 Immune checkpoints

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

6 chemokines

6.1 Variable genes

print_tab(plt = FeaturePlot(object = acc,features = "CCL28"),title = "CCL28",subtitle_num = 3)

CCL28

print_tab(plt = FeaturePlot(object = acc,features = "CCL22"),title = "CCL22",subtitle_num = 3)

CCL22

print_tab(plt = FeaturePlot(object = acc,features = "CCL2"),title = "CCL2",subtitle_num = 3)

CCL2

NA

6.2 lum/myo DEG

CC genes in myo vs luminal deg:

cc_receptor_ligands = rownames(acc)[startsWith(x =rownames(acc), prefix = "CCR") | startsWith(x =rownames(acc), prefix = "CCL")]
deg_genes = cc_receptor_ligands[cc_receptor_ligands %in% rownames(acc_deg)]
acc_deg[deg_genes,]

NA

lumScore_vs_gene = FetchData(object = acc,vars = c("CCL28","lum_or_myo"))
lumScore_vs_gene = lumScore_vs_gene %>% dplyr::filter(lum_or_myo != "NA")
ccl28_plt = ggboxplot(lumScore_vs_gene, 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=12, label = 
                                                   paste("Mean=",mean(x) %>% round(2),
                                                         "\n 0 TPM cells =", (sum(x==0,na.rm = T)/length(x))%>% round(2))), geom="text") +
     theme(legend.position="none")

lumScore_vs_gene = FetchData(object = acc,vars = c("CCL22","lum_or_myo"))
lumScore_vs_gene = lumScore_vs_gene %>% dplyr::filter(lum_or_myo != "NA")
ccl22_plt = ggboxplot(lumScore_vs_gene, 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=12, label = 
                                                   paste("Mean=",mean(x) %>% round(2),
                                                         "\n 0 TPM cells =", (sum(x==0,na.rm = T)/length(x))%>% round(2))), geom="text") +
     theme(legend.position="none")

print_tab(plt = ccl28_plt,title = "CCL28 lum/myo",subtitle_num = 3)

CCL28 lum/myo

print_tab(plt = ccl22_plt,title = "CCL22 lum/myo",subtitle_num = 3)

CCL22 lum/myo

NA

7 CXC

7.1 Variable genes

print_tab(plt = FeaturePlot(object = acc,features = "CXCR4"),title = "receptors",subtitle_num = 3)

receptors

print_tab(plt = FeaturePlot(object = acc,features = c("CXCL1","CXCL2","CXCL3","CXCL14"))
,title = "receptors",subtitle_num = 3)

receptors

NA

7.2 lum/myo DEG

cxc_ligand = rownames(acc)[startsWith(x =rownames(acc), prefix = "CXCL") | startsWith(x =rownames(acc), prefix = "CXCR") ]
acc_deg[cxc_ligand[cxc_ligand %in% rownames(acc_deg)],]

lumScore_vs_gene = FetchData(object = acc,vars = c("CXCL2","lum_or_myo"))
lumScore_vs_gene = lumScore_vs_gene %>% dplyr::filter(lum_or_myo != "NA")
CXCL2_plt = ggboxplot(lumScore_vs_gene, x = "lum_or_myo", y = "CXCL2",
          palette = "jco",
          add = "jitter")+ 
  stat_compare_means(method = "wilcox.test",comparisons = list(c("luminal","myo")))+
  stat_summary(fun.data = function(x) data.frame(y=15, label = 
                                                   paste("Mean=",mean(x) %>% round(2),
                                                         "\n 0 TPM cells =", (sum(x==0,na.rm = T)/length(x))%>% round(2))), geom="text") +
     theme(legend.position="none")

lumScore_vs_gene = FetchData(object = acc,vars = c("CXCL3","lum_or_myo"))
lumScore_vs_gene = lumScore_vs_gene %>% dplyr::filter(lum_or_myo != "NA")
CXCL3_plt = ggboxplot(lumScore_vs_gene, x = "lum_or_myo", y = "CXCL3",
          palette = "jco",
          add = "jitter")+ 
  stat_compare_means(method = "wilcox.test",comparisons = list(c("luminal","myo")))+
  stat_summary(fun.data = function(x) data.frame(y=12, label = 
                                                   paste("Mean=",mean(x) %>% round(2),
                                                         "\n 0 TPM cells =", (sum(x==0,na.rm = T)/length(x))%>% round(2))), geom="text") +
     theme(legend.position="none")

lumScore_vs_gene = FetchData(object = acc,vars = c("CXCL14","lum_or_myo"))
lumScore_vs_gene = lumScore_vs_gene %>% dplyr::filter(lum_or_myo != "NA")
CXCL14_plt = ggboxplot(lumScore_vs_gene, 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) %>% round(2),
                                                         "\n 0 TPM cells =", (sum(x==0,na.rm = T)/length(x))%>% round(2))), geom="text") +
     theme(legend.position="none")


print_tab(plt = CXCL2_plt,title = "CXCL2 lum/myo",subtitle_num = 3)

CXCL2 lum/myo

print_tab(plt = CXCL3_plt,title = "CXCL3 lum/myo",subtitle_num = 3)

CXCL3 lum/myo

print_tab(plt = CXCL14_plt,title = "CXCL14 lum/myo",subtitle_num = 3)

CXCL14 lum/myo

lumScore_vs_gene = FetchData(object = acc,vars = c("CXCL1","lum_or_myo"))
lumScore_vs_gene = lumScore_vs_gene %>% dplyr::filter(lum_or_myo != "NA")
CXCL1_plt = ggboxplot(lumScore_vs_gene, x = "lum_or_myo", y = "CXCL1",
          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) %>% round(2),
                                                         "\n 0 TPM cells =", (sum(x==0,na.rm = T)/length(x))%>% round(2))), geom="text") +
     theme(legend.position="none")


print_tab(plt = CXCL1_plt,title = "CXCL1 lum/myo",subtitle_num = 3)

CXCL1 lum/myo

lumScore_vs_gene = FetchData(object = acc,vars = c("CXCL17","lum_or_myo"))
lumScore_vs_gene = lumScore_vs_gene %>% dplyr::filter(lum_or_myo != "NA")
CXCL17_plt = ggboxplot(lumScore_vs_gene, x = "lum_or_myo", y = "CXCL17",
          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) %>% round(2),
                                                         "\n 0 TPM cells =", (sum(x==0,na.rm = T)/length(x))%>% round(2))), geom="text") +
     theme(legend.position="none")


print_tab(plt = CXCL17_plt,title = "CXCL17 lum/myo",subtitle_num = 3)

CXCL17 lum/myo

NA

8 IL receptor

8.1 Variable genes

print_tab(plt = FeaturePlot(object = acc,features = c("IL17RB","IL17RD","IL17RE","IL17RC","IL17RA")),title = "IL17",subtitle_num = 3)

IL17

print_tab(plt = FeaturePlot(object = acc,features = c("IL1R2","IL1R1","IL1RAP","IL2RA")),title = "IL1-2",subtitle_num = 3)

IL1-2

print_tab(plt = FeaturePlot(object = acc,features = c( "IL10RB","IL11RA","IL27RA","IL13RA1")),title = "IL10-11-13-27",subtitle_num = 3)

IL10-11-13-27

NA

8.2 lum/myo DEG

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

for (gene in genes_in_deg) {
  lumScore_vs_gene = FetchData(object = acc,vars = c(gene,"lum_or_myo"))
lumScore_vs_gene = lumScore_vs_gene %>% dplyr::filter(lum_or_myo != "NA")
plt = 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=10, label = 
                                                   paste("Mean=",mean(x) %>% round(2),
                                                         "\n 0 TPM cells =", (sum(x==0,na.rm = T)/length(x))%>% round(2))), geom="text") +
     theme(legend.position="none")


print_tab(plt = plt,title = paste(gene,"lum/myo"),subtitle_num = 3)
}

IL17RD lum/myo

IL1RAP lum/myo

IL20RA lum/myo

IL2RA lum/myo

NA

9 HLA 1 genes

9.1 UMAP

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

print_tab(plt = FeaturePlot(object = acc,features = c("HLA-A","HLA-B","HLA-C")),title = "classical  HLA 1 genes",subtitle_num = 3)

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) {
  print_tab(plt =  FeaturePlot(object = acc,features = non_classical[i:(i+3)]),title = non_classical[i:(i+3)],subtitle_num = 3)

  i= i+4
  if (i > length(non_classical)) { exit = T}
}

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

NA

9.2 lum/myo DEG

genes_in_deg = hla_class_1[hla_class_1 %in% rownames(acc_deg)]
acc_deg[genes_in_deg,]

for (gene in genes_in_deg) {
  lumScore_vs_gene = FetchData(object = acc,vars = c(gene,"lum_or_myo"))
lumScore_vs_gene = lumScore_vs_gene %>% dplyr::filter(lum_or_myo != "NA")
plt = 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=max(x)+1, label = 
                                                   paste("Mean=",mean(x) %>% round(2),
                                                         "\n 0 TPM cells =", (sum(x==0,na.rm = T)/length(x))%>% round(2))), geom="text") +
     theme(legend.position="none")


print_tab(plt = plt,title = paste(gene,"lum/myo"),subtitle_num = 3)
}

HLA-F lum/myo

HLA-W lum/myo

HLA-L lum/myo

HLA-B lum/myo

NA

10 HLA 1 signature

correlation of non classical hla1

hla = FetchData(object = acc,vars = non_classical)

cor_res = cor(hla)
pheatmap(cor_res)

  geneIds = c("HLA-T","HLA-G","HLA-H","HLA-K","HLA-W","HLA-L","HLA-F","HLA-J")
  genes = paste(geneIds,collapse = ",")
  print("signature genes: " %>% paste(genes))

[1] "signature genes:  HLA-T,HLA-G,HLA-H,HLA-K,HLA-W,HLA-L,HLA-F,HLA-J"

  score <- apply(acc@assays$RNA@data[geneIds,],2,mean)
  acc=AddMetaData(acc,score,col.name = "HLA_signature")
  FeaturePlot(object = acc,features = "HLA_signature")

11 Other immune genes differentially expressed between lum/myo

intersect DEG 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.

acc_deg[intersect(rownames(acc_deg),immune_genes),]

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

12 Retinoic acid pathway

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(2)] #take relevant pathways
pathway_genes = genes %>% filter(gs_name == pathway) %>% pull("gene_symbol") 
print("pathway name:" %>% paste(retinoic_acid_pathways))

[1] "pathway name: REACTOME_SIGNALING_BY_RETINOIC_ACID"

print("genes in pathway:" )

[1] "genes in pathway:"

print(pathway_genes)

 [1] "ADH1A"   "ADH1C"   "ADH4"    "AKR1C3"  "ALDH1A1" "ALDH1A2" "ALDH1A3" "ALDH8A1" "CPT1A"   "CPT1B"   "CRABP1" 
[12] "CRABP2"  "CYP26A1" "CYP26B1" "CYP26C1" "DHRS3"   "DHRS4"   "DHRS4"   "DHRS9"   "DLAT"    "DLD"     "FABP5"  
[23] "PDHA1"   "PDHA2"   "PDHB"    "PDHX"    "PDK1"    "PDK2"    "PDK3"    "PDK4"    "PPARD"   "RARA"    "RARB"   
[34] "RARG"    "RDH10"   "RDH11"   "RDH13"   "RDH13"   "RDH13"   "RDH13"   "RDH13"   "RDH13"   "RDH13"   "RDH13"  
[45] "RDH13"   "RDH13"   "RDH14"   "RDH16"   "RDH5"    "RXRA"    "RXRB"    "RXRB"    "RXRB"    "RXRB"    "RXRB"   
[56] "RXRB"    "RXRG"    "SDR16C5"

12.1 DEG between lum and myo

Positive avg_log2FC indicate that the gene is more highly expressed in the Luminal cells, negative in myo cells

acc_deg[intersect(rownames(acc_deg),pathway_genes),]

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

ATRA pathway

• ALDH1A1, 2 and 3 utilise NAD+ as cofactor to oxidise all-trans-retinal (atRAL) to all-trans-retinoic acid (atRA).

• (FABP5) is a cytosolic, lipid-binding protein thought to bind all-trans-retinoic acid (atRA) and mediate its delivery to retinoic acid receptors (RARs)

• RDH10 is part of membrane bound oxidoreductases that reversibly catalyse the first and rate limiting step in retinoic acid biosynthesis, and use NAD+ as cofactor to the corresponding aldehyde all trans retinal (atRAL)

• atRA can significantly increase the expression of proteins involved in energy metabolism such as PDK via induction of PPARD (Wolf 2010)

This can support the paper from Piero Dalerba, since there are high in the Luminal cells.

• In constrast, DHRS3 may contribute to the reduction of all-trans-retinal (atRAL) to all-trans-retinol (atROL), (Haeseleer et al. 1998, Haeseleer et al. 2002, Kedishvili et al. 2002, Lin et al. 2001, Zhen et al. 2003, Belyaeva et al. 2008). DHRS3 is up in myo cells, thus supports the paper.