Scoliosis-related variants and their annotations were downloaded from ClinVar database (https://www.ncbi.nlm.nih.gov/clinvar/); Data was cleaned to remove those variant samples with variant in intron, and was parsed for annotations, such as gene name, gene ID, variant type and others for each sample. The full coding DNA sequence (CDS) for each variant was generated based on the CDS of gene ID and its variant annotation, and variant CDS traits of dinuleotide (DNT), DCR, codon usage, codonpair and amino acid (AA) were decomposed with a reported decomposer (Li et al., 2022), respectively producing vectors with dimension of 48 (DNT), 1,536 (DCR), 64 (codon usage), 3,721 (codonpair) and 20 (AA) for each sample. The algorithm for counting DNT (Li et al., 2020), DCR (Li et al., 2022), codon, codonpair and AA (Jiang S. et al., 2023) were designed according to Formula 1‐5 respectively. Statistical description of variants was performed based on sample annotation information. f r e q x n y m = ∑ x n y m ∑ i = 1 16 x n y m , ( x , y = T , C , A   o r   G , m = n + 1   f o r   m ≤ 3 , m = n − 2   f o r   m = 4 , n = c o d o n   n t   p o s i t i o n   1 , 2 , o r   3 (1) f r e q w i x j y k z l = ∑ w i x j y k z l ∑ i = 1 256 w i x j y k z l , ( w , x , y , z = T , C , A   o r   G , j   k = i + 1   i + 2 , j     k = j   k   i f   j     k ≤ 3 , e l s e   j   k = j − 3     k − 3 , k = k + 1 , l = l   i f   l ≤ 3   e l s e   l = 3 , i = c o d o n   n t   p o s i t i o n   1 , 2 , o r   3 (2) F r e q C o d o n = c o u n t C o d o n × 64 × 3 C D S   l e n g t h , c o d o n = e a c h   o f   t h e   64   t y p e s   o f   c o d o n s (3) F r e q c o d o n p a i r = c o u n t c o d o n   p a i r × 3721 × 3 C D S   l e n g t h , c o d o n   p a i r = e a c h   p a i r   o f   t h e   64   t y p e s   o f   c o d o n s (4) F r e q A A = c o u n t A A × 20 × 3 C D S   l e n g t h , A A = e a c h   o f   t h e   20   t y p e s   o f   a m i n o   a c i d s (5)

To learn the association of genome information with the scoliosis phenotype, in the gene of FBN1, LAMA2 and SPG11, the clustering and separation of the variant samples were analyzed based on the DCR and other features of these genes. The composition feature vector of DNT, DCR, codon usage, codonpair and AA respectively with a dimension of 48, 1,536, 64, 3,721 and 20 was reduced to two main components with Uniform Manifold Approximation and Projection (UMAP) and were scattered with sample label of “Pathogenic”, “Pathogeniclikely”, “Benign”, “Benignlikely” and “Uncertain”. Hierarchical clustering of these samples was also performed with a python package of sns. clustermap, based on the Euclidean distance of the above-mentioned five types of features with each sample labelled. The components reduced from the compositional features were normalized with the following formula 6. X n o r m a l i z e d = X − X min / X max − X min (6)

A CNN classifier was designed to predict scoliosis genotypes of FBN1 based on genomic compositional features with the variants labelled as “Pathogenic” as positive samples and with the variants labelled as “Benign” or “Benignlikely” as negative samples. A random downsampling was performed to guarantee a sample balance between the two types of data. DCR with 1,536 dimensions were selected to train the classifier with a network structure of CNN. DCR data was then randomly split into training dataset and validation dataset with a 5-fold cross-validation method, then was reshaped into an array with the size of (6, 16, 16) and finally was input into the three-layer 3D-CNN model. The CNN models were set with a convolution kernel with size of (1, 3, 3), with a pooling layer of (1, 2, 2) via average pooling, and with a padding layer of (0, 1, 1), and with a stride of (1, 1, 1). The batch size, learning rate, and training epochs were optimized respectively. 768-dimensioned output from two rounds of convolution of the 1536-dimensioned DCR was linearly transformed for two times, firstly into 192-dimensioned and secondly into 2-dimensioned output, which was finally calculated with Softmax function (formula 7) to output the probability for each of the two scoliosis risks (positive and negative). Detailed parameters are epoch_num = 100, split_size = 0.2 (for training and validation dataset), lr = 0.005 and batch_size = 20. The classification performance of the models was evaluated by with receiver operating characteristic curve (ROC) and the area under curve (AUC), confusion matrix. S o f t m a x : f x i = e x i / ∑ j = 1 J e x j (7)

The variants labelled with “Pathogeniclikely” or “Uncertain”, from ClinVar database were assessed for their scoliosis risk with the trained CNN classifier based on their DCR features. The DCR of FBN1 CDS was reshaped into a dimension of (6, 16, 16), and then were transformed into a tensor, then was input into the loaded model of “DCR-based 3D-CNN for scoliosis. txt”. The prediction of 1 for high scoliosis risk and 0 for low scoliosis risk and the probabilities for the two risk results were finally output. The scoliosis risk was further analyzed in more details by statistically describing distribution of these samples on various variant labels.

To landscape the distribution of scoliosis variants on the 3D structure of FBN1, a reference FBN1 protein (NM_000138.5, NP_000129.3) with truncated N-terminal of 1,100 amino acids, where most of the amino acid variants were located, was utilized to predict its 3D structure. The structure prediction was performed with AlphaFold2 of offline version (Jumper et al., 2021). A virtual environment was first activated with the command “conda activate alphafold” under the same path as the FBN1 fasta file. Then the prediction was performed with the command “python/docker/run_ docker. py--fasta paths = Reference_NM_000138.5_FBN1. fasta--max_template_data = 2020-05–14”. The ranked0_.pdb was taken as the most optimized result and was visualized with PyMOL (version 2.5.7). All variants amino acids were manually labelled in red for the variants in the Pathogenic group, in purple for predicted risk variants in the Pathogeniclike group, and in violet for risk variants in the Uncertain group.

Significance was evaluated with paired t-test for the Principal Component Analysis (PCA)-reduced PCA1 or PCA2 value of the full-connected layer between Pathogenic and Benign groups, and for the probability values of label 1 and label 0 either for negative or positive variants. GraphPad Prism (version 9.0.0) was utilized for statistical analysis and figure plotting. A p-value threshold of 0.05 was taken as statistical significance.

The workflow of this study was set to parse, decompose and predict the genotype with high scoliosis risk, with six successive models. Firstly, full data of variant samples associating with scoliosis or not were downloaded from NCBI and were parsed for their genotypes and annotations (Figure 1A). The distribution of these data on various types of annotation labels was statistically analyzed and the full CDS were generated based on their variant annotations (Figure 1B). Secondly, the genomic compositional traits, such as dinucleotides (DNTs), DCR, and others, for each scoliosis-related gene were calculated based on our previously reported algorithm (Figure 1C), and then were analyzed with unsupervised machine learning methods (Figure 1D). Finally, a Convolutional Neural Network (CNN) classifier based on DCR features of FBN1 was trained with two labels: high and low risk for scoliosis (Figure 1E). It was used to predict the scoliosis risk of variants with unclear risk and was interpreted by optimizing significant amino acids on the 3D structure of FBN1 (Figure 1F).

The distribution of scoliosis-related variants was analyzed. Most of these variants whether with high scoliosis risk (Pathogenic or Pathogeniclikely), with low scoliosis risk (Benign or Benignlikely) or with unknown risk (Uncertain) are located within exon/cDNA, rather than within intron (almost similar value for Total and cDNA, Figure 2A). Most of the variants in cDNA were not synonymous, causing variants in protein level for scoliosis-pathogenic samples (Synonymous/Protein = 0.16 or 0.25 for the samples labelled with pathogenic or pathogeniclikely cDNA, Figure 2A). Counting of scoliosis-related genes indicated that NF1, FBN1, LAMA2 and SPG11 led the top list of genes concerned with scoliosis (Figure 2B). And single nucleotide polymorphism (SNP) dominated the variant type list for most scoliosis-related genes (Figure 2C). Most variant types for these SNPs were the base transition of C > T, G > A and G > T for pathogenic variants (Figure 2D).

To overview the distribution of scoliosis-related variants based on genome decomposition traits, decomposed genome data of DCR and other traits was plotted with scoliosis risk labelled (high scoliosis risk: Pathogenic or Pathogeniclikely), low scoliosis risk: Benign or Benignlikely), post dimension reduction with Uniform Manifold Approximation and Projection (UMAP). A marked clustering of the samples with the same risk label and a long-distance separation of the samples with different labels were observed for DCR features of FBN1 for the two main components of UMAP1 and UMAP2 (Figure 3A). Whereas these samples with the four labels were mixed in the distribution of UMAP1 and UMAP2, either for the trait of DNT (Figure 3B), codon (Figure 3C), codonpair (Figure 3D) or AA (Figure 3E) of FBN1. The intra-risk type clustering and inter-risk type separation were repeatedly observed on the DCR features of LAMA2 and SPG11 (Figures 3F, G), rather than on the AA features of the two genes (Figures 3H, I). The relatively hierarchical clustering (high scoliosis risk of Pathogenic vs. Benign and Benignlikely) was also indicated by DCR of FBN1 (Figure 3J). However, such clustering was not significantly binary on the raw DCR data of FBN1, implying an incapability of the raw DCR data for a binary classification for the two scoliosis risk types and a need of feature optimization of DCR by deep learning.

