Non-negative matrix factorization can factorize a high-dimensional gene expression matrix into two low-dimensional matrices, i.e., a gene-factor matrix and a sample-factor matrix, achieving the purpose of dimension reduction. To integrate and factorize multiple gene expression datasets into a common gene-factor matrix, multiple dataset-specific gene-factor matrices, and sample-factor matrices, iNMF (Yang and Michailidis, 2016) was proposed. The optimization problem is: min W , V 1 , … , V k , H 1 , … H k ∑ i = 1 k ‖ X i − H i ( W + V i ) ‖ F 2 + λ ∑ i = 1 k ‖ H i V i ‖ F 2 , s .t .     W ≥ 0 , V i ≥ 0 , H i ≥ 0 ,   i = 1,2 , … k , where X i ∈ R n i × g denotes each gene expression dataset, g denotes the number of genes, and n i denotes the number of samples in the ith dataset. X i is factorized into three low-dimensional matrices, H i ∈ R n i × m , W ∈ R m × g , and V i ∈ R m × g , where m denotes the number of factors/gene modules. H i is the representation of samples in the low-dimensional space. V i and W are the dataset-specific and shared gene modules, respectively. λ is a regularization parameter.

The regularization of iNMF can make V i sparser to some degree, while to facilitate the analyses of shared and dataset-specific gene modules, we propose integrative regularized non-negative matrix factorization (iRNMF) by imposing a new regularization. The optimization problem is: min W , V 1 , … , V k , H 1 , … H k ∑ i = 1 k ‖ X i − H i ( W + V i ) ‖ F 2 + λ   1 ∑ i = 1 k ‖ H i V i ‖ F 2 + λ 2 ∑ i = 1 k ‖ H i W ‖ F 2 , s .t .     W ≥ 0 , V i ≥ 0 , H i ≥ 0 ,   i = 1,2 , … k , where λ   1 and λ 2 are regularization parameters. The multiplicative updates often used for NMF-like optimization problems do not have a convergence guarantee and may need more iterations; therefore, we applied the block coordinate descent algorithm used in LIGER (Welch et al., 2019). We divided the variables into 2k + 1 blocks (corresponding to H i , V i for each dataset, and W) and performed block coordinate descent, iteratively minimizing the objective with respect to each block, holding the others fixed. We iterated: W = a r g m i n W ≥ 0 | | ( H 1 ⋮ H k λ 2 H 1 ⋮ λ 2 H k ) W − ( X 1 − H 1 V 1 ⋮ X k − H k V k O n 1 × g ⋮ O n k × g ) F 2 | | , H i = a r g m i n H i ≥ 0 | | ( W T + V i T λ 1 V i T λ 2 W T ) H i T − ( X i T O g × n i O g × n i ) F 2 | | , V i = a r g m i n V i ≥ 0 | | ( H i λ 1 H i ) V i − ( X i − H i W O n i × g ) F 2 | | , until convergence. Each of the optimization subproblems mentioned previously requires solving a non-negative least-squares problem, and we used the fast block principal pivoting algorithm (Kim et al., 2014) to solve each of these subproblems.

We downloaded the single-nucleus gene expression data derived from the middle temporal gyrus (MTG) of the human cortex (Hodge et al., 2019) from the Allen Institute for Brain Science. It includes 15,928 nuclei sampled from eight human donor brains, of which 15,206 were from postmortem donors with no known neuropsychiatric or neurological conditions and 722 were from distal and normal tissues of neurosurgical donors. We preprocessed the data with R packages of scatter (McCarthy et al., 2017) and scran (Lun et al., 2016), including the quality control of nuclei and genes, and removing a minority of nuclei assigned to different cell cycle phases by the function of cyclone in scran. Nuclear and mitochondrial genes downloaded from Human MitoCarta2.0 (Calvo et al., 2016) were excluded, and protein-coding genes were retained. After removing the nuclei not assigned to any specific cell type, we obtained the expression level of 17,120 protein-coding genes in 12,246 nuclei. Then, we used scran to obtain 7,011 highly variable protein-coding genes across all nuclei, which were defined as genes with biological components that are significantly greater than zero at a false discovery rate (FDR) of 0.1. After removing the cell types containing less than 20 nuclei, we obtained the gene expression data of nuclei from glutamatergic neuron (Gluta), GABAergic interneuron (GABA), astrocyte (Ast), oligodendrocyte (Oli), and oligodendrocyte precursor cell (OPC), including 8994, 2762, 227, 112, and 133 nuclei, respectively. The gene expression of 7,011 highly variable protein-coding genes in these five cell types was used for analyses.

