DE_LECs

Run Differential expression analysis within LEC subtypes

Standard DE analysis using pseudobulk DE methods (note this differs from the FGCZ default analysis). We want to identify differentially expressed genes within LEC cluster between all groups of skin, YUMMER and YUMM.

Preamble

library(scran)
library(dplyr)
library(tidyr)
library(ggplot2)
library(scuttle)
library(pheatmap)
library(patchwork)
library(stringr)
library(Seurat)
library(ggbeeswarm)
library(EnhancedVolcano)
library(SingleCellExperiment)
library(gridExtra)

Code

# volcano plot function as defined in https://github.com/HelenaLC/TLS-Silina/blob/main/code/geo-02-differential.qmd

.volcano <- \(df, title, fdr = 0.05, lfc = 1, select_lab = NULL) {
  EnhancedVolcano(df, 
    x = "logFC", y = "PValue",
    FCcutoff = lfc, pCutoff = fdr,
    selectLab = select_lab,
    pointSize = 1.7, raster = TRUE,
    title = title, subtitle = NULL,
    lab = rownames(df), labSize = 4, 
    drawConnectors = TRUE, widthConnectors = 0.5) +
  guides(col = guide_legend(override.aes = list(alpha = 1, size = 5))) +
  theme_bw(9) + theme(
    aspect.ratio = 1,
    legend.title = element_blank(),
    panel.grid.minor = element_blank())
}

.heatmap <- \(mtx, de_genes, cd, title, fdr = 0.05, n_lfc = 20) {
  top <- de_genes %>%
    filter(PValue < fdr) |> 
    slice_max(abs(logFC), n = n_lfc)
  mtx_sub <- log1p(mtx[rownames(top),])
  if (length(rownames(top)) < 2){
    return(print("No de genes"))
  }else{
  hm <- pheatmap(mtx_sub, 
      main = title, fontsize = 6,
      col = rev(hcl.colors(51, "RdBu")),
      scale = "row", show_colnames = FALSE, annotation_col = cd)
  hm
  }
}

Data object

seurat<- readRDS(file.path("/Users/thomarin/Documents/PhD/Tumor project/Sequencing experiment/August 2023 first analysis shallow sequencing/rds file/scData_LEC.rds"))

# correct condition assignment!!
seurat$cond <- seurat[[]] |> 
  mutate(
    cond = case_when(
       str_detect(Sample, "YUMM[0-9]") ~ "YUMM",
       str_detect(Sample, "YUMMER") ~ "YUMMER",
       str_detect(Sample, "Skin") ~ "skin"
    )
  ) |> select(cond)

# check assignment
table(seurat$Sample, seurat$cond)

                      
                       skin YUMM YUMMER
  SkinLECs_Leukocytes1  469    0      0
  SkinLECs_Leukocytes2  669    0      0
  TumorYUMM1_1A           0   17      0
  TumorYUMM1_1B           0   18      0
  TumorYUMM2_1A           0   57      0
  TumorYUMM2_1B           0   49      0
  TumorYUMM5_2A           0   44      0
  TumorYUMM5_2B           0   44      0
  TumorYUMM6_2A           0   58      0
  TumorYUMM6_2B           0   72      0
  TumorYUMMER3_1A         0    0     66
  TumorYUMMER3_1B         0    0     44
  TumorYUMMER4_1A         0    0     81
  TumorYUMMER4_1B         0    0     68
  TumorYUMMER7_2A         0    0     86
  TumorYUMMER7_2B         0    0    101
  TumorYUMMER8_2A         0    0    160
  TumorYUMMER8_2B         0    0    135

table(seurat$cond)

  skin   YUMM YUMMER 
  1138    359    741 

# switch to SingleCellExperiment object
sce <- as.SingleCellExperiment(seurat)

DE analysis

cond_de <- function(cond1, cond2){
  # subset sce
  sce_sub <- sce[,sce$cond %in% c(cond1, cond2)]
  # creating pseudobulks
  summed <- aggregateAcrossCells(sce_sub, 
                                 id=colData(sce_sub)[,c("ident", "Sample")])
  # filter min cells
  summed.filt <- summed[,summed$ncells >= 5]
  print(table(summed.filt$cond, summed.filt$ident))
  # de
  design <- model.matrix(~factor(summed.filt$cond), summed.filt$Samples)
  de.results <- pseudoBulkDGE(summed.filt, 
    label=summed.filt$ident,
    design=~factor(cond),
    coef=ncol(design),
    condition=summed.filt$cond 
)
}

