nichenetR
Ehsan Razmara Monash Uni
15/03/2022
library(nichenetr)## Registered S3 method overwritten by 'spatstat.geom':
## method from
## print.boxx clilibrary(Seurat)## Attaching SeuratObjectlibrary(tidyverse)## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──## ✓ ggplot2 3.3.5 ✓ purrr 0.3.4
## ✓ tibble 3.1.6 ✓ dplyr 1.0.7
## ✓ tidyr 1.2.0 ✓ stringr 1.4.0
## ✓ readr 2.1.2 ✓ forcats 0.5.1## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()ligand_target_matrix <- readRDS("/mnt/nectar_volume/home/eraz0001/nichnetR/ligand_target_matrix.rds")
ligand_target_matrix[1:5 , 1:5]## CXCL1 CXCL2 CXCL3 CXCL5 PPBP
## A1BG 3.534343e-04 4.041324e-04 3.729920e-04 3.080640e-04 2.628388e-04
## A1BG-AS1 1.650894e-04 1.509213e-04 1.583594e-04 1.317253e-04 1.231819e-04
## A1CF 5.787175e-04 4.596295e-04 3.895907e-04 3.293275e-04 3.211944e-04
## A2M 6.027058e-04 5.996617e-04 5.164365e-04 4.517236e-04 4.590521e-04
## A2M-AS1 8.898724e-05 8.243341e-05 7.484018e-05 4.912514e-05 5.120439e-05lr_network <- readRDS("/mnt/nectar_volume/home/eraz0001/nichnetR/lr_network.rds")
weighted_networks <- readRDS("/mnt/nectar_volume/home/eraz0001/nichnetR/weighted_networks.rds")
weighted_networks_lr <- weighted_networks$lr_sig %>% inner_join(lr_network %>% distinct(from,to), by = c("from","to"))Mapping human genes to mouse ones
lr_network <- lr_network %>% mutate(from = convert_human_to_mouse_symbols(from), to = convert_human_to_mouse_symbols(to)) %>% drop_na()
colnames(ligand_target_matrix) <- ligand_target_matrix %>% colnames() %>% convert_human_to_mouse_symbols()
rownames(ligand_target_matrix) <- ligand_target_matrix %>% rownames() %>% convert_human_to_mouse_symbols()
ligand_target_matrix <- ligand_target_matrix %>% .[!is.na(rownames(ligand_target_matrix)), !is.na(colnames(ligand_target_matrix))]
weighted_networks_lr <- weighted_networks_lr %>% mutate(from = convert_human_to_mouse_symbols(from), to = convert_human_to_mouse_symbols(to)) %>% drop_na()so <- readRDS("/mnt/nectar_volume/home/eraz0001/KELLY 2020/E11.5/Final_15_clusters.RDS")
DimPlot(so , label = T)
so@meta.data$seurat_clusters %>% table()## .
## 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
## 243 199 194 190 185 177 174 148 145 133 133 118 114 102 40 38#subtract clusters
seuratObj <- subset(so , idents = c("Hoxd13+" , "DCh" , "Osteoblasts", "Osteocytes", "Endothelial Cells", "HtCh", "Chondrocytes"))## Receiver
receiver = c("DCh" , "Osteocytes", "Endothelial Cells", "Osteoblasts", "HtCh", "Chondrocytes")
expressed_genes_receiver = get_expressed_genes(receiver, seuratObj, pct = 0.10)
background_expressed_genes = expressed_genes_receiver %>% .[. %in% rownames(ligand_target_matrix)]
## Sender
sender_celltypes = c("Hoxd13+")
# lapply to get the expressed genes of every sender cell type separately here
list_expressed_genes_sender = sender_celltypes %>% unique() %>% lapply(get_expressed_genes, seuratObj, 0.10)
expressed_genes_sender = list_expressed_genes_sender %>% unlist() %>% unique()DE_table_receiver <- FindMarkers(object = seuratObj, ident.1 = c("DCh" , "Osteocytes", "Endothelial Cells", "Osteoblasts", "HtCh", "Chondrocytes"), ident.2 = "Hoxd13+" , min.pct = 0.10)
DE_table_receiver <- rownames_to_column(DE_table_receiver , "gene")
geneset_oi <- DE_table_receiver %>% filter(p_val_adj <= 0.05 & abs(avg_log2FC) >= 0.25) %>% pull(gene)
geneset_oi <- geneset_oi %>% .[. %in% rownames(ligand_target_matrix)]
dim(ligand_target_matrix)## [1] 17330 644dim(lr_network)## [1] 12156 4ligands = lr_network %>% pull(from) %>% unique()
receptors = lr_network %>% pull(to) %>% unique()expressed_ligands = intersect(ligands,expressed_genes_sender)
expressed_receptors = intersect(receptors,expressed_genes_receiver)
potential_ligands = lr_network %>% filter(from %in% expressed_ligands & to %in% expressed_receptors) %>% pull(from) %>% unique()ligand_activities <- predict_ligand_activities(geneset = geneset_oi, background_expressed_genes = background_expressed_genes, ligand_target_matrix = ligand_target_matrix, potential_ligands = potential_ligands)ligand_activities <- ligand_activities %>% arrange(-pearson) %>% mutate(rank = rank(sort(pearson , decreasing = F)))DimPlot(seuratObj , label = T)
best_upstream_ligands = ligand_activities %>% top_n(20, pearson) %>% pull(test_ligand) %>% unique()
DotPlot(seuratObj, features = best_upstream_ligands %>% rev(), cols = "RdYlBu") + RotatedAxis()
active_ligand_target_links_df = best_upstream_ligands %>%
lapply(get_weighted_ligand_target_links,geneset = geneset_oi, ligand_target_matrix = ligand_target_matrix,
n = 250) %>% bind_rows() %>% drop_na()
active_ligand_target_links = prepare_ligand_target_visualization(ligand_target_df = active_ligand_target_links_df,
ligand_target_matrix = ligand_target_matrix, cutoff = 0.33)dimnames(active_ligand_target_links)## [[1]]
## [1] "Cald1" "Ahnak" "Jun" "Fos" "Lef1" "Ets1"
## [7] "Zfp36l1" "Etv4" "Gapdh" "Hoxc10" "Vcan" "Klf6"
## [13] "Icam1" "Greb1" "Ltbp1" "Zbtb20" "Anxa2" "Wnt5a"
## [19] "Bcl11a" "Hoxa13" "Ranbp1" "Pfn1" "Rcc1" "Pmaip1"
## [25] "Cav1" "Hoxd13" "Tfap2a" "Gja1" "Ldha" "Eno1"
## [31] "Slc3a2" "Ncl" "Ezh2" "E2f1" "Mcm7" "Gsc"
## [37] "Sfrp1" "Maz" "Bmp4" "Cxxc5" "Tgfb2" "Tpm1"
## [43] "Sdc1" "Stmn1" "Gnas" "Osr2" "Lrig3" "Emp1"
## [49] "Cdc45" "Mdk" "Foxc1" "Evx1" "Pcna" "Ddah1"
## [55] "Hsp90ab1" "Msx1" "Igf2" "Prmt1" "Rnaseh2a" "Ptma"
## [61] "Trp53" "Ubc" "Rbbp4" "Rps2" "Set" "Nolc1"
## [67] "Tfap4" "Hnrnpa1" "Slc16a3" "Cdc6" "Hnrnpf" "Srsf1"
## [73] "Hspe1" "Atf4" "Mmp2" "Zwint" "Hat1" "Tomm20"
## [79] "C1qbp" "Fbl" "Eif1ax" "Cops3" "Slc25a5" "Phb2"
## [85] "Crmp1" "Eif3b" "H2afz" "Ahcy" "Eif5a" "Nme1"
## [91] "Anp32a" "Umps" "Cox5a" "Cycs" "Psat1" "Cdca7"
## [97] "Ebna1bp2" "Psma7" "Grwd1" "Erh" "Psma1" "Fam60a"
## [103] "Mef2c"
##
## [[2]]
## [1] "Fbrs" "Ptdss1" "Robo2" "Hsp90b1" "Col4a1" "Yars"
## [7] "Flrt3" "Rps19" "Adam9" "Serpine2" "Efna1" "Efnb1"
## [13] "Ntn3" "Arf1" "Efnb2" "Copa" "Plxnb2" "Vcam1"
## [19] "Nenf" "Pcdh7"order_ligands = intersect(best_upstream_ligands, colnames(active_ligand_target_links)) %>% rev() %>% make.names()
order_targets = active_ligand_target_links_df$target %>% unique() %>% intersect(rownames(active_ligand_target_links)) %>% make.names()
rownames(active_ligand_target_links) = rownames(active_ligand_target_links) %>% make.names() # make.names() for heatmap visualization of genes like H2-T23
colnames(active_ligand_target_links) = colnames(active_ligand_target_links) %>% make.names() # make.names() for heatmap visualization of genes like H2-T23
dimnames(active_ligand_target_links)## [[1]]
## [1] "Cald1" "Ahnak" "Jun" "Fos" "Lef1" "Ets1"
## [7] "Zfp36l1" "Etv4" "Gapdh" "Hoxc10" "Vcan" "Klf6"
## [13] "Icam1" "Greb1" "Ltbp1" "Zbtb20" "Anxa2" "Wnt5a"
## [19] "Bcl11a" "Hoxa13" "Ranbp1" "Pfn1" "Rcc1" "Pmaip1"
## [25] "Cav1" "Hoxd13" "Tfap2a" "Gja1" "Ldha" "Eno1"
## [31] "Slc3a2" "Ncl" "Ezh2" "E2f1" "Mcm7" "Gsc"
## [37] "Sfrp1" "Maz" "Bmp4" "Cxxc5" "Tgfb2" "Tpm1"
## [43] "Sdc1" "Stmn1" "Gnas" "Osr2" "Lrig3" "Emp1"
## [49] "Cdc45" "Mdk" "Foxc1" "Evx1" "Pcna" "Ddah1"
## [55] "Hsp90ab1" "Msx1" "Igf2" "Prmt1" "Rnaseh2a" "Ptma"
## [61] "Trp53" "Ubc" "Rbbp4" "Rps2" "Set" "Nolc1"
## [67] "Tfap4" "Hnrnpa1" "Slc16a3" "Cdc6" "Hnrnpf" "Srsf1"
## [73] "Hspe1" "Atf4" "Mmp2" "Zwint" "Hat1" "Tomm20"
## [79] "C1qbp" "Fbl" "Eif1ax" "Cops3" "Slc25a5" "Phb2"
## [85] "Crmp1" "Eif3b" "H2afz" "Ahcy" "Eif5a" "Nme1"
## [91] "Anp32a" "Umps" "Cox5a" "Cycs" "Psat1" "Cdca7"
## [97] "Ebna1bp2" "Psma7" "Grwd1" "Erh" "Psma1" "Fam60a"
## [103] "Mef2c"
##
## [[2]]
## [1] "Fbrs" "Ptdss1" "Robo2" "Hsp90b1" "Col4a1" "Yars"
## [7] "Flrt3" "Rps19" "Adam9" "Serpine2" "Efna1" "Efnb1"
## [13] "Ntn3" "Arf1" "Efnb2" "Copa" "Plxnb2" "Vcam1"
## [19] "Nenf" "Pcdh7"vis_ligand_target = active_ligand_target_links[order_targets,order_ligands] %>% t()p_ligand_target_network = vis_ligand_target %>%
make_heatmap_ggplot("Prioritized ligands","Predicted target genes", color = "red",
legend_position = "top", x_axis_position = "top",legend_title = "Regulatory potential") +
theme(axis.text.x = element_text(face = "italic")) +
scale_fill_gradient2(low = "whitesmoke", high = "red", breaks = c(0,0.0045,0.0090))## Scale for 'fill' is already present. Adding another scale for 'fill', which
## will replace the existing scale.p_ligand_target_network
lr_network_top = lr_network %>% filter(from %in% best_upstream_ligands & to %in% expressed_receptors) %>% distinct(from,to)
best_upstream_receptors = lr_network_top %>% pull(to) %>% unique()
lr_network_top_df_large = weighted_networks_lr %>% filter(from %in% best_upstream_ligands & to %in% best_upstream_receptors)
lr_network_top_df = lr_network_top_df_large %>% spread("from","weight",fill = 0)
lr_network_top_matrix = lr_network_top_df %>% dplyr::select(-to) %>% as.matrix() %>% magrittr::set_rownames(lr_network_top_df$to)
dist_receptors = dist(lr_network_top_matrix, method = "binary")
hclust_receptors = hclust(dist_receptors, method = "ward.D2")
order_receptors = hclust_receptors$labels[hclust_receptors$order]
class(order_receptors)## [1] "character"dist_ligands = dist(lr_network_top_matrix %>% t(), method = "binary")
hclust_ligands = hclust(dist_ligands, method = "ward.D2")
order_ligands_receptor = hclust_ligands$labels[hclust_ligands$order]
order_receptors = order_receptors %>% intersect(rownames(lr_network_top_matrix))
order_ligands_receptor = order_ligands_receptor %>% intersect(colnames(lr_network_top_matrix))
vis_ligand_receptor_network <- lr_network_top_matrix[order_receptors, order_ligands_receptor]
rownames(vis_ligand_receptor_network) <- order_receptors %>% make.names()
colnames(vis_ligand_receptor_network) <- order_ligands_receptor %>% make.names()p_ligand_receptor_network <- vis_ligand_receptor_network %>% t() %>% make_heatmap_ggplot("Ligands","Receptors", color = "mediumvioletred", x_axis_position = "top",legend_title = "Prior interaction potential")
p_ligand_receptor_network
ligand_pearson_matrix = ligand_activities %>% select(pearson) %>% as.matrix() %>% magrittr::set_rownames(ligand_activities$test_ligand)
vis_ligand_pearson = ligand_pearson_matrix[order_ligands, ] %>% as.matrix(ncol = 1) %>% magrittr::set_colnames("Pearson")
p_ligand_pearson = vis_ligand_pearson %>% make_heatmap_ggplot("Prioritized CAF-ligands","Ligand activity", color = "darkorange",legend_position = "top", x_axis_position = "top", legend_title = "Pearson correlation coefficient\ntarget gene prediction ability)")
p_ligand_pearson