heart_revremod_paper/analysis/singlenuc/nuc_13.4_nichenet.R at main · saezlab/heart_revremod_paper · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
# Copyright (c) [2025] [Ines Rivero Garcia]
# ines.rivero@uni-heidelberg.de

#' In this script we predict target genes in receiver cell from the candidate
#' ligands suggested by LIANA+

suppressPackageStartupMessages({
  library(nichenetr)
  library(Seurat)
  #library(sceasy)
  #library(reticulate)
  library(ggplot2)
  library(tidyr)
  library(dplyr)
  library(purrr)
  library(tibble)
  #library(edgeR)
})

setwd("/mnt/sds-hd/sd22b002/projects/ines/heart_revremod")

# Load adata object as a seurat object
seuratObj = schard::h5ad2seurat('data/CRC_integrated_adata_ann.h5ad')
seuratObj

# Load NicheNet networks
lr_network <- readRDS("data/lr_network_mouse_21122021.rds")
ligand_target_matrix <- readRDS("data/ligand_target_matrix_nsga2r_final_mouse.rds")
weighted_networks <- readRDS("data/weighted_networks_nsga2r_final_mouse.rds")
lr_network <- lr_network %>% distinct(from, to)

# Set conditions and filter seurat object accordingly
condition_oi <- "DB"
condition_reference <- "ORAB"
grouping <- "ConditionSimplified"

Idents(seuratObj) <- "cell_type"
seuratObj
DimPlot(seuratObj, reduction = "Xumap_", split.by = grouping)

# Load differentially expressed genes
degs <- read.csv("results/degs/edgeR_res_condition-simplified.csv")
head(degs)

# Load Liana+ results
liana <- read.csv("results/liana/liana_condition-simplified_DBvORAB_significant_results_candidates.csv", header = TRUE, row.names = 1) %>%
    select(!X)
head(liana)

# Get receiver cells from liana
receivers <- liana %>%
  pull(target) %>%
  unique()
receivers

# Iterate thru receivers to get gene sets, run nichenet and extract ligands
for(cell in receivers){
  print(paste0("Running Nichenet for ", cell, " as receiver."))

  # Get ligands acting on this cell type from liana
  ligands <- liana %>%
    filter(target == cell) %>%
    pull(ligand_complex) %>%
    unique()
  print(paste0("... ", length(ligands), " ligands identified by Liana+ act on ", cell))

  # Define background genes
  expressed_genes_receiver <- get_expressed_genes(cell, seuratObj, pct = 0.1)
  background_genes <- expressed_genes_receiver %>% .[. %in% rownames(ligand_target_matrix)]
  print(paste0("... Number of background genes in receiver: ", length(background_genes)))

  # Get overexpressed genes in receiver
  geneset_oi <- degs %>%
    filter(cell_type == cell & FDR < 0.05) %>%
    pull(gene)
  print(paste0("... Number of differentially expressed genes in receiver: ", length(geneset_oi)))

  # Run Nichenet
  nichenet_activities <- predict_ligand_activities(
    geneset = geneset_oi,
    background_expressed_genes = background_genes,
    ligand_target_matrix = ligand_target_matrix,
    potential_ligands = ligands
  )
  write.table(nichenet_activities,
              paste0("results/nichenet/nichenet_condition-simplified_", condition_oi, "vs", condition_reference, "_signicant_receiver", cell, ".csv"),
              col.names = TRUE,
              row.names = FALSE,
              quote = FALSE,
              sep = ",")

  # Extract potential targets for each ligand
  active_ligand_target_links_df <- ligands %>%
    lapply(get_weighted_ligand_target_links,
           geneset = geneset_oi,
           ligand_target_matrix = ligand_target_matrix,
           n = 200) %>%
    bind_rows() %>%
    drop_na()

  # Add target logFC and FDR
  for(i in 1:nrow(active_ligand_target_links_df)){
    gn <- as.character(active_ligand_target_links_df[i, "target"])
    active_ligand_target_links_df[i, "logFC"] <- degs %>% filter(cell_type == cell & gene == gn & contrast == paste0(condition_oi, "v", condition_reference)) %>% pull(logFC)
    active_ligand_target_links_df[i, "FDR"] <- degs %>% filter(cell_type == cell & gene == gn  & contrast == paste0(condition_oi, "v", condition_reference)) %>% pull(FDR)
  }

  write.table(active_ligand_target_links_df,
              paste0("results/nichenet/nichenet_condition-simplified_", condition_oi, "vs", condition_reference, "_significant_receiver", cell, "_targetgenes.csv"),
              col.names = TRUE,
              row.names = FALSE,
              quote = FALSE,
              sep = ",")
}