Modelling cell type-specific lncRNA regulatory network in autism with Cycle

Abstract

Background

Autism spectrum disorder (ASD) is a class of complex neurodevelopment disorders with high genetic heterogeneity. Long non-coding RNAs (lncRNAs) are vital regulators that perform specific functions within diverse cell types and play pivotal roles in neurological diseases including ASD. Therefore, exploring lncRNA regulation would contribute to deciphering ASD molecular mechanisms. Existing computational methods utilize bulk transcriptomics data to identify lncRNA regulation in all of samples, which could reveal the commonalities of lncRNA regulation in ASD, but ignore the specificity of lncRNA regulation across various cell types.

Results

Here, we present Cycle (Cell type-specific lncRNA regulatory network) to construct the landscape of cell type-specific lncRNA regulation in ASD. We have found that each ASD cell type is unique in lncRNA regulation, and more than one-third and all cell type-specific lncRNA regulatory networks are characterized as scale-free and small-world, respectively. Across 17 ASD cell types, we have discovered 19 rewired and 11 stable modules, along with eight rewired and three stable hubs within the constructed cell type-specific lncRNA regulatory networks. Enrichment analysis reveals that the discovered rewired and stable modules and hubs are closely related to ASD. Furthermore, more similar ASD cell types tend to be connected with higher strength in the constructed cell similarity network. Finally, the comparison results demonstrate that Cycle is a potential method for uncovering cell type-specific lncRNA regulation.

Conclusion

Overall, these results illustrate that Cycle is a promising method to model the landscape of cell type-specific lncRNA regulation, and provides insights into understanding the heterogeneity of lncRNA regulation between various ASD cell types.

Background

Autism spectrum disorder (ASD) refers to a collection of neurodevelopmental disorders exhibiting profound genetic diversity and complexity [1, 2]. Since childhood, ASD individuals have a wide range of difficulties and deficiencies in social interaction and language communication [3]. Despite striking progress in studying ASD has demonstrated that ASD possesses strong genetic heterogeneity and numerous molecules participate in regulating a series of complex biological processes, including neuronal activity [4] and immune response [2], an understanding of the pathobiology of ASD is still largely unclear. Unlocking the underlying pathogenesis of ASD at the molecular regulatory level holds profound implications in early detection and personalized treatment.

Long non-coding RNAs (lncRNAs) comprise a category of non-coding RNAs that are typically longer than 200 nucleotides, which act as regulators to make significant contributions to neurological diseases, e.g. ASD [3, 5]. In the field of neurobiology, previous studies [3, 6] have revealed that numerous lncRNAs exert biological functions specific to cell types, including neuronal differentiation, synaptic development, and plasticity [7]. In addition, lncRNA regulation also exhibits to be tissue-specific [8], and cell developmental-stage specific [9]. Due to the heterogeneity and complexity in the development of ASD, studying cell type-specific or dynamic lncRNA regulation could provide a new perspective for discovering potential therapeutic strategies for ASD.

