hello，大家好，今天周一，我们不分享那些算法思路什么的，分享一点简单的可视化，专门针对cellphoneDB的结果进行可视化，不知道为什么，一直想起吕子乔的语录，谁年轻时没喜欢过几个人渣，😄😄
关于cellphoneDB，大家绝对低估了其价值，之前分享的有关cellphoneDB的重要内容分享在这里，供大家参考。
10X单细胞（10X空间转录组）细胞通讯的一个小小的细节
考虑空间位置的通讯分析手段---CellphoneDB（V3.0）
10X单细胞（10X空间转录组）空间相关性分析和cellphoneDB与NicheNet联合进行细胞通讯分析
关于cellphoneDB，目前升级到了3.0，已经用于分析生态位的通讯分析，当然官网也有很多的介绍，包括考虑复合物等等，也可以用自己定义的配受体数据集，算法上也算是非常经典，对结果文件也有解释，放在下面，供大家查阅。
User-defined receptor-ligand datasets
Our system allows users to create their own lists of curated proteins and complexes. In order to do so, the format of the users’ lists must be compatible with the input files. Users can submit their lists using the Python package version of CellPhoneDB, and then send them via email, the cellphonedb.org form, or a pull request to the CellPhoneDB data repository (https://github.com/Teichlab/cellphonedb-data).
Database structure
Information is stored in an SQLite relational database (https://www.sqlite.org). SQLAlchemy (www.sqlalchemy.org) and Python 3 were used to build the database structure and the query logic. The application is designed to allow analysis on potentially large count matrices to be performed in parallel. This requires an efficient database design, including optimisation for query times, indices and related strategies. All application code is open source and uploaded both to github and the web server.
Statistical inference of receptor-ligand specificity（显著性）
To assess cellular crosstalk between different cell types, we use our repository in a statistical framework for inferring cell–cell communication networks from single-cell transcriptome data. We predict enriched receptor–ligand interactions between two cell types based on expression of a receptor by one cell type and a ligand by another cell type, using scRNA-seq data. To identify the most relevant interactions between cell types, we look for the cell-type specific interactions between ligands and receptors. Only receptors and ligands expressed in more than a user-specified threshold percentage of the cells in the specific cluster are considered significant (default is 0.1).
We then perform pairwise comparisons between all cell types. First, we randomly permute the cluster labels of all cells (1,000 times as a default) and determine the mean of the average receptor expression level in a cluster and the average ligand expression level in the interacting cluster. For each receptor–ligand pair in each pairwise comparison between two cell types, this generates a null distribution. By calculating the proportion of the means which are as or higher than the actual mean, we obtain a p-value for the likelihood of cell-type specificity of a given receptor–ligand complex. We then prioritize interactions that are highly enriched between cell types based on the number of significant pairs, so that the user can manually select biologically relevant ones. For the multi-subunit heteromeric complexes, we require that all subunits of the complex are expressed (using a user-specified threshold), and therefore we use the member of the complex with the minimum average expression to perform the random shuffling.
Cell subsampling for accelerating analyses（下采样）
Technological developments and protocol improvements have enabled an exponential growth of the number of cells obtained from scRNA-seq experiments. Large-scale datasets can profile hundreds of thousands cells, which presents a challenge for the existing analysis methods in terms of both memory usage and runtime. In order to improve the speed and efficiency of our protocol and facilitate its broad accessibility, we integrated subsampling as described in Hie et al.28. This "geometric sketching" approach aims to maintain the transcriptomic heterogeneity within a dataset with a smaller subset of cells. The subsampling step is optional, enabling users to perform the analysis either on all cells, or with other subsampling methods of their choice.
结果文件的解释
Table. Description of the output files means.csv, pvalues.csv, significant_means.csv and relevant_interactions.txt
| Identifier | Definition | Output file | Example |
|---|---|---|---|
| id_cp_interaction | Unique CellPhoneDB identifier for each interaction stored in the database. | means.csv; pvalues.csv; significant_means.csv | CPI-SS096F3E0F2 |
| interacting_pair | Name of the interacting pairs separated by "|". | means.csv; pvalues.csv; significant_means.csv | JAG2|NOTCH4 |
| partner A or B | Identifier for the first interacting partner (A) or the second (B). It could be: UniProt (prefix simple:) or complex (prefix complex:) | means.csv; pvalues.csv; significant_means.csv | simple:Q9Y219 |
| gene A or B | Gene identifier for the first interacting partner (A) or the second (B). The identifier will depend on the input user list. | means.csv; pvalues.csv; significant_means.csv | ENSG00000184916 |
| secreted | True if one of the partners is secreted. | means.csv; pvalues.csv; significant_means.csv | FALSE |
| Receptor A or B | True if the first interacting partner (A) or the second (B) is annotated as a receptor in our database. | means.csv; pvalues.csv; significant_means.csv | FALSE |
| annotation_strategy | Curated if the interaction was annotated by the CellPhoneDB developers. Otherwise, the name of the database where the interaction has been downloaded from. | means.csv; pvalues.csv; significant_means.csv | curated |
| is_integrin | True if one of the partners is integrin. | means.csv; pvalues.csv; significant_means.csv | FALSE |
| rank | Total number of significant p-values for each interaction divided by the number of cell type-cell type comparisons. | significant_means.csv | 0.25 |
| means | Mean values for all the interacting partners: mean value refers to the total mean of the individual partner average expression values in the corresponding interacting pairs of cell types. If one of the mean values is 0, then the total mean is set to 0. | means.csv | 0.53 |
| p.values | p-values for the all the interacting partners: p.value refers to the enrichment of the interacting ligand-receptor pair in each of the interacting pairs of cell types. | pvalues.csv | 0.01 |
| significant_mean | Significant mean calculation for all the interacting partners. If p.value < 0.05, the value will be the mean. Alternatively, the value is set to 0. | significant_means.csv | 0.53 |
| relevant_interactions | Indicates if the interaction is relevant (1) or not (0). If a gene in the interaction is a DEG, and all the participants are expressed, the interaction will be classified as relevant. Alternatively, the value is set to 0. | relevant_interactions.txt | 1 or 0 |
Table. Description of the output file deconvoluted.csv
| Identifier | Definition | Output file | Example |
|---|---|---|---|
| gene_name | Gene identifier for one of the subunits that are participating in the interaction defined in "means.csv" file. The identifier will depend on the input of the user list. | deconvoluted.csv | JAG2 |
| uniprot | UniProt identifier for one of the subunits that are participating in the interaction defined in "means.csv" file. | deconvoluted.csv | Q9Y219 |
| is_complex | True if the subunit is part of a complex. Single if it is not, complex if it is. | deconvoluted.csv | FALSE |
| protein_name | Protein name for one of the subunits that are participating in the interaction defined in "means.csv" file. | deconvoluted.csv | JAG2_HUMAN |
| complex_name | Complex name if the subunit is part of a complex. Empty if not. | deconvoluted.csv | a10b1 complex |
| id_cp_interaction | Unique CellPhoneDB identifier for each of the interactions stored in the database. | deconvoluted.csv | CPI-SS0DB3F5A37 |
| mean | Mean expression of the corresponding gene in each cluster. | 0.9 |
最后软件默认的可视化
可视化进阶，这里介绍一个专门为cellphoneDB分析结果配套的可视化方法，软件在ktplots,这个软件的作者同时开发了BCR的集成分析软件dandelion,大家可以参考我的文章10X单细胞（10X空间转录组）BCR（TCR）数据分析之(11)dandelion
安装
if (!requireNamespace("devtools", quietly = TRUE))
install.packages("devtools")
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
devtools::install_github('zktuong/ktplots', dependencies = TRUE)
library(Seurat)
data(kidneyimmune)
library(ktplots)
plot_cpdb
This function seems like it's the most popular so I moved it up! Please see below for alternative visualisation options.
Generates a dot plot after CellPhoneDB analysis via specifying the query celltypes and genes. The difference compared to the original cellphonedb plot is that this is totally customizable!
The plotting is largely determined by the format of the meta file provided to CellPhoneDB analysis.
For the split.by option to work, the annotation in the meta file must be defined in the following format:
{split.by}_{idents}
so to set up an example vector, it would be something like:
annotation <- paste0(kidneyimmune$Experiment, '_', kidneyimmune$celltype)
To run, you will need to load in the means.txt and pvals.txt from the analysis. If you are using results from cellphonedb version 3, the pvalues.txt is relevant_interactions.txt and also add version3 = TRUE into all the functions below.
# pvals <- read.delim("pvalues.txt", check.names = FALSE)
# means <- read.delim("means.txt", check.names = FALSE)
# I've provided an example dataset
data(cpdb_output)
plot_cpdb(cell_type1 = 'B cell', cell_type2 = 'CD4T cell', scdata = kidneyimmune,
idents = 'celltype', # column name where the cell ids are located in the metadata
split.by = 'Experiment', # column name where the grouping column is. Optional.
means = means, pvals = pvals,
genes = c("XCR1", "CXCL10", "CCL5")) +
small_axis(fontsize = 3) + small_grid() + small_guide() + small_legend(fontsize = 2) # some helper functions included in ktplots to help with the plotting
plot_cpdb(cell_type1 = 'B cell', cell_type2 = 'CD4T cell', scdata = kidneyimmune,
idents = 'celltype', means = means, pvals = pvals, split.by = 'Experiment',
gene.family = 'chemokines') + small_guide() + small_axis() + small_legend(keysize=.5)
plot_cpdb(cell_type1 = 'B cell', cell_type2 = 'CD4T cell', scdata = kidneyimmune,
idents = 'celltype', means = means, pvals = pvals, split.by = 'Experiment',
gene.family = 'chemokines', col_option = "maroon", highlight = "blue") + small_guide() + small_axis() + small_legend(keysize=.5)
plot_cpdb(cell_type1 = 'B cell', cell_type2 = 'CD4T cell', scdata = kidneyimmune,
idents = 'celltype', means = means, pvals = pvals, split.by = 'Experiment',
gene.family = 'chemokines', col_option = viridis::cividis(50)) + small_guide() + small_axis() + small_legend(keysize=.5)
plot_cpdb(cell_type1 = 'B cell', cell_type2 = 'CD4T cell', scdata = kidneyimmune,
idents = 'celltype', means = means, pvals = pvals, split.by = 'Experiment',
gene.family = 'chemokines', noir = TRUE) + small_guide() + small_axis() + small_legend(keysize=.5)
plot_cpdb(cell_type1 = 'B cell', cell_type2 = 'CD4T cell', scdata = kidneyimmune,
idents = 'celltype', means = means, pvals = pvals, split.by = 'Experiment',
gene.family = 'chemokines', default_style = FALSE) + small_guide() + small_axis() + small_legend(keysize=.5)
plot_cpdb2
library(ktplots)
data(kidneyimmune)
data(cpdb_output2)
sce <- Seurat::as.SingleCellExperiment(kidneyimmune)
p <- plot_cpdb2(cell_type1 = 'B cell', cell_type2 = 'CD4T cell',
scdata = sce,
idents = 'celltype', # column name where the cell ids are located in the metadata
means = means2,
pvals = pvals2,
deconvoluted = decon2, # new options from here on specific to plot_cpdb2
desiredInteractions = list(
c('CD4T cell', 'B cell'),
c('B cell', 'CD4T cell')),
interaction_grouping = interaction_annotation,
edge_group_colors = c(
"Activating" = "#e15759",
"Chemotaxis" = "#59a14f",
"Inhibitory" = "#4e79a7",
"Intracellular trafficking" = "#9c755f",
"DC_development" = "#B07aa1",
"Unknown" = "#e7e7e7"
),
node_group_colors = c(
"CD4T cell" = "red",
"B cell" = "blue"),
keep_significant_only = TRUE,
standard_scale = TRUE,
remove_self = TRUE
)
p
# code example but not using the example datasets
library(SingleCellExperiment)
library(reticulate)
library(ktplots)
ad=import('anndata')
adata = ad$read_h5ad('rna.h5ad')
counts <- Matrix::t(adata$X)
row.names(counts) <- row.names(adata$var)
colnames(counts) <- row.names(adata$obs)
sce <- SingleCellExperiment(list(counts = counts), colData = adata$obs, rowData = adata$var)
means <- read.delim('out/means.txt', check.names = FALSE)
pvalues <- read.delim('out/pvalues.txt', check.names = FALSE)
deconvoluted <- read.delim('out/deconvoluted.txt', check.names = FALSE)
interaction_grouping <- read.delim('interactions_groups.txt')
# > head(interaction_grouping)
# interaction role
# 1 ALOX5_ALOX5AP Activating
# 2 ANXA1_FPR1 Inhibitory
# 3 BTLA_TNFRSF14 Inhibitory
# 4 CCL5_CCR5 Chemotaxis
# 5 CD2_CD58 Activating
# 6 CD28_CD86 Activating
test <- plot_cpdb2(cell_type1 = "CD4_Tem|CD4_Tcm|CD4_Treg", # same usage style as plot_cpdb
cell_type2 = "cDC",
idents = 'fine_clustering',
split.by = 'treatment_group_1',
scdata = sce,
means = means,
pvals = pvalues,
deconvoluted = deconvoluted, # new options from here on specific to plot_cpdb2
gene_symbol_mapping = 'index', # column name in rowData holding the actual gene symbols if the row names is ENSG Ids. Might be a bit buggy
desiredInteractions = list(c('CD4_Tcm', 'cDC1'), c('CD4_Tcm', 'cDC2'), c('CD4_Tem', 'cDC1'), c('CD4_Tem', 'cDC2 '), c('CD4_Treg', 'cDC1'), c('CD4_Treg', 'cDC2')),
interaction_grouping = interaction_grouping,
edge_group_colors = c("Activating" = "#e15759", "Chemotaxis" = "#59a14f", "Inhibitory" = "#4e79a7", " Intracellular trafficking" = "#9c755f", "DC_development" = "#B07aa1"),
node_group_colors = c("CD4_Tcm" = "#86bc86", "CD4_Tem" = "#79706e", "CD4_Treg" = "#ff7f0e", "cDC1" = "#bcbd22" ,"cDC2" = "#17becf"),
keep_significant_only = TRUE,
standard_scale = TRUE,
remove_self = TRUE)
compare_cpdb/plot_compare_cpdb
data(covid_sample_metadata)
data(covid_cpdb_meta)
file <- system.file("extdata", "covid_cpdb.tar.gz", package = "ktplots")
# copy and unpack wherever you want this to end up
file.copy(file, ".")
system("tar -xzf covid_cpdb.tar.gz")
It requires 1) an input table like so:
covid_cpdb_meta
sample_id cellphonedb_folder sce_file
1 MH9143325 covid_cpdb/MH9143325/out covid_cpdb/MH9143325/sce.rds
2 MH9143320 covid_cpdb/MH9143320/out covid_cpdb/MH9143320/sce.rds
3 MH9143274 covid_cpdb/MH9143274/out covid_cpdb/MH9143274/sce.rds
4 MH8919226 covid_cpdb/MH8919226/out covid_cpdb/MH8919226/sce.rds
5 MH8919227 covid_cpdb/MH8919227/out covid_cpdb/MH8919227/sce.rds
6 newcastle49 covid_cpdb/newcastle49/out covid_cpdb/newcastle49/sce.rds
7 MH9179822 covid_cpdb/MH9179822/out covid_cpdb/MH9179822/sce.rds
8 MH8919178 covid_cpdb/MH8919178/out covid_cpdb/MH8919178/sce.rds
9 MH8919177 covid_cpdb/MH8919177/out covid_cpdb/MH8919177/sce.rds
10 MH8919176 covid_cpdb/MH8919176/out covid_cpdb/MH8919176/sce.rds
11 MH8919179 covid_cpdb/MH8919179/out covid_cpdb/MH8919179/sce.rds
12 MH9179826 covid_cpdb/MH9179826/out covid_cpdb/MH9179826/sce.rds
where each row is a sample, the location of the out folder generated by cellphonedb and a single-cell object used to generate the cellphonedb result. If cellphonedb was ran with a .h5ad the sce_file would be the path to the .h5ad file.
and 2) a sample metadata data frame for the tests to run:
covid_sample_metadata
sample_id Status_on_day_collection_summary
MH9143325 MH9143325 Severe
MH9143320 MH9143320 Severe
MH9143274 MH9143274 Severe
MH8919226 MH8919226 Healthy
MH8919227 MH8919227 Healthy
newcastle49 newcastle49 Severe
MH9179822 MH9179822 Severe
MH8919178 MH8919178 Healthy
MH8919177 MH8919177 Healthy
MH8919176 MH8919176 Healthy
MH8919179 MH8919179 Healthy
MH9179826 MH9179826 Severe
重点So if I want to compare between Severe vs Healthy, I would specify the function as follows:
## set up the levels
covid_sample_metadata$Status_on_day_collection_summary <- factor(covid_sample_metadata$Status_on_day_collection_summary, levels = c('Healthy', 'Severe'))
out <- compare_cpdb(cpdb_meta = covid_cpdb_meta,
sample_metadata = covid_sample_metadata,
celltypes = c("B_cell", "CD14", "CD16", "CD4", "CD8", "DCs", "MAIT", "NK_16hi", "NK_56hi", "Plasmablast", "Platelets", "Treg", "gdT", "pDC"), # the actual celltypes you want to test
celltype_col = "initial_clustering", # the column that holds the cell type annotation in the sce object
groupby = "Status_on_day_collection_summary") # the column in the sample_metadata that holds the column where you want to do the comparison. In this example, it's Severe vs Healthy
This returns a list of dataframes (for each contrast found) with which you can use to plot the results.
plot_compare_cpdb is a simple function to achieve that but you can always just make a custom plotting function based on what you want.
plot_compare_cpdb(out) # red is significantly increased in Severe compared to Healthy.
If there are multiple contrasts and groups, you can facet the plot by specifying groups = c('group1', 'group2')
# let's mock up a new contrast like this
covid_sample_metadata$Status_on_day_collection_summary <- c(rep('Severe', 3), rep('Healthy', 2), rep('notSevere', 2), rep('Healthy', 4), 'notSevere')
out <- compare_cpdb(cpdb_meta = covid_cpdb_meta,
sample_metadata = covid_sample_metadata,
celltypes = c("B_cell", "CD14", "CD16", "CD4", "CD8", "DCs", "MAIT", "NK_16hi", "NK_56hi", "Plasmablast", "Platelets", "Treg", "gdT", "pDC"), # the actual celltypes you want to test
celltype_col = "initial_clustering", # the column that holds the cell type annotation in the sce object
groupby = "Status_on_day_collection_summary")
plot_compare_cpdb(out, alpha = .1, groups = names(out)) # there's no significant hit at 0.05 in this dummy example
The default method uses a pairwise wilcox.test. Alternatives are pairwise Welch's t.test or a linear mixed model with lmer.
To run the linear mixed effect model, it expects that the input data is paired (i.e an individual with multiple samples corresponding to multiple timepoints):
# just as a dummy example, lets say the samples are matched where there are two samples per individual
covid_sample_metadata$individual <- rep(c("A", "B", "C", "D", "E", "F"), 2)
# actually run it
out <- compare_cpdb(cpdb_meta = covid_cpdb_meta,
sample_metadata = covid_sample_metadata,
celltypes = c("B_cell", "CD14", "CD16", "CD4", "CD8", "DCs", "MAIT", "NK_16hi",
"NK_56hi", "Plasmablast", "Platelets", "Treg", "gdT", "pDC"),
celltype_col = "initial_clustering",
groupby = "Status_on_day_collection_summary",
formula = "~ Status_on_day_collection_summary + (1|individual)", # formula passed to lmer
method = "lmer")
plot_compare_cpdb(out, contrast = 'Status_on_day_collection_summarySevere') # use the colnames(out) to pick the right column.
Specifying cluster = TRUE will move the rows and columns to make it look a bit more ordered.
plot_compare_cpdb(out, contrast = 'Status_on_day_collection_summarySevere', cluster = TRUE)
其他绘图功能
# Note, this conflicts with tidyr devel version
geneDotPlot(scdata = kidneyimmune, # object
genes = c("CD68", "CD80", "CD86", "CD74", "CD2", "CD5"), # genes to plot
idents = "celltype", # column name in meta data that holds the cell-cluster ID/assignment
split.by = 'Project', # column name in the meta data that you want to split the plotting by. If not provided, it will just plot according to idents
standard_scale = TRUE) + # whether to scale expression values from 0 to 1. See ?geneDotPlot for other options
theme(strip.text.x = element_text(angle=0, hjust = 0, size =7)) + small_guide() + small_legend()
correlationSpot(空间相关性通讯分析)
library(ggplot2)
scRNAseq <- Seurat::SCTransform(scRNAseq, verbose = FALSE) %>% Seurat::RunPCA(., verbose = FALSE) %>% Seurat::RunUMAP(., dims = 1:30, verbose = FALSE)
anchors <- Seurat::FindTransferAnchors(reference = scRNAseq, query = spatial, normalization.method = "SCT")
predictions.assay <- Seurat::TransferData(anchorset = anchors, refdata = scRNAseq$label, dims = 1:30, prediction.assay = TRUE, weight.reduction = spatial[["pca"]])
spatial[["predictions"]] <- predictions.assay
Seurat::DefaultAssay(spatial) <- "predictions"
Seurat::DefaultAssay(spatial) <- 'SCT'
pa <- Seurat::SpatialFeaturePlot(spatial, features = c('Tnfsf13b', 'Cd79a'), pt.size.factor = 1.6, ncol = 2, crop = TRUE) + viridis::scale_fill_viridis()
Seurat::DefaultAssay(spatial) <- 'predictions'
pb <- Seurat::SpatialFeaturePlot(spatial, features = 'Group1-3', pt.size.factor = 1.6, ncol = 2, crop = TRUE) + viridis::scale_fill_viridis()
p1 <- correlationSpot(spatial, genes = c('Tnfsf13b', 'Cd79a'), celltypes = 'Group1-3', pt.size.factor = 1.6, ncol = 2, crop = TRUE) + scale_fill_gradientn( colors = rev(RColorBrewer::brewer.pal(12, 'Spectral')),limits = c(-1, 1))
p2 <- correlationSpot(spatial, genes = c('Tnfsf13b', 'Cd79a'), celltypes = 'Group1-3', pt.size.factor = 1.6, ncol = 2, crop = TRUE, average_by_cluster = TRUE) + scale_fill_gradientn(colors = rev(RColorBrewer::brewer.pal(12, 'Spectral')),limits = c(-1, 1)) + ggtitle('correlation averaged across clusters')
cowplot::plot_grid(pa, pb, p1, p2, ncol = 2)
好像改进的余地还挺多。
StackedVlnPlot
features <- c("CD79A", "MS4A1", "CD8A", "CD8B", "LYZ", "LGALS3", "S100A8", "GNLY", "NKG7", "KLRB1", "FCGR3A", "FCER1A", "CST3")
StackedVlnPlot(kidneyimmune, features = features) + theme(axis.text.x = element_text(angle = 90, hjust = 1, size = 8))
rainCloudPlot(data = kidneyimmune@meta.data, groupby = "celltype", parameter = "n_counts") + coord_flip()
生活很好，有你更好，百度文库出现了大量抄袭我的文章，对此我深表无奈，我写的文章，别人挂上去赚钱，抄袭可耻，挂到百度文库的人更可耻