To determine the number of factors/gene modules m, we used the same way with LIGER, applying Kullback–Leibler (KL) divergence as a criterion. When the number of factors is too low, factors will include many genes and samples will load on many factors, with the distribution of factor loadings for a particular sample approaching a uniform distribution (Welch et al., 2019). As the number of factors approaches the true number of gene modules, each sample will generally load on only a few factors. Therefore, we calculated the KL divergence, compared to a uniform distribution, of the factor loadings for each sample and plotted the median across samples as a function of m to select the saturation point of the curve as the optimal m. We also considered the mean squared error (MSE) between X i and the reconstructed data X i ^ , i.e., ∑ i = 1 k 1 / ( ni × g ) | | X i − H i ( W + V i ) | | F 2 , to help to determine the optimal m. To select the regularization parameters λ   1   and λ 2 , we applied the alignment metric (Butler et al., 2018) as a criterion, which LIGER also used, and plotted the alignment metric as a function of a combination of λ   1 and λ 2 to choose the point at which the alignment metric reaches the minimum value. λ   1   and λ 2 can be a value among 0.01, 0.1, 1, 10, and 100.

We used iRNMF to analyze gene expression datasets of multiple cell types derived from human MTG. After obtaining the cell-type-specific and shared gene module matrices V i   and   W , for each gene module, we calculated the z-scores of genes, and the genes whose z-scores are larger than one were regarded as module genes. The modules with no less than 20 module genes were reported. The gene modules significantly enriched with ASD genes were regarded as ASD-associated gene modules. ASD candidate genes were downloaded from the Simons Foundation Autism Research Initiative (SFARI), version of 2 September 2021. We identified ASD-associated gene modules by hypergeometric tests and performed the correction for multiple testing by the Bonferroni method (Rupert, 2012). Gene Ontology analysis was performed using the R package of clusterProfiler (Yu et al., 2012), with background genes set at the genes in the analyzed expression matrix. The GO term whose FDR-adjusted p-value < 0.1 and the number of genes in the term is not less than ten was reported.

We proposed integrative regularized non-negative matrix factorization (iRNMF) to learn homogeneous and heterogeneous gene expression patterns across multiple datasets. Single-nucleus gene expression datasets of multiple cell types of human MTG (Hodge et al., 2019) were analyzed using iRNMF, involving glutamatergic neuron (Gluta), GABAergic interneuron (GABA), astrocyte (Ast), oligodendrocyte (Oli), and oligodendrocyte precursor cell (OPC) denoted by X 1 , … ,   X 5 , X i ∈ R n i × g ,   i = 1,2 , … , 5 , where g denotes the number of genes and n i denotes the number of samples in the ith cell type. iRNMF decomposed each gene expression dataset, X i , into three low-dimensional matrices, including the representation of samples in the low-dimensional space H i ∈ R n i × m and the cell type-specific and shared gene module matrices V i ∈ R m × g and W ∈ R m × g , respectively, where m denotes the number of factors/gene modules. As we study the shared and cell type-specific gene expression patterns across cells, we mainly focus on V i and W . Based on V i and W , for each gene module, we first calculated the z-scores of genes and determined the module genes as those with z-score > 1. The modules with no less than 20 module genes were reported. The gene modules determined from W were regarded as shared gene modules, and those determined from V i were regarded as cell-type-specific gene modules. Then, we identified the gene modules significantly enriched with SFARI ASD candidate genes using hypergeometric tests. The gene modules whose Bonferroni-adjusted hypergeometric test p-values < 0.1 were reported as ASD-associated gene modules. By analyzing the shared and cell-type-specific ASD-associated gene modules, we study the shared and cell-type-specific dysregulated gene expression across different cells in ASD.

To show the effectiveness of iRNMF, we compared iRNMF with LIGER (which only imposes regularization on ∑ i = 1 k | | H i V i | | F 2 but not on ∑ i = 1 k | | H i W | | F 2 ). First, we needed to determine the parameter values for LIGER and iRNMF. KL divergence was used to determine the optimal number of gene modules m, and the alignment metric (Butler et al., 2018) was used to determine the regularization parameters. For LIGER and iRNMF, we plotted the median of KL divergence across samples as a function of m to select the saturation point of the curve and also considered mean squared error (MSE), ∑ i = 1 k 1 / ( ni × g ) | | X i − H i ( W + V i ) | | F 2 (Supplementary Figure S1A). Thus, m was set to 100 for both LIGER and iRNMF. For LIGER, the regularization parameter λ was set to 1, which makes the alignment metric reach the minimum value (Supplementary Figure S1B). For iRNMF, we plotted the alignment metric as a function of a combination of λ   1 and λ 2 . The parameter values λ   1 = 0.01 and λ 2 = 10 can make the alignment metric reach the minimum value, while we noticed that λ   1 = 1 and λ 2 = 10 can give the second smallest alignment metric (Supplementary Figure S1C). The regularization parameter of LIGER was determined as 1, which is actually our λ   1 ; to better compare iRNMF with LIGER and analyze the effectiveness of the newly added constraint, we chose λ   1 = 1 and λ 2 = 10 instead. The reason why we chose the minimum alignment instead of the maximum alignment as a criterion is that the input datasets are from different cell types and they could not be aligned together. The alignment metric measures the uniformity of mixing for multiple samples in the aligned latent space, which should be high when datasets share underlying cell types and low when datasets do not share cognate populations (Butler et al., 2018; Welch et al., 2019). As our analyzed datasets are from different cell types, we used the minimum alignment to determine the regularization parameters.