Recently, computational methods are promising ways to decipher the function of lncRNAs in modulating ASD-related biological processes. By using bulk transcriptomics data, computational methods for identifying lncRNA regulation can be grouped into three primary categories: sequence-based methods that rely on nucleic acid sequence characteristics, expression-based methods focusing on variations in lncRNA expression levels, and integration-based methods that combine multiple sources of data. Sequence-based methods calculate the binding energy of RNA base pairs to infer lncRNA-target pairs. A prime example is LncTar [10], which utilizes the nearest neighbour thermodynamic model to compute the binding free energy of lncRNA-RNA pairs. Expression-based methods encompass a diverse range of approaches, including statistical methods [11, 12], deep learning methods [13, 14], and causal inference approaches [15]. These methods utilize gene expression profiles to derive and establish lncRNA-target correlation or causality pairs. Alternatively, integration-based methods [16, 17] combine different types of data (e.g., sequence information and expression profiles), thereby enhancing the precision and reliability of lncRNA target prediction. The major limitation of the above methods using bulk transcriptomics data is that they ignore the heterogeneity of lncRNA regulation across various samples (cell lines or tissues). As single-cell and single-nucleus RNA sequencing technology continues to evolve, inferring lncRNA regulation with single-cell or cell type resolution opens a way to explore lncRNA regulation specific to unique cells or cell types in ASD. Regarding cell-specific gene regulation, CSN (Cell-Specific Network) method [18] pioneers the construction of cell-specific networks using single-cell transcriptome data. Subsequently, as an improvement of CSN, c-CSN [19], loc-CSN [20], and p-CSN [21] are also presented to infer conditional, local, and partial cell-specific networks, respectively. Specifically, for exploring cell-specific miRNA regulation, CSmiR [22] has also been developed to investigate single-cell level modulation of miRNA expression. In terms of regulation specific to individual cell types, scHumanNet [23] aims to generate specialized gene regulatory networks (GRNs) for individual cell types by leveraging the information contained in the HumanNet reference interactome [24] and single-cell expression data. Recently, scMTNI [25] integrates single-cell multi-omics datasets to build GRNs specific to cell types across cell lineages. However, these cell-specific or cell type-specific regulation approaches primarily prioritize transcription factor or miRNA regulation over lncRNA regulation. To infer lncRNA regulation specific to biological conditions, CDSlncR [9] could infer lncRNA regulatory networks corresponding to distinct developmental states of the brain neocortex. Given that the pathogenesis of ASD involves a series of cell types and biological processes regulated by lncRNAs, thus it is crucial to study cell type-specific lncRNA regulation in ASD.

To explore the dynamic lncRNA regulation across various ASD cell types, we present a novel method, Cycle (Cell type-specific lncRNA regulatory network), to model cell type-specific lncRNA regulatory networks in ASD. Cycle has two main contributions as follows. Firstly, instead of considering all types of interactions among genes (lncRNAs and mRNAs), Cycle concentrates on identifying interactions between lncRNAs and mRNAs. Secondly, taking the diversity and specificity of cells and cell types into consideration, Cycle identifies lncRNA regulatory networks specific to each cell type.

We have applied Cycle into single-nucleus RNA-sequencing (snRNA-seq) data of ASD brain tissues [26] for modelling the landscape of cell type-specific lncRNA regulation in ASD. Our research has found that each ASD cell type is unique in lncRNA regulation. Notably, over one-third of the cell type-specific lncRNA regulatory networks are scale-free, and all of them exhibit to be small-world. Among 17 ASD cell types, we have inferred 19 rewired modules and 11 stable modules, along with eight rewired hubs and three stable hubs based on the identified cell type-specific lncRNA regulatory networks. Importantly, the discovered rewired and stable modules, and stable hubs are closely associated with ASD. Additionally, ASD cell types that are more similar tend to be strongly connected in the constructed cell similarity network. Finally, our comparison results suggest that Cycle is a promising approach in elucidating cell type-specific lncRNA regulation.

Methods

The flowchart of Cycle

Cycle includes three main components (Fig. 1). Firstly, Cycle conducts data pre-processing, including gene annotation, feature selection, and data splitting to acquire highly expressed lncRNAs and mRNAs in 17 ASD cell types. For each ASD cell type, Cycle further identifies lncRNA regulatory networks specific to it. In total, 17 cell type-specific lncRNA regulatory networks are modelled. Based on the constructed cell type-specific lncRNA regulatory networks, Cycle further infers the rewired and stable modules and hubs. Finally, Cycle performs four types of downstream analyses, including network topological analysis, uniqueness analysis, cell similarity network construction, and enrichment analysis. The details of each component will be described in the following.

Fig. 1

Workflow of Cycle. Firstly, Cycle extracts the matched lncRNA and mRNA expression data by using gene annotation information from HGNC (HUGO Gene Nomenclature Committee), and further retains the highly expressed lncRNAs and mRNAs for each cell type. In total, we have obtained 17 cell type-specific expression data of highly expressed lncRNAs and mRNAs. Secondly, Cycle models cell type-specific lncRNA regulatory networks for 17 ASD cell types. Furthermore, Cycle identifies the rewired and stable modules, and hub lncRNAs by utilizing the constructed lncRNA regulatory networks. Finally, Cycle conducts four types of downstream analyses, including network topological analysis, uniqueness analysis, cell similarity network construction, and enrichment analysis. Created with BioRender.com

Single-nucleus RNA-sequencing data in ASD

