We first conducted a systematic literature search for murine in vivo EC bulk RNAseq studies in PubMed, the NCBI GEO database, and the ArrayExpress database. It resulted in 19 RNA studies for EC under various conditions. They include 71 EC conditions. Each condition has multiple replicated samples, and in total, there are 203 samples. The raw sequence data for each condition, including the raw data for its exact control group, were obtained from the NCBI Short Read Archive (SRA) or ArrayExpress database.

The raw sequence data obtained from SRA and ArrayExpress were preprocessed with the Galaxy online server (Jalili et al., 2020) (https://usegalaxy.eu/, version: 20.09) using a standardized procedure for all datasets. The detailed procedure is described in the Galaxy RNA-seq analysis instruction (https://training.galaxyproject.org/training-material/topics/transcriptomics/tutorials/rna-seq-reads-to-counts/tutorial.html).

The sequence data were uploaded in two ways: for the data available in SRA, the SRA-tools (Leinonen et al., 2011) (version: 2.10.8) in Galaxy were used to upload these datasets reads in the FASTA/Q format from the NCBI; for the other datasets from ArrayExpress, the ArrayExpress FTP download links were used. The FASTQ sequence files were then aligned to the referenced mouse genome assembly (GRCm38/mm10) obtained from the UCSC Genome Browser database (Navarro Gonzalez et al., 2021) using the HISAT2 tool (Kim et al., 2015) (version: 2.1.0) on Galaxy. The gene annotation file GTF (2020, ncbiRefSeq, mm10) was also obtained from the UCSC Genome Browser database, which was consistent with the genome sequence file. The alignment bam files were then input to the featureCounts tool (Liao et al., 2014) (version: 2.0.1, with default parameters) to get the raw read counts for each genes (feature count files). In total, 203 samples were quantified and their count data were processed in R (version: 4.0.3) for downstream analysis.

In order to compare the EC expression level among different samples in different conditions, all the raw count data were normalized using rpkm function in the edgeR package (version: 3.32.0). The FPKM values for each sample were calculated, and then, the average expressions and standard deviations for each of the 32 conditions (71 bars in the FPKM plot) were calculated in R. The result for each gene was visualized in bar-plot using the ggplot2 package (version: 3.3.2).

The gene expression raw count files were imported into the limma package (version: 3.46.0) in R, and the voom function was used to compare the gene expression between two groups (treated versus control) with the default parameters. To remove low-expression genes in each sample, the genes which were detected in only one sample were filtered out. To visualize the differential expression profiles among the 40 comparison groups, the fold changes and the standard deviations for each gene were visualized in bar-plots.

To search for the genes with similar expression profiles with a query gene, the corr.test function from psych package (version: 2.0.12) was applied. The correlation coefficient and the p values were calculated. The sorted result was stored in a table and is available for download through our ECO database. In addition, to better illustrate the correlation result, we chose the 10 most correlated genes to the query gene and generated a heatmap with the pheatmap package (version: 1.0.12).

Our ECO database, an interactive web application, is built mainly using the R Shiny package (version: 1.6.0), as well as the other auxiliary packages including shinythemes (version: 1.2.0), ggplot2 (version: 3.3.3), and ggh4x (version: 0.1.2.1). The ECO database is available for free at https://heomics.shinyapps.io/ecodb/.

In order to construct a comprehensive omnibus of mouse in vivo EC RNAseq profiles, we performed literature mining and identified 19 currently available RNA studies (Supplementary Table S1), which cover EC in a variety of pathological alterations, genetic modifications, and other stimulated conditions. In these studies, freshly isolated ECs were analyzed with RNA sequencing. In total, there are 203 samples covering 71 in vivo conditions from 10 organs. These data composite the base for ECO, and they were processed as shown in the workflow (Figure 1). First, their raw data were obtained from the GEO database or ArrayExpress database, respectively. The sequence data were aligned to a standardized mouse genome assembly (GRCm38/mm10), and gene expression in each sample was quantified using the Galaxy analysis platform (Jalili et al., 2020). The gene expression levels in each condition were then summarized (average FPKM and standard deviation) and available for bar-plot visualization in the ECO database (https://heomics.shinyapps.io/ecodb). Also, the gene expression in each condition was compared with its respective control by differential expression analysis, and log scaled fold change (logFC) and p values were calculated, which are also illustrated with the bar-plot in the database. Besides the display of expression profiles in ECs, ECO can further identify the genes which showed similar expression profiles with the queried gene by using correlation analysis. The results were shown both as a heatmap and table.

ECs in different tissues have heterogeneous phenotypes for their distinct physiological needs (Kalucka et al., 2020). For instance, brain ECs form tight junctions and express active transporters to restrict diffusion, known as the blood–brain barrier (BBB) (Daneman and Prat, 2015). In contrast, ECs in the kidney are associated with fenestrae to allow efficient passage of high-volume fluids and formation of urine (Dumas et al., 2021). EC profiles from 10 organs, including the brain, lung, bone, kidney aorta, liver, eye, muscles, lymph node, and embryo, were cataloged in ECO. Users can access and download the expression of the gene of their interest in ECs of different organs in ECO by simple one-click of FPKM button. Also, users can input a customized gene list to analyze their overall gene expression enrichment pattern in a heatmap.

We use Slc2a1 as an example to explore the inter-organ heterogeneity of a given gene. Slc2a1, encoding Glut1, which is highly expressed in BBB ECs but not peripheral ECs and facilitates glucose transport over BBB (Zheng et al., 2010). When we access Slc2a1 expression by pressing the FPKM button after entering the gene symbol in the query interface, we get the normalized data bar-plot visualization for 32 sub-groups from 10 organs. As expected, Slc2a1 is highly expressed in ECs from the brain, but almost absent in other organs (Figure 2).

ECs participate in the regulation of multiple processes including angiogenesis, coagulation, and inflammation. Endothelial dysfunction is associated with many pathological alterations and aggravates progression of multiple life-threatening diseases including cancers, cardiovascular disease, diabetes mellitus, and renal disorders. In ECO, we collected EC transcriptomes from eleven mice disease models (cerebral cavernous malformation (CCM), epilepsy, experimental autoimmune encephalomyelitis (EAE), stroke, traumatic brain injury (TBI), diabetic nephropathy, Alport syndrome, liver cancer, non-alcoholic steatohepatitis (NASH), experimental autoimmune uveitis (EAU), and facioscapulohumeral muscular dystrophy (FSHD)), seven gene-modified animal models (Jnk1/2/3 EC-specific deficient, Cpt1a EC-specific deficient, Tsc2 mesenchyme cell-specific deficient, Zmpste24 deficient, adrenomedullin (AD) EC-specific deficient, Tankl stroma cell-deficient, and EC-specific Notch1 mutants), and two clinically relevant treatments (VEGF stimulations and chemo/radiotreatment) (Supplementary Table S1). The users can access the alteration of the genes of their interest in response to the abovementioned conditions compared to their control by clicking the logFC button. The result is illustrated in a bar-plot with 40 columns; each column represents the log2 scaled fold change, and its statistical significance (p value range) is indicated by asterisks (Figure 3).

We use Sele as an example to demonstrate the exploration of its regulation in different pathological conditions, genetic modifications, and treatments in vivo. E-selectin, encoded by Sele, is upregulated in ECs in response to pro-inflammatory signals, promoting the rolling and adherence of immune cells to ECs for their diapedesis (Jubeli et al., 2012). Inflammation is closely linked in the EC dysfunction in multiple diseases (Steyers and Miller, 2014). As shown in Figure 3, ECO provides a comprehensive portrait for Sele in different pathological conditions. Sele was upregulated in ECs from eight disease models including CCM, EAE, stroke, TBI, epilepsy, AOD, diabetic nephropathy, Alport syndrome, and NASH, highlighting the broad role of Sele in multiple disease progressions (Silva et al., 2017).

The correlation analysis identifies the genes which have similar expression profiles, and those co-expressed genes may have similar function. In ECO, it provides correlation analysis for the query gene to all other genes in all cataloged EC groups, as well as in individual organs which have relatively large number of samples. This facilitates uncovering the function of novel or not much characterized genes based on correlation analysis.

For example, we used a gene named C330027C09Rik as an example. C330027C09Rik did not yet have a clear gene name at the time of the gene assembly from the Ensembl database and was named after the full-length cDNA sequences from the RIKEN project (Hayashizaki, 2003). Among the top correlated genes, a list of well-known cell cycle-related genes appears, for example, Mki67 and Cenpf, indicating that this gene maybe related with cell cycle (Figure 4). Interestingly, in the NCBI gene database, C330027C09Rik has been formally named as cell proliferation-regulating inhibitor of protein phosphatase 2A (Cip2a) (https://www.ncbi.nlm.nih.gov/gene/?term=C330027C09Rik). This confirmed the prediction from the correlation analysis.