Next, we compared iRNMF with LIGER based on cell representation. For each cell type, we calculated the Pearson correlation between flatten X i and X i ^ = H i ( W + V i ) (Supplementary Figure S2A) to evaluate the reconstruction. Also, we calculated sample–sample distance matrices using X i and X i ^ and then flatten the distance matrices to calculate their Pearson correlation (Supplementary Figure S2B). Both correlations of iRNMF are slightly better than those of LIGER. Then, for each cell, we calculated the Pearson correlation between the gene expression levels of this cell in X i and X i ^ (Supplementary Figure S2C). We found that iRNMF is better than LIGER in the cell type GABA, and in the other four cell types, iRNMF and LIGER are evenly matched. To compare the low-dimensional H i obtained from LIGER and iRNMF, we combined all H i as H and performed cell clustering based on H to check if different cell types are distinguishable. We performed K-means based on H and calculated the clustering indexes, including ARI (adjusted Rand index), FMI (Fowlkes and Mallows index), JC (Jaccard coefficient), NMI (normalized mutual information), PUR (purity), and SC (silhouette coefficient) (Supplementary Figure S2D). It can be noted that the clustering performances of iRNMF are better than those of LIGER when being faced with datasets of different cell types.

Lastly, we compared the gene modules obtained using LIGER and iRNMF. We calculated gene–gene correlation matrices using X i and X i ^ and then flatten the correlation matrices to calculate their Pearson correlation (Supplementary Figure S2E). We also calculated gene–gene correlation matrices using X i and low-dimensional gene representation W + V i and then calculated their Pearson correlation (Supplementary Figure S2F). It can be seen that both correlations of iRNMF are better than those of LIGER. Then, for each gene, we calculated the Pearson correlation between the expression levels of this gene in X i and X i ^ (Supplementary Figure S2G). The correlations of iRNMF are significantly higher than those of LIGER in four cell types, except for Gluta, in which iRNMF and LIGER are evenly matched. Moreover, we expected that different modules should represent distinct biological functions and should not overlap too much. To evaluate the distinct biological functions of gene modules, we adopted the evaluation way as in Cleary et al (2017), using the number of uniquely enriched gene sets. For each gene module of W and V i , we tested it for enrichment in GO terms and considered its top five significant GO terms. Then, we identified the uniquely enriched GO terms of each module, which are the terms enriched in at most one module of this considered cell type, and calculated the average number of unique gene sets per module (Supplementary Figure S2H). It can be seen that in all cell types, the number of uniquely enriched GO terms of iRNMF is larger than that of LIGER. The comparisons indicate that iRNMF is effective, and the obtained low-dimensional matrices are helpful for subsequent gene module analysis.

Among all shared gene modules determined from W, 46 are significantly enriched with ASD genes (Supplementary Table S1). For the top ten shared ASD-associated gene modules, we list their Bonferroni-adjusted p-values, top three z-score genes (Figure 1A1), and top one significant GO term (Figure 1A2). Some top genes are ASD genes, including PDE1C and MKX in W_M42, GPC6 in W_M10, and NTNG1 in W_M26. The top one significant GO term is all related to synapses, whose dysregulation has been proven to be associated with ASD. Then, we checked which kinds of GO terms are the most common among all GO terms of all shared ASD-associated gene modules and found that the top ten common GO terms are also associated with the functioning of synapses, appearing in all shared ASD-associated gene modules (Figure 1B1). Next, we focused on the modules which have module-specific gene functions, by removing the repeated GO terms between gene modules. There are 36 shared ASD-associated gene modules with module-specific GO terms (Supplementary Table S2). The top ten modules with module-specific gene functions are also the ones shown in Figure 1A1, and their top one module-specific GO term is shown in Figure 1B2. The top three modules most significantly enriched with ASD genes, W_M78, W_M91, and W_M42, are related to cortical actin cytoskeleton, heparan sulfate proteoglycan metabolic process, and regulation of sodium ion transmembrane transport, respectively. The actin cytoskeleton has been associated with ASD and provides a strategy for ASD treatment by targeting actin regulators (Duffney et al., 2015; Hlushchenko et al., 2018). The lacking of heparin sulfate, a proteoglycan involved in a variety of neurodevelopmental processes, has been correlated with ASD (Irie et al., 2012; Pérez et al., 2016). Ion channels, including sodium, calcium, and potassium, are implicated in the etiology of ASD (Daghsni et al., 2018). It can be seen that the identified gene modules are meaningful.