The raw snRNA-seq data of ASD [26] is obtained from the Sequence Read Archive (SRA) with accession number PRJNA434002 (https://ncbi.nlm.nih.gov/bioproject/434002), and the analyzed data is from https://autism.cells.ucsc.edu. As a pre-processing step, we obtain the matched lncRNA and mRNA expression data by utilizing gene annotation information from HGNC (HUGO Gene Nomenclature Committee, https://www.genenames.org/), and select lncRNAs and mRNAs whose expression levels are higher than the average expression level across all cells. In total, we have retained 813 lncRNAs and 5,133 mRNAs highly expressed in 52,003 ASD cells. Similar to [26], we annotate the 52,003 ASD cells into 17 cell types based on the expression of known cell type markers. These 17 cell types include oligodendrocyte precursor cells (OPC), oligodendrocytes, microglia cells, fibrous astrocytes (ASTFB), protoplasmic astrocytes (ASTPP), layer 2/3 excitatory neurons (L2/3), layer 4 excitatory neurons (L4), layer 5/6 corticofugal projection neurons (L5/6), layer 5/6 cortico-cortical projection neurons (L5/6-CC), SV2C interneurons (IN-SV2C), somatostatin interneurons (IN-SST), VIP interneurons (IN-VIP), parvalbumin interneurons (IN-PV), endothelial cells, NRGN-expressing neurons (NeuNRGN-I), NRGN-expressing neurons (NeuNRGN-II), and maturing neurons (Neu-mat). Various cell types play different roles in the brain region (see Additional file 1 for the detailed information of 17 ASD cell types).

Identification of cell type-specific lncRNA regulatory networks

For each cell type, modelling cell type-specific networks is grounded upon the identified cell-specific regulatory networks. Hence, Cycle is firstly used to identify cell-specific lncRNA regulatory networks. Here, Cycle adapts CSN [18] with local strategy [20] to quantitatively estimate the correlation strength of lncRNA-mRNA relationship pairs in each cell.

In cell k, $l_{k}$ and $m_{k}$ are the expression values of lncRNA $lncR_{l}$ and mRNA $mR_{m}$, respectively, and $\rho_{lm}^{\left( k \right)}$ is calculated as the interaction strength between $lncR_{l}$ and $mR_{m}$ in the following:

$$\rho_{lm}^{\left( k \right)} = \frac{{n_{lm}^{\left( k \right)} }}{N} - \frac{{n_{l}^{\left( k \right)} }}{N} \cdot \frac{{n_{m}^{\left( k \right)} }}{N}$$

(1)

where N is the number of cells for ASD snRNA-seq data, $n_{l}^{\left( k \right)}$ and $n_{m}^{\left( k \right)}$ are the neighbourhood number of $l_{k}$ and $m_{k}$ in the bins of cell k for $lncR_{l}$ and $mR_{m}$, respectively, and $n_{lm}^{\left( k \right)}$ is the neighbourhood number of $\left( {l_{k} ,m_{k} } \right)$ in the interaction bin of cell k.

Owing to the specificity and heterogeneity of cells, self-adaptive window size $B_{l}^{\left( k \right)}$ and $B_{m}^{\left( k \right)}$ of bins in cell $k$ are iteratively generated based on local standard deviations as follows:

$$B_{l}^{\left( k \right)} \left( t \right) = \left\{ {i:\left| {l_{i} - l_{k} } \right| \le St. Dev\left\{ {l_{i} :i \in B_{m}^{\left( k \right)} \left( {t - 1} \right)} \right\}} \right\}$$

(2)

$$B_{m}^{\left( k \right)} \left( t \right) = \left\{ {i:\left| {m_{i} - m_{k} } \right| \le St. Dev\left\{ {m_{i} :i \in B_{l}^{\left( k \right)} \left( {t - 1} \right)} \right\}} \right\}$$

(3)

where initial $B_{l} \left( 0 \right)$ and $B_{m} \left( 0 \right)$ are a quantile range, and $t$ starts from $1,2,...$, until convergence is achieved. If the convergence is not achieved during the iterations, $B_{l}^{\left( k \right)} \left( 1 \right)$ and $B_{m}^{\left( k \right)} \left( 1 \right)$ will be adopted as window sizes in practice for cell $k$.

$$z_{lm}^{\left( k \right)} = \frac{{\rho_{lm}^{\left( k \right)} - \mu_{lm}^{\left( k \right)} }}{{\sigma_{lm}^{\left( k \right)} }} = \frac{{\sqrt {n - 1} \cdot \left( {n \cdot n_{lm}^{\left( k \right)} - n_{l}^{\left( k \right)} n_{m}^{\left( k \right)} } \right)}}{{\sqrt {n_{l}^{\left( k \right)} n_{m}^{\left( k \right)} \left( {n - n_{l}^{\left( k \right)} } \right)\left( {n - n_{m}^{\left( k \right)} } \right)} }}$$

(4)

where $\mu_{lm}^{\left( k \right)} = 0$ represents the mean value of $\rho_{lm}^{\left( k \right)}$, and $\sigma_{lm}^{\left( k \right)} = \sqrt {\frac{{n_{l}^{\left( k \right)} n_{m}^{\left( k \right)} \left( {n - n_{l}^{\left( k \right)} } \right)\left( {n - n_{m}^{\left( k \right)} } \right)}}{{n^{4} \left( {n - 1} \right)}}}$ denotes the standard deviation of $\rho_{lm}^{\left( k \right)}$. Each $z_{lm}^{\left( k \right)}$ value has a corresponding p-value, and a smaller p value (e.g., p < 0.01) indicates a higher strength between $lncR_{l}$ and $mR_{m}$ in cell k.

For each cell, we only focus on the lncRNA-mRNA interactions with statistical significance (e.g., p-value less than 0.01). If a significant lncRNA-mRNA interaction exists in more than 90% of total cells of a cell type, the lncRNA-mRNA interaction is regarded as one of a collection of lncRNA-mRNA interactions for the cell type. By integrating all of the lncRNA-mRNA interactions belonging to individual cell types, Cycle constructs 17 cell type-specific lncRNA-mRNA regulatory networks.

Network topological analysis

Topological analysis contributes to exploring the characteristics and organization of biological networks including lncRNA regulatory networks. Degree and density are two widely used metrics to characterize a biological network. If the node degree distribution of a cell type-specific lncRNA regulatory network adheres to a power-law distribution with a p value of a Kolmogorov–Smirnov test [27] larger than 0.05, the network tends to be a scale-free network. If the characteristic density of a cell type-specific lncRNA regulatory network is higher than that of its corresponding random networks at a significance level (e.g., 0.05), the network is regarded as a small-world network. Here, for each cell type-specific lncRNA regulatory network, we generate 1,000 random networks by randomizing the lncRNA-mRNA interactions. We utilize the Student's t-test for statistically quantifying the differences between the constructed cell type-specific lncRNA regulatory networks and their corresponding random networks. In this work, the igraph R package [28] is applied to analyze the topological attributes of the constructed cell type-specific lncRNA regulatory networks.

Hub lncRNA inference

Hub lncRNAs with high connectivity play key pivot roles in a cell type-specific lncRNA regulatory network. Rather than inferring hubs as those with a node degree exceeding a giving value, we assume that the node degree of lncRNAs follows the Poisson distribution [29,30,31]. For each lncRNA, we calculate the p value of it as follows:

$$p\left( {x \ge k} \right) = 1 - p\left( {x < k} \right) = 1 - \mathop \sum \limits_{i = 0}^{k - 1} \frac{{e^{ - \lambda } \lambda^{i} }}{i!}$$

(5)

where $\lambda = np$, $p = \frac{m}{{A_{n}^{2} }}$, $n$ is the number of lncRNAs, m is the number of lncRNA-mRNA pairs in a lncRNA–mRNA regulatory network, and $A_{n}^{2}$ is the number of all possible lncRNA-mRNA interactions. In this work, a lncRNA with a p value less than 0.05 is considered as a hub lncRNA.

Stable and dynamic analysis

We perform stable and dynamic analysis to reveal the commonality and heterogeneity of lncRNA regulation between different ASD cell types. Previous studies [15, 32] have shown that lncRNA regulation is 'on' in some biological conditions but is 'off' in other biological conditions. Here, lncRNA-mRNA interactions or hub lncRNAs existing in at least 90% of ASD cell types are considered as the stable lncRNA regulatory network or hub lncRNAs, and lncRNA-mRNA interactions or hub lncRNAs existing in only one ASD cell type are viewed as the rewired lncRNA regulatory network or hub lncRNAs. To further identify highly connected functional modules within the stable and rewired lncRNA regulatory networks, we have applied the Markov Cluster (MCL) algorithm [33] to discover the stable and rewired lncRNA regulatory modules. In each module, the number of genes (lncRNAs and mRNAs) should be at least three.

Uniqueness of cell type-specific lncRNA regulation

To know whether lncRNA regulation in each ASD cell type is unique or not, we calculate the differences of lncRNA-mRNA regulatory networks (or hub lncRNAs) across ASD cell types. Before computing differences between each pair of ASD cell types, we need to calculate the similarity of lncRNA-mRNA regulatory networks (or hub lncRNAs) across ASD cell types. To assess the similarity of lncRNA-mRNA regulatory networks (or hub lncRNAs) between two ASD cell types, we employ the Simpson model [34] to analyze lncRNA-mRNA interactions (or hub lncRNAs) obtained from ASD cell types $i$ and $j$, respectively. This approach generates a similarity value $sim_{ij}$ that measures the similarity between cell type-specific lncRNA regulatory networks (or hub lncRNAs). The dissimilarity or difference $dif_{ij}$ between cell type-specific lncRNA regulatory networks (or hub lncRNAs) is then computed as described below.

$$dif_{ij} = 1 - sim_{ij}$$

(6)

$$sim_{ij} = \frac{{\left| {LR_{i} } \right| \cap \left| {LR_{j} } \right|}}{{min\left( {\left| {LR_{i} } \right|,\left| {LR_{j} } \right|} \right)}}$$

(7)

where $LR_{i}$ and $LR_{j}$ are lncRNA-mRNA interactions or hub lncRNAs existing in ASD cell types $i$ and $j$, $\left| {LR_{i} } \right| \cap \left| {LR_{j} } \right|$ represents the intersection number of lncRNA-mRNA interactions or hub lncRNAs between $LR_{i}$ and $LR_{j}$, and $min\left( {\left| {LR_{i} } \right|,\left| {LR_{j} } \right|} \right)$ is the smaller number of lncRNA-mRNA interactions or hub lncRNAs between $LR_{i}$ and $LR_{j}$. A larger value of $dif_{ij}$ denotes a higher uniqueness between ASD cell types $i$ and $j$.

Enrichment analysis

To understand the potential biological functions of the stable and rewired lncRNA-mRNA regulatory modules, we conduct enrichment analysis with miRspongeR [35] and clusterProfiler [36] R packages. The databases used for functional enrichment analysis include Gene Ontology (GO) [37], Kyoto Encyclopedia of Genes and Genomes (KEGG) [38], and Reactome Pathway database (Reactome) [39]. Additionally, three disease databases including Disease Ontology (DO) [40], DisGeNET [41], and Network of Cancer Genes (NCG) [42] are also considered for disease enrichment analysis. With regard to hub lncRNAs, we employ RNAenrich [43], a powerful comprehensive web server for ncRNA enrichment, to explore potential pathways, biological processes, and diseases in which they participate. In this work, the enriched KEGG, GO, Reactome, DO, DisGeNET, or NCG term with an adjusted p value < 0.05 (adjusted by the Benjamini–Hochberg approach) is viewed as a significant term.

Results

In this section, we show the application of Cycle in uncovering cell type-specific lncRNA regulation in ASD. The ASD dataset used and code for the reproducibility of the analysis are available at https://github.com/chenchenxiong/Cycle.

The landscape of cell type-specific lncRNA regulation in ASD

To investigate the lncRNA regulation across 17 ASD cell types, we have identified the landscape of cell type-specific lncRNA regulation by following the workflow of Cycle (Fig. 1). We have discovered that the number of lncRNA-mRNA interactions and hub lncRNAs tends to be various across 17 ASD cell types (Fig. 2a). In the case of lncRNA-mRNA interactions, L4 and microglia cells obtain the largest and least number of interactions, respectively. In the case of hub lncRNAs, L2/3 and ASTFB have the largest and least number of hubs, respectively. Network topological analysis displays that 6 out of 17 (~ 35.29%) cell type-specific lncRNA-mRNA regulatory networks adhere to a power law distribution, and all cell type-specific lncRNA-mRNA networks display higher densities compared to their corresponding random networks (Fig. 2b and Additional file 2). These results indicate that over one-third of these cell type-specific lncRNA regulatory networks tend to be scale-free, and all of these cell type-specific lncRNA regulatory networks exhibit to be small-world.

Fig. 2

The lncRNA regulation landscape across 17 cell types. a The number of lncRNA-mRNA interactions and hub lncRNAs within each cell type. b The topological properties of lncRNA-mRNA regulatory networks for each cell type

Each ASD cell type is unique in lncRNA regulation

To understand whether the identified lncRNA regulation for each ASD cell type is unique, we calculate the dissimilarity or difference between each pair of ASD cell types. We have found that the identified lncRNA regulatory networks and hub lncRNAs between any pairs of 17 ASD cell types are various, indicating the uniqueness of each cell type (Fig. 3a and b). For the identified lncRNA regulatory networks, nearly half of pairs between 17 cell types (~ 46.32%) have a difference value with more than 0.500. Specifically, the highest difference value (0.951) between 17 cell types is between microglia and NeuNRGN-II (Fig. 3a). For hub lncRNAs, more than one-third pairs between 17 cell types (~ 36.76%) have a difference value with more than 0.500, and the highest difference value (0.844) between 17 cell types also exists between L5/6-CC and Neu-mat (Fig. 3b). These results have suggested that each ASD cell type is unique in lncRNA regulatory networks and hub lncRNAs.

Fig. 3

Uniqueness of lncRNA regulation across ASD cell types. a, b) Uniqueness of lncRNA-mRNA regulatory networks and hub lncRNAs in each ASD cell type. c The radar chart of lncRNA-mRNA interactions and hub lncRNAs

