Skip to content

Commit

Permalink
Add sender-ligand DE calculation to the seurat wrapper function
Browse files Browse the repository at this point in the history
  • Loading branch information
@browaeysrobin
browaeysrobin committed Jan 11, 2022
1 parent 241b622 commit 0bf85d6
Show file tree
Hide file tree
Showing 20 changed files with 161 additions and 123 deletions.
99 changes: 93 additions & 6 deletions R/application_prediction.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -818,8 +818,9 @@ nichenet_seuratobj_aggregate = function(receiver, seurat_obj, condition_colname,
## sender
if (length(sender) == 1){
if (sender == "all"){
list_expressed_genes_sender = Idents(seurat_obj) %>% unique() %>% lapply(get_expressed_genes, seurat_obj, expression_pct, assay_oi)
names(list_expressed_genes_sender) = Idents(seurat_obj) %>% unique()
sender_celltypes = Idents(seurat_obj) %>% levels()
list_expressed_genes_sender = sender_celltypes %>% lapply(get_expressed_genes, seurat_obj, expression_pct, assay_oi)
names(list_expressed_genes_sender) = sender_celltypes
expressed_genes_sender = list_expressed_genes_sender %>% unlist() %>% unique()

} else if (sender == "undefined") {
Expand Down Expand Up @@ -1031,6 +1032,64 @@ nichenet_seuratobj_aggregate = function(receiver, seurat_obj, condition_colname,
lr_network_top_df_large_strict = lr_network_top_df_large_strict %>% rename(ligand = from, receptor = to)
}

# DE analysis for each sender cell type -- of course only possible when having sender cell types
if (length(sender) > 1 | (length(sender) == 1 & sender != "undefined")){
DE_table_all = Idents(seurat_obj) %>% levels() %>% intersect(sender_celltypes) %>% lapply(get_lfc_celltype, seurat_obj = seurat_obj, condition_colname = condition_colname, condition_oi = condition_oi, condition_reference = condition_reference, expression_pct = expression_pct, celltype_col = NULL) %>% reduce(full_join, by = "gene") # use this if cell type labels are the identities of your Seurat object -- if not: indicate the celltype_col properly
DE_table_all[is.na(DE_table_all)] = 0

# Combine ligand activities with DE information
ligand_activities_de = ligand_activities %>% select(test_ligand, pearson) %>% rename(ligand = test_ligand) %>% left_join(DE_table_all %>% rename(ligand = gene), by = "ligand")
ligand_activities_de[is.na(ligand_activities_de)] = 0

# make LFC heatmap
lfc_matrix = ligand_activities_de %>% select(-ligand, -pearson) %>% as.matrix() %>% magrittr::set_rownames(ligand_activities_de$ligand)
rownames(lfc_matrix) = rownames(lfc_matrix) %>% make.names()

order_ligands = order_ligands[order_ligands %in% rownames(lfc_matrix)]
vis_ligand_lfc = lfc_matrix[order_ligands,]
vis_ligand_lfc = vis_ligand_lfc %>% as.matrix(ncol = length(Idents(seurat_obj) %>% levels() %>% intersect(sender_celltypes)))
colnames(vis_ligand_lfc) = vis_ligand_lfc %>% colnames() %>% make.names()

p_ligand_lfc = vis_ligand_lfc %>% make_threecolor_heatmap_ggplot("Prioritized ligands","LFC in Sender", low_color = "midnightblue",mid_color = "white", mid = median(vis_ligand_lfc), high_color = "red",legend_position = "top", x_axis_position = "top", legend_title = "LFC") + theme(axis.text.y = element_text(face = "italic"))
p_ligand_lfc

# change colors a bit to make them more stand out
p_ligand_lfc = p_ligand_lfc + scale_fill_gradientn(colors = c("midnightblue","blue", "grey95", "grey99","firebrick1","red"),values = c(0,0.1,0.2,0.25, 0.40, 0.7,1), limits = c(vis_ligand_lfc %>% min() - 0.1, vis_ligand_lfc %>% max() + 0.1))
p_ligand_lfc


# ligand expression Seurat dotplot
real_makenames_conversion = lr_network$from %>% unique() %>% magrittr::set_names(lr_network$from %>% unique() %>% make.names())
order_ligands_adapted = real_makenames_conversion[order_ligands]
names(order_ligands_adapted) = NULL
rotated_dotplot = DotPlot(seurat_obj %>% subset(idents = sender_celltypes), features = order_ligands_adapted, cols = "RdYlBu") + coord_flip() + theme(legend.text = element_text(size = 10), legend.title = element_text(size = 12)) # flip of coordinates necessary because we want to show ligands in the rows when combining all plots

# combined plot
figures_without_legend = cowplot::plot_grid(
p_ligand_pearson + theme(legend.position = "none", axis.ticks = element_blank()) + theme(axis.title.x = element_text()),
rotated_dotplot + theme(legend.position = "none", axis.ticks = element_blank(), axis.title.x = element_text(size = 12), axis.text.y = element_text(face = "italic", size = 9), axis.text.x = element_text(size = 9, angle = 90,hjust = 0)) + ylab("Expression in Sender") + xlab("") + scale_y_discrete(position = "right"),
p_ligand_lfc + theme(legend.position = "none", axis.ticks = element_blank()) + theme(axis.title.x = element_text()) + ylab(""),
p_ligand_target_network + theme(legend.position = "none", axis.ticks = element_blank()) + ylab(""),
align = "hv",
nrow = 1,
rel_widths = c(ncol(vis_ligand_pearson)+6, ncol(vis_ligand_lfc) + 7, ncol(vis_ligand_lfc) + 8, ncol(vis_ligand_target)))

legends = cowplot::plot_grid(
ggpubr::as_ggplot(ggpubr::get_legend(p_ligand_pearson)),
ggpubr::as_ggplot(ggpubr::get_legend(rotated_dotplot)),
ggpubr::as_ggplot(ggpubr::get_legend(p_ligand_lfc)),
ggpubr::as_ggplot(ggpubr::get_legend(p_ligand_target_network)),
nrow = 1,
align = "h", rel_widths = c(1.5, 1, 1, 1))

combined_plot = cowplot::plot_grid(figures_without_legend, legends, rel_heights = c(10,5), nrow = 2, align = "hv")
combined_plot

} else {
rotated_dotplot = NULL
p_ligand_lfc = NULL
}

return(list(
ligand_activities = ligand_activities,
top_ligands = best_upstream_ligands,
Expand All @@ -1039,6 +1098,8 @@ nichenet_seuratobj_aggregate = function(receiver, seurat_obj, condition_colname,
ligand_target_matrix = vis_ligand_target,
ligand_target_heatmap = p_ligand_target_network,
ligand_target_df = active_ligand_target_links_df,
ligand_expression_dotplot = rotated_dotplot,
ligand_differential_expression_heatmap = p_ligand_lfc,
ligand_activity_target_heatmap = combined_plot,
ligand_receptor_matrix = vis_ligand_receptor_network,
ligand_receptor_heatmap = p_ligand_receptor_network,
Expand Down Expand Up @@ -1356,8 +1417,9 @@ nichenet_seuratobj_cluster_de = function(seurat_obj, receiver_affected, receiver
## sender
if (length(sender) == 1){
if (sender == "all"){
list_expressed_genes_sender = Idents(seurat_obj) %>% unique() %>% lapply(get_expressed_genes, seurat_obj, expression_pct, assay_oi)
names(list_expressed_genes_sender) = Idents(seurat_obj) %>% unique()
sender_celltypes = Idents(seurat_obj) %>% levels()
list_expressed_genes_sender = sender_celltypes %>% lapply(get_expressed_genes, seurat_obj, expression_pct, assay_oi)
names(list_expressed_genes_sender) = sender_celltypes
expressed_genes_sender = list_expressed_genes_sender %>% unlist() %>% unique()

} else if (sender == "undefined") {
Expand Down Expand Up @@ -1572,6 +1634,18 @@ nichenet_seuratobj_cluster_de = function(seurat_obj, receiver_affected, receiver

}

# ligand expression Seurat dotplot
if (length(sender) > 1 | (length(sender) == 1 & sender != "undefined")){
real_makenames_conversion = lr_network$from %>% unique() %>% magrittr::set_names(lr_network$from %>% unique() %>% make.names())
order_ligands_adapted = real_makenames_conversion[order_ligands]
names(order_ligands_adapted) = NULL
rotated_dotplot = DotPlot(seurat_obj %>% subset(idents = sender_celltypes), features = order_ligands_adapted, cols = "RdYlBu") + coord_flip() + theme(legend.text = element_text(size = 10), legend.title = element_text(size = 12)) # flip of coordinates necessary because we want to show ligands in the rows when combining all plots

} else {
rotated_dotplot = NULL
}


return(list(
ligand_activities = ligand_activities,
top_ligands = best_upstream_ligands,
Expand All @@ -1580,6 +1654,7 @@ nichenet_seuratobj_cluster_de = function(seurat_obj, receiver_affected, receiver
ligand_target_matrix = vis_ligand_target,
ligand_target_heatmap = p_ligand_target_network,
ligand_target_df = active_ligand_target_links_df,
ligand_expression_dotplot = rotated_dotplot,
ligand_activity_target_heatmap = combined_plot,
ligand_receptor_matrix = vis_ligand_receptor_network,
ligand_receptor_heatmap = p_ligand_receptor_network,
Expand Down Expand Up @@ -1748,8 +1823,9 @@ nichenet_seuratobj_aggregate_cluster_de = function(seurat_obj, receiver_affected
## sender
if (length(sender) == 1){
if (sender == "all"){
list_expressed_genes_sender = Idents(seurat_obj) %>% unique() %>% lapply(get_expressed_genes, seurat_obj, expression_pct, assay_oi)
names(list_expressed_genes_sender) = Idents(seurat_obj) %>% unique()
sender_celltypes = Idents(seurat_obj) %>% levels()
list_expressed_genes_sender = sender_celltypes %>% lapply(get_expressed_genes, seurat_obj, expression_pct, assay_oi)
names(list_expressed_genes_sender) = sender_celltypes
expressed_genes_sender = list_expressed_genes_sender %>% unlist() %>% unique()

} else if (sender == "undefined") {
Expand Down Expand Up @@ -1973,6 +2049,16 @@ nichenet_seuratobj_aggregate_cluster_de = function(seurat_obj, receiver_affected

}

# ligand expression Seurat dotplot
if (length(sender) > 1 | (length(sender) == 1 & sender != "undefined")){
real_makenames_conversion = lr_network$from %>% unique() %>% magrittr::set_names(lr_network$from %>% unique() %>% make.names())
order_ligands_adapted = real_makenames_conversion[order_ligands]
names(order_ligands_adapted) = NULL
rotated_dotplot = DotPlot(seurat_obj %>% subset(idents = sender_celltypes), features = order_ligands_adapted, cols = "RdYlBu") + coord_flip() + theme(legend.text = element_text(size = 10), legend.title = element_text(size = 12)) # flip of coordinates necessary because we want to show ligands in the rows when combining all plots

} else {
rotated_dotplot = NULL
}
return(list(
ligand_activities = ligand_activities,
top_ligands = best_upstream_ligands,
Expand All @@ -1981,6 +2067,7 @@ nichenet_seuratobj_aggregate_cluster_de = function(seurat_obj, receiver_affected
ligand_target_matrix = vis_ligand_target,
ligand_target_heatmap = p_ligand_target_network,
ligand_target_df = active_ligand_target_links_df,
ligand_expression_dotplot = rotated_dotplot,
ligand_activity_target_heatmap = combined_plot,
ligand_receptor_matrix = vis_ligand_receptor_network,
ligand_receptor_heatmap = p_ligand_receptor_network,
Expand Down
25 changes: 11 additions & 14 deletions vignettes/seurat_wrapper.Rmd
View file
Original file line number Diff line number Diff line change
Expand Up @@ -102,6 +102,8 @@ As described in the main vignette, the pipeline of a basic NicheNet analysis con

All these steps are contained in one of three following similar single functions: `nichenet_seuratobj_aggregate`, `nichenet_seuratobj_cluster_de` and `nichenet_seuratobj_aggregate_cluster_de`.

In addition to these steps, the function `nichenet_seuratobj_aggregate` that is used for the analysis when having two conditions will also calculate differential expression of the ligands in the sender cell type. Note that this ligand differential expression is not used for prioritization and ranking of the ligands!

## NicheNet analysis on Seurat object: explain differential expression between two conditions

In this case study, the receiver cell population is the 'CD8 T' cell population, whereas the sender cell populations are 'CD4 T', 'Treg', 'Mono', 'NK', 'B' and 'DC'. The above described functions will consider a gene to be expressed when it is expressed in at least a predefined fraction of cells in one cluster (default: 10%).
Expand Down Expand Up @@ -146,21 +148,19 @@ nichenet_output$top_ligands

These ligands are expressed by one or more of the input sender cells. To see which cell population expresses which of these top-ranked ligands, you can run the following:

```{r, fig.width=12}
DotPlot(seuratObj, features = nichenet_output$top_ligands %>% rev(), cols = "RdYlBu") + RotatedAxis()
```{r, fig.width=8}
nichenet_output$ligand_expression_dotplot
```

As you can see, most op the top-ranked ligands seem to be mainly expressed by dendritic cells and monocytes.

It could also be interesting to see whether some of these ligands are differentially expressed after LCMV infection.

```{r, fig.width=12}
DotPlot(seuratObj, features = nichenet_output$top_ligands %>% rev(), split.by = "aggregate") + RotatedAxis()
```{r, fig.width=4, fig.height=8}
nichenet_output$ligand_differential_expression_heatmap
```

```{r}
VlnPlot(seuratObj, features = c("Il15", "Cxcl10","Cxcl16"), split.by = "aggregate", pt.size = 0, combine = FALSE)
```
As you can see, most op the top-ranked ligands seem also to be upregulated themselves in monocytes after viral infection. This is not a prerequisite to be top-ranked (cf: ranking only determined based on enrichment of target genes among DE genes in the receiver, CD8T cells), but is nice additional "evidence" that these ligands might indeed be important.


#### Inferred active ligand-target links
Expand Down Expand Up @@ -203,9 +203,9 @@ DotPlot(seuratObj %>% subset(idents = "CD8 T"), features = nichenet_output$top_t
VlnPlot(seuratObj %>% subset(idents = "CD8 T"), features = c("Zbp1","Ifit3","Irf7"), split.by = "aggregate", pt.size = 0, combine = FALSE)
```

To visualize both ligand activities and target genes of ligands, run the following command
To visualize ligand activities, expression, differentil expression and target genes of ligands, run the following command

```{r, fig.width = 12}
```{r, fig.width = 16}
nichenet_output$ligand_activity_target_heatmap
```

Expand Down Expand Up @@ -257,21 +257,18 @@ nichenet_output$background_expressed_genes %>% length()

Instead of focusing on multiple sender cell types, it is possible that you are only interested in doing the analyis for one sender cell type, such as dendritic cells in this case.

```{r}
```{r, fig.width=14}
nichenet_output = nichenet_seuratobj_aggregate(seurat_obj = seuratObj, receiver = "CD8 T", condition_colname = "aggregate", condition_oi = "LCMV", condition_reference = "SS", sender = "DC", ligand_target_matrix = ligand_target_matrix, lr_network = lr_network, weighted_networks = weighted_networks, organism = "mouse")
nichenet_output$ligand_activity_target_heatmap
DotPlot(seuratObj, features = nichenet_output$top_ligands %>% rev(), cols = "RdYlBu") + RotatedAxis()
```

Instead of focusing on one or multiple predefined sender cell types, it is also possible that you want to consider all cell types present as possible sender cell. This also includes the receiver cell type, making that you can look at autocrine signaling as well.

```{r, fig.width=8}
```{r, fig.width=14}
nichenet_output = nichenet_seuratobj_aggregate(seurat_obj = seuratObj, receiver = "CD8 T", condition_colname = "aggregate", condition_oi = "LCMV", condition_reference = "SS", sender = "all", ligand_target_matrix = ligand_target_matrix, lr_network = lr_network, weighted_networks = weighted_networks, organism = "mouse")
nichenet_output$ligand_activity_target_heatmap
DotPlot(seuratObj, features = nichenet_output$top_ligands %>% rev(), cols = "RdYlBu") + RotatedAxis()
```

In some cases, it could be possible that you don't have data of potential sender cells. If you still want to predict possible upstream ligands that could have been responsible for the observed differential expression in your cell type, you can do this by following command. This will consider all possible ligands in the NicheNet databases for which a receptor is expressed by the receiver cell of interest.
Expand Down
Loading

0 comments on commit 0bf85d6

Please sign in to comment.