Labeled variant call dataset from previous work was used to test the MultiCapsNet model (Ainscough et al., 2018). This dataset contains more than 41,000 samples, which are assembled to train models for automating somatic variant refinement. Each sample in the dataset is manually labeled as one of four tags by the reviewer: “somatic”, “ambiguous’, “germline”, and “fail”, which represent the confidence of a variant call by upstream somatic variant caller. As in previous work, we merged the variant calls labeled as “germline” and “fail” into a class named “fail”. The number of instances in each class are around 10,000, 13,000, 18,000 for “ambiguous”, “fail”, and “somatic”. There are 71 features that are associated with each sample, including cancer types, reviewers, tumor read depth, normal read depth, and so on. According to the data sources and data attributes, we divided these 71 features into eight groups (Supplementary Table S1). Group 1 contains nine cancer types, and is encoded as one-hot encoding vector. We call group 1 as “Disease” because it indicates the disease to which each variant call belongs. Group 2 contains four reviewers, and is encoded as one-hot encoding vector. We call group 2 as “Reviewer”. Group 3 contains information of “normal VAF”, “normal depth”, “normal other bases count”, and is called as “Normal_pro”, short for “Normal properties”. Group 4 contains 13 features that describe reference reads in normal, including base quality, mapping quality, numbers of mismatches, numbers of minus and plus strand, and so on. We call group 4 as “Normal_ref”. Group 5 contains 13 features extracted from variant reads in normal, also including base quality, mapping quality, numbers of mismatches, numbers of minus and plus strand, and so on. We call group 5 as “Normal_var”. The last three groups contain features drawn from tumor instead of normal in previous three groups. As same as Group 3, 4, and 5, we label group 6, 7, and 8 as “Tumor_pro”, “Tumor_ref”, and “Tumor_var” respectively.