In addition, we have further explored stable and rewired lncRNA-mRNA regulatory networks and hubs across various cell types. As a result, we have found that 115,108 lncRNA-mRNA interactions and eight hub lncRNAs (CPVL-AS2, LINC00343, LINC01202, LINC01619, LINC01811, LINC02301, LINC03013, and SYNPO2L-AS1) only exist in one cell type, and 63 lncRNA-mRNA interactions and three hub lncRNAs (ANKRD17-DT, LINC01572, and MIRLET7BHG) exist in at least 90% cell types (Fig. 3c). In total, we have obtained 115,108 rewired interactions, 63 stable interactions, eight rewired hubs, and three stable hubs across 17 ASD cell types. Overall, the number of rewired interactions or hubs is larger than that of stable interactions or hubs, indicating that lncRNA regulation tends to be dynamic across ASD cell types.

Rewired and stable lncRNA regulatory modules and hub lncRNAs are closely associated with ASD

To identify highly connected functional modules within the identified stable and rewired lncRNA regulatory networks, we have discovered 19 rewired and 11 stable lncRNA regulatory modules by applying MCL algorithm [33]. To reveal the biological significance of the identified rewired and stable modules, we further conduct functional and disease enrichment analysis of them. Enrichment analysis results indicate that most of the rewired and stable lncRNA regulatory modules are functional and significantly enriched in ASD-related biological processes and pathways (Tables 1 and 2, Additional file 3). For example, a large number of enriched terms, including neurodevelopment [44], synaptic transmission [45], cellular signalling [46], and immune response [47] are closely related to the pathogenesis of ASD.