skin_yumm <- cond_de("skin", "YUMM")

      
       0 1 2 3 4 5 6 7
  skin 2 2 2 2 2 0 2 1
  YUMM 6 4 5 3 3 6 1 5

#yummer_yumm <- cond_de("YUMMER", "YUMM")
yumm_yummer <- cond_de("YUMM", "YUMMER")

        
         0 1 2 3 4 5 6 7
  YUMM   6 4 5 3 3 6 1 5
  YUMMER 8 5 8 8 6 7 5 2

skin_yummer <- cond_de("skin","YUMMER")

        
         0 1 2 3 4 5 6 7
  skin   2 2 2 2 2 0 2 1
  YUMMER 8 5 8 8 6 7 5 2

#DE skin_yumm
is.de <- decideTestsPerLabel(skin_yumm, threshold=0.05)
summarizeTestsPerLabel(is.de)

  -1    0   1    NA
0 79 1216  56 11727
1 13  538  11 12516
2 64  773  26 12215
3 52  444  32 12550
4  8  623   8 12439
6 59 1663  42 11314
7 44  749 131 12154

#DE yumm_yummer
is.de <- decideTestsPerLabel(yumm_yummer, threshold=0.05)
summarizeTestsPerLabel(is.de)

  -1    0  1    NA
0  1 1155 15 11907
1  1  453  2 12622
2  6  999 11 12062
3  0 1985 13 11080
4  1  722  4 12351
5  3  529 10 12536
6  9 1403  1 11665
7  0  670  4 12404

#DE skin_yummer
is.de <- decideTestsPerLabel(skin_yummer, threshold=0.05)
summarizeTestsPerLabel(is.de)

   -1    0   1    NA
0 193 2059 168 10658
1  29 1262  26 11761
2 216 1755 152 10955
3 284 2497 297 10000
4  10  992  18 12058
6 114 1566  93 11305
7  44  456  94 12484

#sort

Plot results

skin vs yumm

#first condition (skin) = baseline, compared to second condition
plot_list <- lapply(names(skin_yumm), \(.) {
  df = skin_yumm[[.]]
  .volcano(df = df, title = ., fdr = 0.05, lfc = 1)
  }
)

wrap_plots(plot_list, nrow = 4) + 
  plot_layout(guides = "collect") &
  theme(legend.position = "top")

skin vs yumm selected genes

#first condition (skin) = baseline, compared to second condition
plot_list <- lapply(names(skin_yumm), \(.) {
  df = skin_yumm[[.]]
  .volcano(df = df, title = ., fdr = 0.05, lfc = 1, select_lab = c("Lyve1", "Cxcl9", "Cxcl10", "Sema3a", "CD112", "Ccl21a", "Reln", "Cxcl12", "H2-Aa", "H2-Ab1", "Cd74", "Ctss", "Ccr2", "Il10", "Alcam", "Esam", "Cd274", "Ifng","Tnfa"))
  #.volcano(df = df, title = ., fdr = 0.05, lfc = 1)
  }
)

wrap_plots(plot_list, nrow = 4) + 
  plot_layout(guides = "collect") &
  theme(legend.position = "top")

yumm vs yummer

plot_list <- lapply(names(yumm_yummer), \(.) {
  df = yumm_yummer[[.]]
  .volcano(df = df, title = ., fdr = 0.05, lfc = 1)
  }
)

wrap_plots(plot_list, nrow = 4) + 
  plot_layout(guides = "collect") &
  theme(legend.position = "top")

yumm vs yummer selected genes

#first condition (skin) = baseline, compared to second condition
plot_list <- lapply(names(yumm_yummer), \(.) {
  df = yumm_yummer[[.]]
  .volcano(df = df, title = ., fdr = 0.05, lfc = 1, select_lab = c("Lyve1", "Cxcl9", "Cxcl10", "Sema3a", "CD112", "Ccl21a", "Reln", "Cxcl12", "H2-Aa", "H2-Ab1", "Cd74", "Ctss", "Ccr2", "Il10", "Alcam", "Esam", "Cd274", "Ifng","Tnfa"))
  #.volcano(df = df, title = ., fdr = 0.05, lfc = 1)
  }
)

