From 2e598492458fff64f668a90cd12d873e8d67c7a6 Mon Sep 17 00:00:00 2001 From: browaeysrobin Date: Tue, 19 Nov 2024 12:03:14 +0100 Subject: [PATCH] new scenario --- vignettes/DE_MultiNicheNet.Rmd | 568 +++++++++++++++++++++++++++++++++ 1 file changed, 568 insertions(+) create mode 100644 vignettes/DE_MultiNicheNet.Rmd diff --git a/vignettes/DE_MultiNicheNet.Rmd b/vignettes/DE_MultiNicheNet.Rmd new file mode 100644 index 0000000..d1f3209 --- /dev/null +++ b/vignettes/DE_MultiNicheNet.Rmd @@ -0,0 +1,568 @@ +--- +title: "DE analysis through MultiNicheNet" +author: "Robin Browaeys" +package: "`r BiocStyle::pkg_ver('multinichenetr')`" +output: + BiocStyle::html_document +output_dir: "/Users/robinb/Work/multinichenetr/vignettes" +vignette: > + %\VignetteIndexEntry{DE analysis through MultiNicheNet} + %\VignetteEngine{knitr::rmarkdown} + %\VignetteEncoding{UTF-8} +--- + + + + + +```{r setup, include = FALSE} +knitr::opts_chunk$set( + collapse = TRUE, + # comment = "#>", + warning = FALSE, + message = FALSE +) +library(BiocStyle) +``` + +In this vignette, you can learn how to perform a MultiNicheNet analysis to compare cell-cell communication between conditions of interest. A MultiNicheNet analysis can be performed if you have multi-sample, multi-condition/group single-cell data. We strongly recommend having at least 4 samples in each of the groups/conditions you want to compare (3 samples absolute minimum). With less samples, the benefits of performing a pseudobulk-based DE analysis are less clear. + +As input you need a SingleCellExperiment object containing at least the raw count matrix and metadata providing the following information for each cell: the **group**, **sample** and **cell type**. + +As example expression data of interacting cells, we will here use scRNAseq data of immune cells in MIS-C patients and healthy siblings from this paper of Hoste et al.: [TIM3+ TRBV11-2 T cells and IFNγ signature in patrolling monocytes and CD16+ NK cells delineate MIS-C](https://rupress.org/jem/article/219/2/e20211381/212918/TIM3-TRBV11-2-T-cells-and-IFN-signature-in) [![DOI](https://zenodo.org/badge/DOI/10.5281/zenodo.6362434.svg)](https://doi.org/10.5281/zenodo.6362434). +MIS-C (multisystem inflammatory syndrome in children) is a novel rare immunodysregulation syndrome that can arise after SARS-CoV-2 infection in children. We will use MultiNicheNet to explore immune cell crosstalk enriched in MIS-C compared to healthy siblings and adult COVID-19 patients. + +In this vignette, we will first prepare the MultiNicheNet core analysis, then run the several steps in the MultiNicheNet core analysis, and finally interpret the output. + +# Preparation of the MultiNicheNet core analysis + +```{r load-libs, message = FALSE, warning = FALSE} +library(SingleCellExperiment) +library(dplyr) +library(ggplot2) +library(nichenetr) +library(multinichenetr) +``` + +## Load NicheNet's ligand-receptor network and ligand-target matrix + +MultiNicheNet builds upon the NicheNet framework and uses the same prior knowledge networks (ligand-receptor network and ligand-target matrix, currently v2 version).+ + +The Nichenet v2 networks and matrices for both mouse and human can be downloaded from Zenodo [![DOI](https://zenodo.org/badge/DOI/10.5281/zenodo.7074291.svg)](https://doi.org/10.5281/zenodo.7074291). + +We will read these object in for human because our expression data is of human patients. +Gene names are here made syntactically valid via `make.names()` to avoid the loss of genes (eg H2-M3) in downstream visualizations. + +```{r} +organism = "human" +``` + +```{r, results='hide'} +options(timeout = 240) + +if(organism == "human"){ + + lr_network_all = + readRDS(url( + "https://zenodo.org/record/10229222/files/lr_network_human_allInfo_30112033.rds" + )) %>% + mutate( + ligand = convert_alias_to_symbols(ligand, organism = organism), + receptor = convert_alias_to_symbols(receptor, organism = organism)) + + lr_network_all = lr_network_all %>% + mutate(ligand = make.names(ligand), receptor = make.names(receptor)) + + lr_network = lr_network_all %>% + distinct(ligand, receptor) + + ligand_target_matrix = readRDS(url( + "https://zenodo.org/record/7074291/files/ligand_target_matrix_nsga2r_final.rds" + )) + + colnames(ligand_target_matrix) = colnames(ligand_target_matrix) %>% + convert_alias_to_symbols(organism = organism) %>% make.names() + rownames(ligand_target_matrix) = rownames(ligand_target_matrix) %>% + convert_alias_to_symbols(organism = organism) %>% make.names() + + lr_network = lr_network %>% filter(ligand %in% colnames(ligand_target_matrix)) + ligand_target_matrix = ligand_target_matrix[, lr_network$ligand %>% unique()] + +} else if(organism == "mouse"){ + + lr_network_all = readRDS(url( + "https://zenodo.org/record/10229222/files/lr_network_mouse_allInfo_30112033.rds" + )) %>% + mutate( + ligand = convert_alias_to_symbols(ligand, organism = organism), + receptor = convert_alias_to_symbols(receptor, organism = organism)) + + lr_network_all = lr_network_all %>% + mutate(ligand = make.names(ligand), receptor = make.names(receptor)) + lr_network = lr_network_all %>% + distinct(ligand, receptor) + + ligand_target_matrix = readRDS(url( + "https://zenodo.org/record/7074291/files/ligand_target_matrix_nsga2r_final_mouse.rds" + )) + + colnames(ligand_target_matrix) = colnames(ligand_target_matrix) %>% + convert_alias_to_symbols(organism = organism) %>% make.names() + rownames(ligand_target_matrix) = rownames(ligand_target_matrix) %>% + convert_alias_to_symbols(organism = organism) %>% make.names() + + lr_network = lr_network %>% filter(ligand %in% colnames(ligand_target_matrix)) + ligand_target_matrix = ligand_target_matrix[, lr_network$ligand %>% unique()] + +} +``` + +## Read in SingleCellExperiment Object + +In this vignette, we will load in a subset of the scRNAseq data of the MIS-C [![DOI](https://zenodo.org/badge/DOI/10.5281/zenodo.8010790.svg)](https://doi.org/10.5281/zenodo.8010790). For the sake of demonstration, this subset only contains 3 cell types. These celltypes are some of the cell types that were found to be most interesting related to MIS-C according to Hoste et al. + +If you start from a Seurat object, you can convert it easily to a SingleCellExperiment object via `sce = Seurat::as.SingleCellExperiment(seurat_obj, assay = "RNA")`. + +Because the NicheNet 2.0. networks are in the most recent version of the official gene symbols, we will make sure that the gene symbols used in the expression data are also updated (= converted from their "aliases" to official gene symbols). Afterwards, we will make them again syntactically valid. + +```{r, results='hide'} +sce = readRDS(url( + "https://zenodo.org/record/8010790/files/sce_subset_misc.rds" + )) +sce = alias_to_symbol_SCE(sce, "human") %>% makenames_SCE() # this can be removed +``` + +## Prepare the settings of the MultiNicheNet cell-cell communication analysis + +In this step, we will formalize our research question into MultiNicheNet input arguments. + +### Define in which metadata columns we can find the **group**, **sample** and **cell type** IDs + +In this case study, we want to study differences in cell-cell communication patterns between MIS-C patients (M), their healthy siblings (S) and adult patients with severe covid (A). The meta data columns that indicate this disease status(=group/condition of interest) is `MIS.C.AgeTier`. + +Cell type annotations are indicated in the `Annotation_v2.0` column, and the sample is indicated by the `ShortID` column. +If your cells are annotated in multiple hierarchical levels, we recommend using a relatively high level in the hierarchy. This for 2 reasons: 1) MultiNicheNet focuses on differential expression and not differential abundance, and 2) there should be sufficient cells per sample-celltype combination (see later). + +```{r} +sample_id = "ShortID" +group_id = "MIS.C.AgeTier" +celltype_id = "Annotation_v2.0" +``` + +__Important__: It is required that each sample-id is uniquely assigned to only one condition/group of interest. See the vignettes about paired and multifactorial analysis to see how to define your analysis input when you have multiple samples (and conditions) per patient. + +If you would have batch effects or covariates you can correct for, you can define this here as well. However, this is not applicable to this dataset. Therefore we will use the following NA settings: + +```{r} +covariates = NA +batches = NA +``` + +__Important__: for categorical covariates and batches, there should be at least one sample for every group-batch combination. If one of your groups/conditions lacks a certain level of your batch, you won't be able to correct for the batch effect because the model is then not able to distinguish batch from group/condition effects. + +__Important__: The column names of group, sample, cell type, batches and covariates should be syntactically valid (`make.names`) + +__Important__: All group, sample, cell type, batch and covariate names should be syntactically valid as well (`make.names`) (eg through `SummarizedExperiment::colData(sce)$ShortID = SummarizedExperiment::colData(sce)$ShortID %>% make.names()`) + +### Define the contrasts of interest. + +Here, we want to compare each patient group to the other groups, so the MIS-C (M) group vs healthy control siblings (S) and adult COVID19 patients (A) (= M vs S+A) and so on. We want to know which cell-cell communication patterns are specific for the M vs A+S group, the A vs M+S group and the S vs A+M group. + +To perform this comparison, we need to set the following contrasts: + +```{r} +contrasts_oi = c("'M-(S+A)/2','S-(M+A)/2','A-(S+M)/2'") +``` + +__Very Important__ Note the format to indicate the contrasts! This formatting should be adhered to very strictly, and white spaces are not allowed! Check `?get_DE_info` for explanation about how to define this well. The most important points are that: +*each contrast is surrounded by single quotation marks +*contrasts are separated by a comma without any white space +*all contrasts together are surrounded by double quotation marks. + +If you compare against two groups, you should divide by 2 (as demonstrated here), if you compare against three groups, you should divide by 3 and so on. + +For downstream visualizations and linking contrasts to their main condition, we also need to run the following: +This is necessary because we will also calculate cell-type+condition specificity of ligands and receptors. + +```{r} +contrast_tbl = tibble(contrast = + c("M-(S+A)/2","S-(M+A)/2", "A-(S+M)/2"), + group = c("M","S","A")) +``` + +If you want to compare only two groups (eg M vs S), you can use the following: +`contrasts_oi = c("'M-S','S-M'") ` +`contrast_tbl = tibble(contrast = c("M-S","S-M"), group = c("M","S"))` + +Other vignettes will demonstrate how to formalize different types of research questions. + +### Define the sender and receiver cell types of interest. + +If you want to focus the analysis on specific cell types (e.g. because you know which cell types reside in the same microenvironments based on spatial data), you can define this here. If you have sufficient computational resources and no specific idea of cell-type colocalzations, we recommend to consider all cell types as potential senders and receivers. Later on during analysis of the output it is still possible to zoom in on the cell types that interest you most, but your analysis is not biased to them. + +Here we will consider all cell types in the data: + +```{r} +senders_oi = SummarizedExperiment::colData(sce)[,celltype_id] %>% unique() +receivers_oi = SummarizedExperiment::colData(sce)[,celltype_id] %>% unique() +sce = sce[, SummarizedExperiment::colData(sce)[,celltype_id] %in% + c(senders_oi, receivers_oi) + ] +``` + +In case you would have samples in your data that do not belong to one of the groups/conditions of interest, we recommend removing them and only keeping conditions of interst: + +```{r} +conditions_keep = c("M", "S", "A") +sce = sce[, SummarizedExperiment::colData(sce)[,group_id] %in% + conditions_keep + ] +``` + +```{r} +metadata_combined = SummarizedExperiment::colData(sce) %>% tibble::as_tibble() + +if(!is.na(batches)){ + grouping_tbl = metadata_combined[,c(sample_id, group_id, batches)] %>% + tibble::as_tibble() %>% distinct() + colnames(grouping_tbl) = c("sample","group",batches) +} else { + grouping_tbl = metadata_combined[,c(sample_id, group_id)] %>% + tibble::as_tibble() %>% distinct() + colnames(grouping_tbl) = c("sample","group") +} +``` + +# Running the MultiNicheNet core analysis + +Now we will run the core of a MultiNicheNet analysis, here until the DE analysis step. This analysis consists of the following steps: + +* 1. Cell-type filtering: determine which cell types are sufficiently present +* 2. Gene filtering: determine which genes are sufficiently expressed in each present cell type +* 3. Pseudobulk expression calculation: determine and normalize per-sample pseudobulk expression levels for each expressed gene in each present cell type +* 4. Differential expression (DE) analysis: determine which genes are differentially expressed + +After these steps, the output can be further explored as we will demonstrate in the "Downstream analysis of the MultiNicheNet output" section. + +In this vignette, we will demonstrate these steps one-by-one, which offers the most flexibility to the user to assess intermediary results. Other vignettes will demonstrate the use of the `multi_nichenet_analysis` wrapper function. + +## Cell-type filtering: determine which cell types are sufficiently present + +In this step we will calculate and visualize cell type abundances. This will give an indication about which cell types will be retained in the analysis, and which cell types will be filtered out. + +Since MultiNicheNet will infer group differences at the sample level for each cell type (currently via Muscat - pseudobulking + EdgeR), we need to have sufficient cells per sample of a cell type, and this for all groups. In the following analysis we will set this minimum number of cells per cell type per sample at 10. Samples that have less than `min_cells` cells will be excluded from the analysis for that specific cell type. + +```{r} +min_cells = 10 +``` + +We recommend using `min_cells = 10`, except for datasets with several lowly abundant cell types of interest. For those datasets, we recommend using `min_cells = 5`. + +```{r} +abundance_info = get_abundance_info( + sce = sce, + sample_id = sample_id, group_id = group_id, celltype_id = celltype_id, + min_cells = min_cells, + senders_oi = senders_oi, receivers_oi = receivers_oi, + batches = batches + ) +``` + +First, we will check the cell type abundance diagnostic plots. + +### Interpretation of cell type abundance information + +The first plot visualizes the number of cells per celltype-sample combination, and indicates which combinations are removed during the DE analysis because there are less than `min_cells` in the celltype-sample combination. + +```{r} +abundance_info$abund_plot_sample +``` + +The red dotted line indicates the required minimum of cells as defined above in `min_cells`. We can see here that some sample-celltype combinations are left out. For the DE analysis in the next step, only cell types will be considered if there are at least two samples per group with a sufficient number of cells. But as we can see here: all cell types will be considered for the analysis and there are no condition-specific cell types. + +__Important__: Based on the cell type abundance diagnostics, we recommend users to change their analysis settings if required (eg changing cell type annotation level, batches, ...), before proceeding with the rest of the analysis. If too many celltype-sample combinations don't pass this threshold, we recommend to define your cell types in a more general way (use one level higher of the cell type ontology hierarchy) (eg TH17 CD4T cells --> CD4T cells) or use `min_cells = 5` if this would not be possible. + +### Cell type filtering based on cell type abundance information + +Running the following block of code can help you determine which cell types are condition-specific and which cell types are absent. + +```{r} +abundance_df_summarized = abundance_info$abundance_data %>% + mutate(keep = as.logical(keep)) %>% + group_by(group_id, celltype_id) %>% + summarise(samples_present = sum((keep))) + +celltypes_absent_one_condition = abundance_df_summarized %>% + filter(samples_present == 0) %>% pull(celltype_id) %>% unique() +# find truly condition-specific cell types by searching for cell types +# truely absent in at least one condition + +celltypes_present_one_condition = abundance_df_summarized %>% + filter(samples_present >= 2) %>% pull(celltype_id) %>% unique() +# require presence in at least 2 samples of one group so +# it is really present in at least one condition + +condition_specific_celltypes = intersect( + celltypes_absent_one_condition, + celltypes_present_one_condition) + +total_nr_conditions = SummarizedExperiment::colData(sce)[,group_id] %>% + unique() %>% length() + +absent_celltypes = abundance_df_summarized %>% + filter(samples_present < 2) %>% + group_by(celltype_id) %>% + count() %>% + filter(n == total_nr_conditions) %>% + pull(celltype_id) + +print("condition-specific celltypes:") +print(condition_specific_celltypes) + +print("absent celltypes:") +print(absent_celltypes) +``` +Absent cell types will be filtered out, condition-specific cell types can be filtered out if you as a user do not want to run the alternative workflow for condition-specific cell types in the optional step 8 of the core MultiNicheNet analysis. + +```{r} +analyse_condition_specific_celltypes = FALSE +``` + +```{r} +if(analyse_condition_specific_celltypes == TRUE){ + senders_oi = senders_oi %>% setdiff(absent_celltypes) + receivers_oi = receivers_oi %>% setdiff(absent_celltypes) +} else { + senders_oi = senders_oi %>% + setdiff(union(absent_celltypes, condition_specific_celltypes)) + receivers_oi = receivers_oi %>% + setdiff(union(absent_celltypes, condition_specific_celltypes)) +} + +sce = sce[, SummarizedExperiment::colData(sce)[,celltype_id] %in% + c(senders_oi, receivers_oi) + ] +``` + +## Gene filtering: determine which genes are sufficiently expressed in each present cell type + +Before running the DE analysis, we will determine which genes are not sufficiently expressed and should be filtered out. +We will perform gene filtering based on a similar procedure as used in `edgeR::filterByExpr`. However, we adapted this procedure to be more interpretable for single-cell datasets. + +For each cell type, we will consider genes expressed if they are expressed in at least a `min_sample_prop` fraction of samples in the condition with the lowest number of samples. By default, we set `min_sample_prop = 0.50`, which means that genes should be expressed in at least 2 samples if the group with lowest nr. of samples has 4 samples like this dataset. + +```{r} +min_sample_prop = 0.50 +``` + +But how do we define which genes are expressed in a sample? For this we will consider genes as expressed if they have non-zero expression values in a `fraction_cutoff` fraction of cells of that cell type in that sample. By default, we set `fraction_cutoff = 0.05`, which means that genes should show non-zero expression values in at least 5% of cells in a sample. + +```{r} +fraction_cutoff = 0.05 +``` + +We recommend using these default values unless there is specific interest in prioritizing (very) weakly expressed interactions. In that case, you could lower the value of `fraction_cutoff`. We explicitly recommend against using `fraction_cutoff > 0.10`. + +Now we will calculate the information required for gene filtering with the following command: + +```{r} +frq_list = get_frac_exprs( + sce = sce, + sample_id = sample_id, celltype_id = celltype_id, group_id = group_id, + batches = batches, + min_cells = min_cells, + fraction_cutoff = fraction_cutoff, min_sample_prop = min_sample_prop) +``` +Now only keep genes that are expressed by at least one cell type: + +```{r} +genes_oi = frq_list$expressed_df %>% + filter(expressed == TRUE) %>% pull(gene) %>% unique() +sce = sce[genes_oi, ] +``` + +## Pseudobulk expression calculation: determine and normalize per-sample pseudobulk expression levels for each expressed gene in each present cell type + +After filtering out absent cell types and genes, we will continue the analysis by calculating the different prioritization criteria that we will use to prioritize cell-cell communication patterns. + +First, we will determine and normalize per-sample pseudobulk expression levels for each expressed gene in each present cell type. The function `process_abundance_expression_info` will link this expression information for ligands of the sender cell types to the corresponding receptors of the receiver cell types. This will later on allow us to define the cell-type specicificy criteria for ligands and receptors. + +```{r} +abundance_expression_info = process_abundance_expression_info( + sce = sce, + sample_id = sample_id, group_id = group_id, celltype_id = celltype_id, + min_cells = min_cells, + senders_oi = senders_oi, receivers_oi = receivers_oi, + lr_network = lr_network, + batches = batches, + frq_list = frq_list, + abundance_info = abundance_info) +``` + +Normalized pseudobulk expression values per gene/celltype/sample can be inspected by: + +```{r} +abundance_expression_info$celltype_info$pb_df %>% head() +``` + +An average of these sample-level expression values per condition/group can be inspected by: + +```{r} +abundance_expression_info$celltype_info$pb_df_group %>% head() +``` + +Inspecting these values for ligand-receptor interactions can be done by: + +## Differential expression (DE) analysis: determine which genes are differentially expressed + +In this step, we will perform genome-wide differential expression analysis of receiver and sender cell types to define DE genes between the conditions of interest (as formalized by the `contrasts_oi`). Based on this analysis, we later can define the levels of differential expression of ligands in senders and receptors in receivers, and define the set of affected target genes in the receiver cell types (which will be used for the ligand activity analysis). + +We will apply pseudobulking followed by EdgeR to perform multi-condition multi-sample differential expression (DE) analysis (also called 'differential state' analysis by the developers of Muscat). + +```{r} +DE_info = get_DE_info( + sce = sce, + sample_id = sample_id, group_id = group_id, celltype_id = celltype_id, + batches = batches, covariates = covariates, + contrasts_oi = contrasts_oi, + min_cells = min_cells, + expressed_df = frq_list$expressed_df) +``` + +### Check DE results + +Check DE output information in table with logFC and p-values for each gene-celltype-contrast: + +```{r} +DE_info$celltype_de$de_output_tidy %>% head() +``` +Evaluate the distributions of p-values: + +```{r} +DE_info$hist_pvals +``` + +These distributions look fine (uniform distribution, except peak at p-value <= 0.05), so we will continue using these regular p-values. In case these p-value distributions look irregular, you can estimate empirical p-values as we will demonstrate in another vignette. + +```{r} +empirical_pval = FALSE +``` + +```{r} +if(empirical_pval == TRUE){ + DE_info_emp = get_empirical_pvals(DE_info$celltype_de$de_output_tidy) + celltype_de = DE_info_emp$de_output_tidy_emp %>% select(-p_val, -p_adj) %>% + rename(p_val = p_emp, p_adj = p_adj_emp) +} else { + celltype_de = DE_info$celltype_de$de_output_tidy +} +``` + +# Interpreting the MultiNicheNet analysis output + +## Necessary code to be able to use downstream plotting function without having run the entire MultiNicheNet pipeline + +```{r} +sample_prioritization_tbl = abundance_expression_info$abundance_data_receiver +sample_prioritization_tbl$n_cells_receiver[is.na(sample_prioritization_tbl$n_cells_receiver)] = 0 +sample_prioritization_tbl$keep_receiver[is.na(sample_prioritization_tbl$keep_receiver)] = 0 + +prioritization_tables = list(sample_prioritization_tbl = sample_prioritization_tbl) +``` +## Visualize top DE genes for a cell type of interest + +Finally, we provide some visualizations to just inspect the DE results that were generated during the MultiNicheNet analysis. + +```{r, fig.width=6, fig.height=16.5} +group_oi = "M" +celltype_oi = "M_Monocyte_CD16" + +DE_genes = DE_info$celltype_de$de_output_tidy %>% + inner_join(contrast_tbl) %>% + filter(group == group_oi) %>% + arrange(p_val) %>% + filter( + cluster_id == celltype_oi & + logFC > 2 & + p_val <= 0.05 & + contrast == contrast_tbl %>% filter(group == group_oi) %>% pull(contrast)) %>% + pull(gene) %>% unique() + +p_target = make_DEgene_dotplot_pseudobulk( + genes_oi = DE_genes, + celltype_info = abundance_expression_info$celltype_info, + prioritization_tables = prioritization_tables, + celltype_oi = celltype_oi, + grouping_tbl = grouping_tbl) + + +# create option for facetting cell types + +p_target$pseudobulk_plot + ggtitle("DE genes (pseudobulk expression)") +p_target$singlecell_plot + ggtitle("DE genes (single-cell expression)") +``` + +Among these DE genes, you may be most interested in ligands or receptors + +Ligands: +```{r, fig.width=7, fig.height=7} +group_oi = "M" +celltype_oi = "M_Monocyte_CD16" + +DE_genes = DE_info$celltype_de$de_output_tidy %>% + inner_join(contrast_tbl) %>% + filter(group == group_oi) %>% + arrange(p_val) %>% + filter( + cluster_id == celltype_oi & + logFC > 1 & + p_val <= 0.05 & + contrast == contrast_tbl %>% filter(group == group_oi) %>% pull(contrast)) %>% + pull(gene) %>% unique() +DE_genes = DE_genes %>% intersect(lr_network$ligand) + +p_target = make_DEgene_dotplot_pseudobulk( + genes_oi = DE_genes, + celltype_info = abundance_expression_info$celltype_info, + prioritization_tables = prioritization_tables, + celltype_oi = celltype_oi, + grouping_tbl = grouping_tbl) + +p_target$pseudobulk_plot + ggtitle("DE ligands (pseudobulk expression)") +p_target$singlecell_plot + ggtitle("DE ligands (single-cell expression)") +``` + +Receptors: +```{r, fig.width=7, fig.height=7} +group_oi = "M" +celltype_oi = "M_Monocyte_CD16" + +DE_genes = DE_info$celltype_de$de_output_tidy %>% + inner_join(contrast_tbl) %>% + filter(group == group_oi) %>% + arrange(p_val) %>% + filter( + cluster_id == celltype_oi & + logFC > 1 & + p_val <= 0.05 & + contrast == contrast_tbl %>% filter(group == group_oi) %>% pull(contrast)) %>% + pull(gene) %>% unique() +DE_genes = DE_genes %>% intersect(lr_network$receptor) + +p_target = make_DEgene_dotplot_pseudobulk( + genes_oi = DE_genes, + celltype_info = abundance_expression_info$celltype_info, + prioritization_tables = prioritization_tables, + celltype_oi = celltype_oi, + grouping_tbl = grouping_tbl) + +p_target$pseudobulk_plot + ggtitle("DE receptors (pseudobulk expression)") +p_target$singlecell_plot + ggtitle("DE receptors (single-cell expression)") +``` +