Table 1 Enrichment analysis of rewired lncRNA regulatory modules related to ASD

Table 2 Enrichment analysis of stable lncRNA regulatory modules related to ASD

Based on the rewired and stable lncRNA-mRNA regulatory networks, we further infer eight rewired hubs and three stable hubs, and conduct enrichment analysis of them with RNAenrich [43]. Enrichment analysis reveals that the rewired hubs are not significantly enriched in any functional terms, but the stable hub lncRNAs are significantly enriched in 184 KEGG, 3074 GO, and 364 Reactome terms. Particularly, several GO, KEGG or Reactome terms are closely associated with ASD (Table 3, Additional file 3).

Table 3 Enrichment analysis of stable hub lncRNAs related to ASD

Cell similarity network

To understand the similarity of each pair of ASD cell types, we have further constructed a cell similarity network by using the inferred cell type-specific lncRNA-mRNA interactions and hub lncRNAs in each ASD cell type. In this work, if the similarity value of a cell–cell pair is larger than the median value of similarity, the cell–cell pair is considered to be a link in the cell similarity network. As a result, we have found that L4 is similar with the largest number of other ASD cell types, while ASTFB is similar with the least number of other ASD cell types (Fig. 4).

Fig. 4

Cell similarity network among 17 ASD cell types by using lncRNA-mRNA interactions and hub lncRNAs. A larger circle indicates that the cell type is similar to a larger number of other cell types