In light of the high association between high-dimensional DCR features and scoliosis risk, however their nonlinear distribution depending on scoliosis type, a deep learning predictor based on the 1,536-dimension DCR trait was trained for scoliosis risk classification for the leading gene of scoliosis risk, FBN1. A CNN network (Li et al., 2023; Li et al., 2022) was utilized for the binary classification based on randomly sampled FBN1 data. Firstly, a binary classifier of Convolutional Neural Networks (CNN) with labels of high risk (Pathogenic) and low risk (Benign or Benignlikely) was built based on DCR of FBN1. The 1536-dimension DCR were sequentially subject to two rounds of convolution, two times of linear transformation, and one-time Softmax transformation, and then outputted two classification labels of 1 (high) and 0 (low). The CNN classifier was quickly converged with optimized parameters, indicating an early decrease and following long micro-movement of training loss (Figure 4A). The fully connected layer of 768 dimensions was reduced with Principal Component Analysis (PCA) to visualize the concentration of key information of DCR by deep learning. It was indicated that the key difference between high and low risk data concentrated on the first component, showing a significant difference of PCA1 (p < 0.001, Figure 4B), whereas the insignificant difference of PCA2 (Figure 4C), between the two groups. Then the prediction performance was evaluated with an independent validation dataset with a confusion matrix and Receiver Operating Characteristic Curve_Area under Curve (ROC_AUC). A right angle-like ROC and an AUC value of more than 0.92 (Figure 4D) and a confusion matrix with 100% of accuracy for low risk data and 92.72% of accuracy for high risk data were obtained on independent training data (Figure 4E). Thus, DCR of FBN1 was classifiable and predictable for high or low scoliosis risk, post convolutional transformation.