The mouse scRNA-seq is measured by Microwell-Seq (Han et al., 2018). We downloaded scRNA-seq data and the annotation information through the link provided by the authors (https://figshare.com/s/865e694ad06d5857db4b). Then we use the annotation information to select parts of data from whole dataset. The cell types we chose include “Cartilage cell”, “Secretory alveoli cell”, “ Epithelial cell ”, “Kupffer cell”, “Muscle cell”, “Dendritic cell”, “ Spermatocyte”, and the number of instances in each cell type are 527, 1,195, 1,219, 356, 626, 717, 353. Moreover, we only use the genes contained in prior knowledge (Lin et al., 2017) to fit the model structure, and set the default value to zero when the downloaded scRNA-seq data does not contain this gene (Han et al., 2018).

A SCENIC example dataset was used to compare the performances of MultiCapsNet and SCENIC (https://scenic.aertslab.org/examples/). The dataset (sceMouseBrain.RData) contains seven cell types of mouse cortex and hippocampus (Zeisel et al., 2015) [“astrocytes_ependymal” (224), “endothelial_mural” (235), “interneurons” (290), “microglia” (98), “oligodendrocytes” (820), “pyramidal_CA1” (939), and “pyramidal_SS” (399)].

In the architecture of our multiCapsNet model, there are l neural networks corresponding to l input modular data. u i =   tanh ( W p i x i )             i ∈ [ 1,2 … , l ] (1) x _i represents i’s input modular data. W p i represents weight matrices of neural networks with dimension (n, r _i ), where the r _i is the length of the input modular data x _i . The output u i of each neural network i   ( i ∈ [ 1,2 … , l ] ) is a vector with length n, viewed as “primary capsule” in the model. The inputs standardization part converts the modular data with different type and length into real valued vectors with equal length n (n = 8 by default).

The standardized information is subsequently delivered through primary capsule to the capsule in the final layer by “dynamic routing” (Supplementary Figure S1). Each capsule in the final layer, named “type capsule”, corresponds to each cell type. They are denoted as vectors v j , where j ∈ [ 1,2 … , k ] , k is the number of cell types and m is the length of vectors. The capsule network module is implemented in Keras (https://github.com/bojone/Capsule).

Prior to the “dynamic routing” process, the primary capsules are multiplied by weight matrices W i j to produce “prediction vectors” u ^ j | i . u ^ j | i   = W i j u i (2)

Then the iterative dynamic routing begins. Firstly, the “coupling coefficients” c i j is calculated by formula: c i j   = exp ( b i j ) ∑ k ⁡ exp ( b i k ) (3) Where b i j is an intermediate parameter with initial value of zero, representing the inner product of the prediction vector and type capsule vector.

In order to compute the b i j for next round iteration, the weighted sum s j over all k prediction vectors u ^ j | i is calculated by formula: s j   =   n o r m a l i z e ( ∑ i c i j u ^ j | i ) (4)

Secondly b i j is computed by the dot product of u ^ j | i and s j as the last step of one round dynamic routing process. b i j   = u ^ j | i . s j (5)

After several rounds of dynamic routing, the type capsule v j is calculated by a non-linear “squashing” function: v j   = ‖ s j ‖ 2 0.5 + ‖ s j ‖ 2 s j ‖ s j ‖ (6)

The following pseudocode illustrates the implementation of MultiCapsNet.

1) f o r   a l l   p r i m a r y   c a p s u l e   i : u i =   A c t i v a t i o n   F u n c t i o n ( W p i x i )

The implementation of MultiCapsNet can be found in https://github.com/wanglf19/MultiCapsNet.

In the somatic variant refinement task, the eight groups mentioned above in the section of “Datasets and data preprocessing” correspond to eight input sources. Therefore, there are eight neural networks corresponding to eight groups of input modular data (l = 8). u i =   tanh ( W p i x i )           i ∈ [ 1,2 … , 8 ] (7)

After the input standardization part, the input data x _i is converted into a primary capsule u _i having the same length. Next, the standardized information stored in the primary capsules would be delivered to the final layer capsules by “dynamic routing”. The capsules in the final layer, which corresponds to labels of variant calls, is called “label capsule”. In capsnet, the non-linear “squashing” function ensure that short vectors get shrunk to almost zero length and long vectors get shrunk to a length slightly below 1 (Sabour et al., 2017). The length of the label capsule represents the probability that a variant call is either “ambiguous”,“fail”, or “somatic” (Figure 2). To evaluate the performance of the model, we use the “area under the curve” (AUC) score as previous (Ainscough et al., 2018) and prediction accuracy.

The MultiCapsNet could integrate prior knowledge into its structure. In brief, PPI information store in BIOGRID (Stark et al., 2006) and HPRD (Keshava Prasad et al., 2009), and PDI coming from DREM 2.0 (Schulz et al., 2012), are used as prior knowledge for specifying network connections between the inputs and the primary capsules (Figure 4A), just as previous work used this prior knowledge to specify network connections between the inputs and neurons (Lin et al., 2017). For example, the prior knowledge indicates that Gene₁,…, Gene_n are regulated by a TF (colored with green), so there are connections between Gene₁,…, Gene_n and primary capsule representing corresponding TF (green connection); the prior knowledge indicate that Gene₂,…, Gene_n are regulated by a TF (colored with blue), then there are connections between Gene₂,…, Gene_n and primary capsule representing corresponding TF (blue connection); and the prior knowledge indicates that Gene₂, Gene₃,…, are in a subnetwork of PPI network (colored with red), then there are connections between Gene₂, Gene₃,…, and primary capsule representing corresponding PPI subnetwork (red connection). Although there is only one input source, namely scRNA-seq data, the input source can be decomposed into several parts by integrating prior knowledge, and each part is connected to a primary capsule. Therefore, we also took a single input source integrated with prior knowledge as an input from multiple sources, each of which is associated with a TF or a PPI subnetwork (Figure 4B).

In total there are 696 input modular data, with 348 TF-targets relationships extracted from PDI information and 348 PPI subnetworks. Therefore, there are 696 neural networks corresponding to 696 modular data (l = 696). u i = tanh ( W p i x i )           i ∈ [ 1,2 … , 696 ] (8)

After the input standardization part, the input data x _i is converted into primary capsule u _i with the same length. Next, the standardized information stored in the primary capsules would be delivered to the final layer capsules by “dynamic routing”. The capsules in the final layer, which correspond to cell types, is called “type capsule”.

The SCENIC is a workflow for simultaneous reconstruction of gene regulatory networks and identification of cell states using scRNA-seq data (Aibar et al., 2017). The workflow consists of three modules (R/bioconductor packages): GENIE3 (GRNboost), RcisTarget, AUCell. The first two modules were responsible to find potential TF-targets relationships based on co-expression and subsequently select the highly confident TF-target regulation according to TF-motif enrichment analysis. After that, several potential TF-target relationships across all cell types, called regulons, were identified in the dataset. The AUCell would score the activity of these regulons in each single cell. Finally, the unsupervised method is used to cluster cell, identify cell types and states based on the scores of the regulongs, which are used as features for each cell. In our model, we utilized the regulon information identified by the first two modules of SCENIC as the prior knowledge to specify the connections between input and primary capsules (Supplementary Figure S2A). The dataset, intermediate results and the output of SCENIC for a mouse brain example were downloaded from the website (https://scenic.aertslab.org/examples/). The regulon information was extracted from the intermediate result file (regulons_asGeneSet.Rds).

In total there are 253 regulons, which specify TFs and their target genes. Therefore, there are 253 neural networks corresponding to 253 modular data (l = 253). u i = tanh ( W p i x i )           i ∈ [ 1,2 … , 253 ] (9)

After the input standardization part, the input data x _i is converted into primary capsule u _i with same length. Next, the standardized information stored in the primary capsules would be delivered to the final layer capsules by “dynamic routing”. The capsules in the final layer, which correspond to cell types, is called “type capsule”.

In scCapsNet, we showed that the average coupling coefficients represent the contribution of the primary capsules to the final layer type capsules for each cell type (Wang et al., 2020). Similarly, in the multiCapsNet model, the type (label) capsule v j derives from a weighted sum of prediction vectors u ^ j | i . The weights are the coupling coefficients c i j and the magnitude of these coefficients could roughly be regarded as the contribution of the primary capsules u i to the type capsules v j . Each sample (single cell, somatic variant) generates its own coupling coefficients. The average coupling coefficients for samples with same type (label) are calculated by the formular: c i j t y p e   a v e r a g e   =   ∑ t y p e c i j t y p e ∑ t y p e 1   (10)

Therefore, each classification category (cell type/variant call label) corresponds to an average coupling coefficients matrix ( c i j t y p e   a v e r a g e ), called type average coupling coefficients, with rows representing type capsules and columns representing primary capsules. The type average coupling coefficients matrix could be plotted as heatmap for visualization of data. For each classification category (cell type/variant call label), the corresponding type average coupling coefficients matrix contain an effective type capsule row, which is the row whose type is consistent with this classification category. For example, the effective type capsule row in the type average coupling coefficients matrix ( c i j t y p e   2     a v e r a g e ) is the row c i 2 t y p e   2     a v e r a g e . In this row, the magnitude of each element could be regarded as the importance score of the corresponding primary capsule to this classification category. The effective type capsule rows of all classification categories ( c i 1 t y p e   1     a v e r a g e ,     c i 2 t y p e   2     a v e r a g e ,     c i 3 t y p e   3     a v e r a g e …) could be organized into a new matrix, visually represented as an overall heatmap.

A neural network with sigmoid activation function was implemented in Keras. The random forest and nearest-neighbour are implemented with the Python package “scikit-learn”. The comparison transformers model was originally used for IMDB movie review sentiment classification dataset. This transformer model contains the embedding layer for embedding the words into vectors and the Multi-head attention layer (https://github.com/bojone/attention/). We replace the embedding layer with our data standardization layer, and retain the Multi-head attention layer for classification.

The variant call dataset (Please refer to Datasets and data preprocessing in METHODS section for the details) was randomly divided into training set and validation set with a ratio of 9:1. Our MultiCapsNet model performs well in the classification of variant call (Figure 4). The results show that the AUC of the MultiCapsNet model is 0.94, 0.99, and 0.97, respectively, in the classification categories of “ambiguous”, “fail”, and “somatic” (Figure 4A). These AUC scores are similar with those obtained by the Multi-head Attention model (0.93, 0.98, 0.96), feed forward neural network (0.93, 0.99, 0.96), and random forest (0.96, 0.99, 0.98) (Ainscough et al., 2018). Meanwhile, the average prediction accuracy of the MultiCapsNet model is around 0.873, similar to those obtained by the Multi-head Attention model (0.866), and slightly lower than that of feed forward neural network (0.887), and random forest (0.895).

In MultiCapsNet, the coupling coefficient c _ij is viewed as important scores, which is the weight that measure the contribution of each primary capsule to the final layer type capsule. Each input would generate its own coupling coefficient, and the type average coupling coefficient is the average over all the inputs with same classification category. After MultiCapsNet model training, the type average coupling coefficients for each variant label (“ambiguous”, “fail”, and “somatic”) were calculated and visualized as heatmaps (Supplementary Figure S3A) (Please refer to “METHODS” section for the detailed calculation formula of type average coupling coefficients). In each type average coupling coefficient, the most important row, named as “effective type capsule row”, is the row whose type is consistent with this classification category. The overall heatmap is assembled with the “effective type capsule row” which describes the importance scores of all the data sources for distinct category classification (Supplementary Figure S3B). Therefore, the overall heatmap also shows the contribution of each data source to the recognition of each variant labels (“ambiguous”, “fail”, and “somatic”). For example, the data source of “Disease” has the contribution to the classification of “somatic” category and the “Reviewer” source contributes to the classification of “ambiguous” category. The “Tumor_var” source is the most important one for the classification of all the three categories (Supplementary Figure S3B). Over 9 repetitions, the values of each row in 9 overall heatmap are averaged to determine the importance scores of each data sources for the classification of all the categories in MultiCapsNet model (Figure 3B). In feed forward neural network model, the feature importance is measured by average change of AUC after randomly shuffling individual features. Based on the step of features grouping, we added the feature importance scores belonging to the same group together, and take these values as importance of data sources (each group) in feed forward neural networks model (Figure 3B). Then, we calculated the correlation between the data source importance scores obtained by our MultiCapsNet model and those provided by feed forward neural network model. Although our MultiCapsNet model is substantially different from the previous feed forward neural network, and the source importance measuring methods are also different, there is very high correlation between them (Pearson Correlation Coefficient = 0.876) (Figure 4B). Both models indicate that tumor variant group is very important for variant call classification.

The dataset is a portion of mouse scRNA-seq data measured by Microwell-Seq, which consists of nearly 5,000 cells of seven types and 9,437 genes (Please refer to METHODS section for the details). The MultiCapsNet model that integrates prior knowledge (Figure 4) was trained and tested by using this dataset. The average validation accuracy and F1 score are around 97%, comparable with those generated by the feed forward neural network, Multi-head Attention model and random forest (Supplementary Figures S4A, B). After training, the average coupling coefficients, which represent the contribution of the primary capsules (TF/PPI) to the type capsules (Cell type), were calculated and visualized as heatmaps for each cell type (Figure 5A). In each heatmaps, we should clearly observe that the high value elements in the average coupling coefficients (dark line in the plot) are exclusively located in the effective type capsule row. Then, the corresponding type capsule row was selected from each heat map in Figure 5A, and organized into an overall heatmap (Figure 5B).

We repeat the training process 9 times and generate nine overall heatmaps accordingly. Based on the average value of the nine overall heatmaps, the top 10 relevant TFs/PPI subnetwork was generated (Supplementary Table S2). Most of the top 10 relevant TFs/PPI subnetwork were specific to one cell type, and many of them have been reported to be associated with corresponding cell types previously (Figure 5C). For example, Gata1 and Gata2 are top contributors for dendritic cell recognition. Previous work indicated that Gata1 regulates dendritic cell development and survival (Gutiérrez et al., 2007), Gata2 regulates dendritic cell differentiation (Onodera et al., 2016). Srf and Yy1 are ranked as the top contributors for muscle cell recognition by the model. However, Srf is required for skeletal muscle growth and maturation (Li et al., 2005), Yy1 is associated with increased smooth muscle specific gene expression (Favot et al., 2005). FoxA2 and FoxA3 are ranked as top contributors for Cartilage cell recognition, and FoxA2 and FoxA3 are necessary to promote high-level expression of several hypertrophic chondrocyte markers (Ionescu et al., 2012). The model reports Rxrg, Rara, Rarg, Rarb, Rxra, and Rxrb as top contributors for Kupffer cell recognition. Previous research report RA receptor (RAR) and retinoid X receptor (RXR) were expressed by Kupffer cells (Ulven et al., 1998; Ohata et al., 2000). Pgr is ranked as a top contributor for secretory alveoli cell recognition, and the progesterone receptor (Pgr) knockout mouse demonstrated that Pg is required for alveolar morphogenesis (Oakes et al., 2006). Topors is ranked as a top contributor for spermatocyte recognition. Previous work indicates dtopors, the Drosophila homolog of the mammalian Topors, plays a structural role in spermatocyte lamina that is critical for multiple aspects of meiotic chromosome transmission (Matsui et al., 2011).

To further demonstrate the effectiveness of our MultiCapsNet model to reveal cell type related TFs from scRNA-seq data, we compare it with established single-cell regulatory network inference methods: SCENIC (Single-cell regulatory network inference and clustering) (Supplementary Figure S2A). The scRNA-seq data from mouse cortex and hippocampus were used to evaluate these two methods (Please refer to METHODS section for the details).

After MultiCapsNet training, the average coupling coefficients in the overall heatmap would indicate the most relevant TFs associated with each cell type (Supplementary Figure S5). We repeated the experiment 9 times, the average validation accuracy was 97%, and the average F1 score was around 95%, which were comparable to the results generated by feed forward neural network, Multi-head Attention model and random forest (Supplementary Figures S4C, D). According to the average value of nine overall heatmaps, the top 30 relevant TFs could be generated (Figure 6A left; Supplementary Table S3 top). The original regulon may contain TFs that label the 253 regulons. In order to eliminate the influence caused by the expression of those labeling TF, the potential TF-target relationships that exclude the labeling TF in the set of target genes are also made (Supplementary Figure S2B). We also repeated the training process of MultiCapsNet that integrated with those new potential TF-target relationships. After training, the top 30 relevant TFs could also be generated according to the average value of the nine overall heatmaps (Figure 6A right; Supplementary Table S3 bottom). The results show that the inclusion or exclusion of labeling TF has little influence on prediction accuracy and interpretability of the model. The overlap rates of top 30 most relevant TF of each cell type (around top 10% of total TFs) between model including labeling TF and that excluding labeling TF are very high, around 90% for every cell type (Figure 6B).

Many high score TFs predicted by MultiCapsNet are consistent with that reported by SCENIC (Aibar et al., 2017). For example, in both methods, Rorb is identified as a relevant TF for astrocytes; Ets1, Elk3, and Gata2 are identified as relevant TFs for endothelial-mural cells; Zmat4, Dlx5, Dlx2, and Dlx1 are identified as relevant TFs for interneurons; Maf, Rel, Cebpa, Cebpb, Nfatc2, Prdm1, Nfkb1, and Stat6 are identified as relevant TFs for microglia; Sox10 and Sox8 are identified as relevant TFs for oligodendrocytes. Besides the TFs listed above, MultiCapsNet also detected several high confidence cell type relevant TFs that are also found by SCENIC. For example, Rfx3 shows a high association with both pyramidal SS and CA1 cells. Previous studies reported that downstream target of Rfx3 displayed cytosolic expression in pyramidal neurons (Remnestål, 2015) and Rfx3 expresses in cortical pyramidal neurons (Benadiba et al., 2012). Neurod2 is also identified as a relevant TF for both pyramidal SS and CA1 cells. Previous studies reported that Neurod2 coordinates synaptic innervation and cell intrinsic properties to control excitability of cortical pyramidal neurons (Chen et al., 2016). Cux1 has been identified as a relevant TF for pyramidal SS cells, and Cux1 has been reported as a restricted molecular marker for the upper layer (II-IV) pyramidal neurons in murine cerebral cortex (Li et al., 2010). smarca4 has been identified as relevant TF for pyramidal CA1 cells, and Brg1/smarca4 deficiency leads to mouse pyramidal neuron degeneration (Deng et al., 2015). Ezh2 has been suggested as a relevant TF for oligodendrocytes, and the expression of Ezh2 in OPCs (oligodendrocytes precursor cells), even up to the stage of pre-myelinating immature oligodendrocytes, remains high (Copray et al., 2009) (Figure 6C). Furthermore, the MultiCapsNet found that Rpp25 is strongly associated with interneurons which SCENIC did not, and Rpp25 has been reported up-regulated in GABAergic interneuron (Fukumoto et al., 2018).