In comparison with the other method

In this section, we compare Cycle with the other three methods, including CDSlncR [9], LncRNA2Target v3.0 [82], and WGCNA [83] in the identification of cell type-specific lncRNA regulation. CDSlncR [9] is the first method to investigate cell type-specific lncRNA regulation. LncRNA2Target [82] predicts lncRNA-mRNA interactions by analyzing the knockdown and overexpression of lncRNAs. WGCNA (Weighted Gene Co-expression Network Analysis) [83] is a co-expression-based prediction method that can be utilized to identify lncRNA-mRNA interactions. In the WGCNA analysis, we set the scale-free topology model fit index ($R^{2}$) to 0.85, and filter out lncRNA-mRNA pairs whose topological overlap measure (TOM) similarity values are greater than the median TOM similarity value of all lncRNA-mRNA pairs within each cell type. To ensure fairness, the p value cutoff of CDSlncR, LncRNA2Target, and Cycle is set to be the same. We focus on comparing the number of validated lncRNA-mRNA interactions predicted by Cycle, CDSlncR [9], LncRNA2Target [82], and WGCNA [83]. These experimentally validated lncRNA-mRNA interactions are collected from NPInter v5.0 [84], LncTarD v2.0 [85] and LncRNA2Target v3.0 [82]. For each ASD cell type, Cycle performs the best due to obtaining the largest number of validated lncRNA-mRNA interactions (Fig. 5). The comparison result indicates that Cycle is a promising approach to model cell type-specific lncRNA regulation in ASD.