Among all cell-type-specific gene modules, we identified 11, 25, 29, 45, and 14 cell-type-specific ASD-associated gene modules for Ast, GABA, Gluta, Oli, and OPC, respectively (Supplementary Table S1). We list the top ten significant gene modules along with their Bonferroni-adjusted p-values, top three z-score genes, and top one significant GO term (Figure 2). Noted that for the two kinds of neurons, GABA and Gluta, the cell-type-specific ASD-associated gene modules are more significantly enriched with ASD genes and more top three genes are ASD genes, compared with glial cells. Many of the top GO terms of cell-type-specific ASD-associated gene modules are related to synapses, while different gene functions may still be dysregulated in different cell types. For instance, gamma-tubulin complex, cadherin binding, and protein tyrosine kinase activity are associated with Ast-specific ASD-associated gene modules; regulation of microtubule cytoskeleton organization, phosphatase binding, and desmosome are significant in OPC-specific ASD-associated gene modules. These may indicate that different gene functions may be dysregulated by ASD in different cells, demonstrating the cell heterogeneity of ASD.

Then, we checked which kinds of GO terms are the most common among all cell-type-specific ASD-associated gene modules in each cell type. Indeed, the functioning of synapses is important across all cell types (Figure 3A). Next, we focused on the modules which have module-specific gene functions. There are 7, 23, 24, 24, and 8 cell-type-specific ASD-associated gene modules left in Ast, GABA, Gluta, Oli, and OPC, respectively (Supplementary Table S2). We reported the top ten, along with their top three z-score genes (Figure 3B), and top one GO term (Figure 3C). In Ast, locomotory behavior, integral component of the postsynaptic membrane, and cadherin binding are functions specific to the top three modules, Ast_M99, Ast_M63, and Ast_M39. For GABA, it can be noted that inhibitory extracellular ligand-gated ion channel activity is specific to GABA_M7. On the contrary, the modulation of excitatory postsynaptic potential and ionotropic glutamate receptor signaling pathway are specific to Gluta_M10 and Gluta_M99, respectively. These gene functions are obviously associated with particular cell types. Neurons communicate with one another at synapses using two types of signals, electrical and chemical signals. At an electrical synapse, ions flow directly between cells. At a chemical synapse, neurotransmitters pass messages from the presynaptic to the postsynaptic neuron. The major excitatory and inhibitory neurotransmitters in brains are glutamate and GABA (gamma-aminobutyric acid), respectively. For Oli, regulation of dendrite morphogenesis and regulation of gliogenesis are specific to Oli_M82 and Oli_M97. For OPC, protein homooligomerization and endoplasmic reticulum unfolded protein response are specific to the top two modules, OPC_M34 and OPC_M54, respectively. The analysis of module-specific gene functions and top genes of cell-type-specific ASD-related gene modules can facilitate the development of more targeted biomarkers and treatments for ASD.

Next, we further examined the modules with both cell-type-specific and module-specific gene functions, which are those GO terms that only appear in one module of one cell type. In GABA, Gluta, Oli, and OPC, there are 14, 18, 1, and 1 cell-type-specific ASD-associated gene modules that have both cell-type-specific and module-specific gene functions (Supplementary Table S3). It can be noted that more modules have cell type-specific and module-specific gene functions in neuronal cells, emphasizing the neurons are mainly affected by ASD. For the cell types with more than one cell-type-specific ASD-associated gene modules, we show the top ten modules along with their Bonferroni-adjusted p-values, the enriched top one GO term, and the top three z-score genes (Figure 4). Among the top three genes, CTNNA2 in GABA_M42, ZBTB2 in GABA_M40, and SLC9A9 in GABA_M75 are ASD genes. MKX and PDE1C in Gluta_M42, GPC6 and CUX2 in Gluta_M10, and PXDN and NFIA in Gluta_M99 are ASD genes. These gene modules may need more attention. We note that different kinds of gene functions are specific to ASD-associated modules of different cell types. GABA-specific ASD-associated gene modules are responsible for inhibitory extracellular ligand-gated ion channel activity and forebrain neuron differentiation, and so on. Gluta-specific ASD-associated gene modules are responsible for nerve growth factors, excitatory postsynaptic potential, and ionotropic glutamate receptor signaling pathway, and so on. Oli_M60 and OPC_M12 have a cell-type-specific and module-specific function, regulation of bone mineralization and lipid transporter activity, respectively (Supplementary Table S3). These results indicate that in different cell types, different kinds of gene functions may be specifically dysregulated in ASD, highlighting the cell heterogeneity of ASD.