wrap_plots(plot_list, nrow = 4) + 
  plot_layout(guides = "collect") &
  theme(legend.position = "top")

skin vs yummer

plot_list <- lapply(names(skin_yummer), \(.) {
  df = skin_yummer[[.]]
  .volcano(df = df, title = ., fdr = 0.05, lfc = 1, select_lab = c("Lyve1"))
  .volcano(df = df, title = ., fdr = 0.05, lfc = 1)
  }
)

wrap_plots(plot_list, nrow = 4) + 
  plot_layout(guides = "collect") &
  theme(legend.position = "top")

skin vs yummer selected

#first condition (skin) = baseline, compared to second condition
plot_list <- lapply(names(skin_yummer), \(.) {
  df = skin_yummer[[.]]
  .volcano(df = df, title = ., fdr = 0.05, lfc = 1, select_lab = c("Lyve1", "Cxcl9", "Cxcl10", "Sema3a", "CD112", "Ccl21a", "Reln", "Cxcl12", "H2-Aa", "H2-Ab1", "Cd74", "Ctss", "Ccr2", "Il10", "Alcam", "Esam", "Cd274", "Ifng","Tnfa"))
  #.volcano(df = df, title = ., fdr = 0.05, lfc = 1)
  }
)

wrap_plots(plot_list, nrow = 4) + 
  plot_layout(guides = "collect") &
  theme(legend.position = "top")

Heatmap plots

skin vs yumm

for (. in names(skin_yumm)) {
    cat("####", ., "\n")
    sub <- subset(x = seurat, idents = .)
    sub <- subset(x = sub, subset = cond %in% c("YUMM", "skin"))
    mtx <- GetAssayData(object = sub, slot = 'data')
    top <- as.data.frame(skin_yumm[[.]]) %>%
    filter(PValue < 0.05) |> 
    slice_max(abs(logFC), n = 40)
    mtx_sub <- log1p(mtx[rownames(top),])
    cd <- data.frame(sub[[]] |> select('cond'))
    hm <- pheatmap(mtx_sub, 
        main = ., fontsize = 6,
        col = rev(hcl.colors(51, "RdBu")),
        scale = "row", 
        show_colnames = FALSE, 
        cluster_cols = FALSE,
        annotation_col = cd)
    print(hm); cat("\n\n")
}

#### 0 

#### 1 

#### 2 

#### 3 

#### 4 

#### 6 

#### 7 

yumm vs yummer

for (. in names(yumm_yummer)) {
    cat("####", ., "\n")
    sub <- subset(x = seurat, idents = .)
    sub <- subset(x = sub, subset = cond %in% c("YUMM", "YUMMER"))
    mtx <- GetAssayData(object = sub, slot = 'data')
    top <- as.data.frame(yumm_yummer[[.]]) %>%
    filter(PValue < 0.05) |> 
    slice_max(abs(logFC), n = 40)
    mtx_sub <- log1p(mtx[rownames(top),])
    cd <- data.frame(sub[[]] |> select('cond'))
    hm <- pheatmap(mtx_sub, 
        main = ., fontsize = 6,
        col = rev(hcl.colors(51, "RdBu")),
        scale = "row", 
        show_colnames = FALSE, 
        cluster_cols = FALSE,
        annotation_col = cd)
    print(hm); cat("\n\n")
}

#### 0 

#### 1 

#### 2 

#### 3 

#### 4 

#### 5 

#### 6 

#### 7 

skin vs yummer

for (. in names(skin_yummer)) {
    cat("####", ., "\n")
    sub <- subset(x = seurat, idents = .)
    sub <- subset(x = sub, subset = cond %in% c("YUMMER", "skin"))
    mtx <- GetAssayData(object = sub, slot = 'data')
    top <- as.data.frame(skin_yummer[[.]]) %>%
    filter(PValue < 0.05) |> 
    slice_max(abs(logFC), n = 40)
    mtx_sub <- log1p(mtx[rownames(top),])
    cd <- data.frame(sub[[]] |> select('cond'))
    hm <- pheatmap(mtx_sub, 
        main = ., fontsize = 6,
        col = rev(hcl.colors(51, "RdBu")),
        scale = "row", 
        show_colnames = FALSE, 
        cluster_cols = FALSE,
        annotation_col = cd)
    print(hm); cat("\n\n")
}

#### 0 

#### 1 

#### 2 

#### 3 

#### 4 

#### 6 

#### 7