Fig. 5

Comparison results in the number of validated lncRNA-mRNA interactions between Cycle, CDSlncR, LncRNA2Target and WGCNA

Discussion

ASD is a set of complex neurodevelopmental disorders that manifest with varying symptoms among individuals. Exploring the regulatory mechanisms of lncRNAs within and between different ASD cell types helps to elucidate the etiology and ontogeny of ASD. In this work, we develop a novel approach (Cycle) to model cell type-specific lncRNA regulatory networks in ASD. For each ASD cell type, we have shown that the lncRNA regulation tends to be unique. The rewired and stable lncRNA regulatory modules and hub lncRNAs are significantly enriched in several ASD-related terms or pathways. In addition, cell similarity networks can help to know which cell types are similar to the least or largest number of other ASD cell types. In comparison with other methods, Cycle performs the best in inferring cell type-specific lncRNA regulation.

In future, Cycle can be further improved in the following four aspects. Firstly, Cycle currently focuses on lncRNA regulation specific to ASD cell types. In future, it is necessary to study other condition-specific lncRNA regulation, e.g., sex-specific or region-specific lncRNA regulation. Secondly, Cycle only infers the association/correlation rather than causal relationships between lncRNAs and mRNAs. In future, we will conduct cell type-specific lncRNA causal regulation research. Thirdly, the enrichment analysis highly depends on the incomplete annotated databases and prior knowledge. Consequently, the enrichment analysis results are biased, which is a common issue of existing computational methods, including Cycle. To alleviate this problem, it is necessary to integrate more functional annotation information in future. Finally, the competing endogenous RNA (ceRNA) hypothesis [86] suggests that lncRNAs can modulate gene expression by acting as ceRNAs, thus it is strongly needed to identify cell type-specific lncRNA-related ceRNA networks for comprehensively understanding lncRNA regulation.

Conclusion

Taken together, Cycle is useful for modelling the landscape of cell types-specific lncRNA regulation in ASD and contributes to elucidating the heterogeneity of lncRNA regulation underlying various ASD cell types.