The trained CNN classifier was utilized to predict the scoliosis risk of the variant labelled with “Pathogeniclikely” or “Uncertain” from the ClinVar database. 179 variants (Supplementary Table S1) were predicted as high risk (label 1) and the other 2044 variants were low risk (label 0) based on their genomic DCR (indicated as positive and negative respectively, Figure 5A), with significant probability difference (p < 0.001 between label 0 and 1, respectively for both positive and negative groups, Figure 5B). There were 118 variants with false stop in FBN1 CDS and 61 variants without false stop (Figure 5C). Deletion and duplication were the main variant types for the scoliosis risk variants (Figure 5D). In more detail, the frequency of various SNP types for all scoliosis risk variants and for predicted scoliosis risk variants were compared.

Interestingly, there were relatively high levels of SNP types of C > T, C > A, G > T and T > A were observed in either the Pathogenic group from the ClinVar database (Figure 5E) or the predicted scoliosis risk group (Figure 5F). Moreover, the frequency vector for these SNP types were highly similar between the two groups, indicating a linear correlation with high slope signification (Figure 5G). Such linear correlation was also observed for low scoliosis risk prediction in the frequency vector for these SNP types between the two groups (Figure 5H). Detailed percentages for high (label 1) and low (label 0) scoliosis risk variants for all five groups (Uncertain, Pathogeniclikely, Pathogenic, Benignlikely and Benign) were plotted (Figure 5I). Therefore, the DCR embedding of FBN1 and the trained CNN classifier worked well in predicting scoliosis-risk variants, based on the interpretability of the high similarity in variant between pathogenic variants and predicted variants.

We further interpreted the dependence of the predictor on protein sequence or structure. Firstly, the frequency of variants for all variants in each of the four groups was analyzed. It indicated a marked biased variant distribution on the FBN1 sequence, with most of the pathogenic variants located in the N-terminal, followed by pathogeniclikely variant, with benign and uncertain variants in the C-terminal (Figure 6A). The variant distribution within the 2D structures of FBN1 was also analyzed. The 3D structure of the human FBN1 protein with Alphafold2 and mutated amino acid with the 2D structure of loop, helix and sheet was counted. It was shown that most of these variants were within the loop structure of FBN1, without marked differences among the four groups (Figure 6B). The full view of the 3D structure of full-length FBN1 was visualized in a cartoon (Figure 6C) and surface form (Figure 6D), indicating regular repeated units of beta-sheet and loop/alpha-helix. A truncated FBN1 with the 1,100 amino acids in the N-terminal indicated the repeated sheet-loop/helix units in more detail (Figure 6E). All the amino acid variants in FBN1 for the Pathogenic (Figure 6F), the Pathogeniclikely (Figure 6G) and the Uncertain group (Figure 6H) distributed uniformly, without any dependence on 2D structure. Taken together, the DCR-based CNN predictor for scoliosis-risk prediction was not dependent on the protein sequence or protein structure of FBN1.