PlantPhoneDB: A manually curated pan-plant database of ligand-receptor pairs infers cell-cell communication

Abstract

Ligand-receptor pairs play important roles in cell-cell communication for multicellular organisms in response to environmental cues. Recently, the emergence of single-cell RNA-sequencing (scRNA-seq) provides unprecedented opportunities to investigate cellular communication based on ligand-receptor expression. However, so far, no reliable ligand-receptor interaction database is available for plant species. In this study, we developed PlantPhoneDB (https://jasonxu.shinyapps.io/PlantPhoneDB/), a pan-plant database comprising a large number of high-confidence ligand-receptor pairs manually curated from seven resources. Also, we developed a PlantPhoneDB R package, which not only provided optional four scoring approaches that calculate interaction scores of ligand-receptor pairs between cell types but also provided visualization functions to present analysis results. At the PlantPhoneDB web interface, the processed datasets and results can be searched, browsed, and downloaded. To uncover novel cell-cell communication events in plants, we applied the PlantPhoneDB R package on GSE121619 dataset to infer significant cell-cell interactions of heat-shocked root cells in Arabidopsis thaliana. As a result, the PlantPhoneDB predicted the actively communicating AT1G28290-AT2G14890 ligand-receptor pair in atrichoblast-cortex cell pair in Arabidopsis thaliana. Importantly, the downstream target genes of this ligand-receptor pair were significantly enriched in the ribosome pathway, which facilitated plants adapting to environmental changes. In conclusion, PlantPhoneDB provided researchers with integrated resources to infer cell-cell communication from scRNA-seq datasets.

Keywords: cell-cell communication; ligand-receptor interactions; plants; signalling pathway; single-cell transcriptomics.

Figure 1

Statistics of PlantPhoneDB and summary of scRNA‐seq datasets were analysed. (a) The number of ligand‐receptor pairs curated from plant.MAP, Interactome v2.0, IntAct, BioGRID, Text‐mining from literatures, STRING, and Orthologs resources in Arabidopsis thaliana. And 3514 unique ligand‐receptor pairs are used to infer cell–cell communication. (b) The number of ligands, receptors and ligand‐receptor pairs identified in 5 plant species, including Arabidopsis thaliana, Oryza sativa, Populus alba x Populus glandulosa, Solanum lycopersicum, and Zea mays. (c) PlantPhoneDB includes 29 scRNA‐seq datasets information, covering ~560 000 cells of 15 tissues from 5 plant species. FAIL, PASS, and pending datasets are indicated as blue, black, and purple bar, respectively. PASS datasets indicate scRNA‐seq datasets with ≥1000 high‐quality cells; pending datasets indicate PASS datasets without available the expression matrix or datasets are too large to be analysed on our laptop. The rest of scRNA‐seq datasets were considered to be FAIL datasets (<1000 high‐quality cells). The number of cells (recorded by original paper) for each dataset is shown inside the parenthesis.

Figure 2

Functions of PlantPhoneDB web interface and a visualization example of AT2G43610 gene expression across cell types. (a) Overview of PlantPhoneDB. Seven modules are displayed on the navigation bar. Number of ligand‐receptor pairs from five plant species are collected in PlantPhoneDB. The detailed information of scRNA‐seq datasets and resources are showed in the box. (b) Visualization example of GSE114615 dataset using cellxgene software. The detailed meta‐information for each dataset was displayed on the left, such as annotated cell identities and treatment conditions. On the right, the UMAP plot of GSE114615 dataset with cells coloured by trichoblast, atrichoblast, lateral root, meristem, endodermis, and cortex cell types. Each dot represents one cell. (c) A violin plot shows the distribution of AT2G43610 expression level across six different cell types.

Figure 3

Benchmarks the performance of 10 classifiers across 7 PbmcBench datasets. (a) Heatmap shows the median F1‐scores of 10 classifiers for 42 pairwise train‐test combination across different protocols (inter‐datasets model). Datasets on the top of the heatmap are used as training datasets, and testing datasets are indicated on the bottom of the heatmap. The inter‐datasets model indicated trained scRNA‐seq dataset from one sequencing protocol was used to predict cell type of scRNA‐seq dataset from another sequencing protocol. (b) Median F1‐scores of 10 classifiers within a dataset of different protocols (intra‐dataset model), including 10 × v2, 10 × v3, CEL_Seq, Drop_Seq, inDrop, Seq_Well and Smart_Seq2 protocol. The intra‐dataset model indicated trained scRNA‐seq dataset from one sequencing protocol was used to predict cell type of scRNA‐seq dataset from the same sequencing protocol. (c) Evaluates mean computation time and mean F1‐scores of each classifier. Barplot indicates mean running time of each classifier (left); line plot indicates mean F1‐scores (right).

Figure 4

Comparison of four scoring approaches (LRscore, WeightProduct, Average, and Product). (a) The number of ligand‐receptor pairs identified using four scoring approaches on 8 k 10 × PBMCs dataset. Rows represent cells expressing the receptors and columns represent cells expressing the ligands. Low and high number of ligand‐receptor pair are showed by purple and yellow, respectively. (b) The number of ligand‐receptor pairs identified using four scoring approaches on 3 k 10 × PBMCs dataset.

Figure 5

Significant cell–cell interactions of heat‐shocked root cells in Arabidopsis thaliana. (a) UMAP plot of GSE121619 dataset with cells coloured by atrichoblast, cortex, endodermis, lateral root, meristem, pericycle, phloem, trichoblast, and xylem cell types. (b) The mean expression of signature genes for each cell type annotated by MAESTRO software. Low and high gene expression levels are showed by blue and read, respectively. (c) Preference of each cell type under heat‐shock stress. R_O/E above 1 indicates enrichment. (d) Chord diagram of cell–cell communication between pairwise cell types. The line width indicates the number of significant ligand‐receptor pairs. (e) Top 10 ligand‐receptor pairs with P‐value <0.05 show different regulatory pattern. Columns are scaled by max ligand‐receptor expression.

Figure 6

Comparison of the number of interactions of cell pairs between two rice cultivar datasets (93–11 and Nipponbare). (a) UMAP visualization of GSE146035 dataset, including 10 968 cells and 12 564 cells from two cultivar Nipponbare (Japonica) and 93–11 (Indica), respectively. Each dot represents one cell. (b) The mean expression of known marker genes for each cell type from two rice cultivars. (c) Difference of cell–cell interactions of each cell type in two rice cultivars, accounting for total cell number. (d) Identification of significant ligand‐receptor pairs between pairwise cell types in rice cultivar 93–11. (e) Identification of significant ligand‐receptor pairs between pairwise cell types in rice cultivar Nipponbare.