Integrating AI and genomics: predictive CNN models for schizophrenia phenotypes

2.1 Convolutional Neural Network (CNN) model architecture

This study utilizes a Convolutional Neural Network Architecture, a deep learning model tailored for processing data with grid-like topology, commonly employed for image and pattern recognition tasks [27]. CNNs are highly efficient for feature extraction due to their hierarchical structure, low number of trainable parameters compared to fully connected networks, and adaptability in learning abstract data representations [28]. This architecture has been expanded to applications beyond images, including genomic sequence data, which exhibits a regular grid structure, making CNNs well-suited for this task [29].

The proposed model is structured with multiple layers, each playing a crucial role in identifying complex patterns within genomic data. The model architecture comprises two convolutional layers, each seamlessly integrated with a max-pooling layer to reduce dimensionality and minimize overfitting effectively. This configuration enables the model to capture intricate and broader genomic patterns while preserving a compact and efficient design. As shown in Figure 1, the CNN starts with a 2-dimensional (2D) convolutional layer, ideal for handling the structured input data despite the dataset’s one-dimensional nature since 2D convolutions allow for better spatial feature extraction. A 2D convolutional layer was chosen to capture the spatial patterns in the genomic data, even though our dataset is one-dimensional, meaning that the DNA sequences are represented as linear strings. The 2D approach is ideal for this case because, in the 1D one, the model would only be able to scan along the length of the sequence, which would limit the model’s ability to detect the interactions between more distant regions; however, by using a 2D convolution, we can reshape the 1D data by organizing the genomic sequences into smaller segments, allowing the model to learn patterns across both the sequence and different features (segments) simultaneously, ending up learning both local dependencies and higher order spatial relationships [30].

Figure 1:

Convolutional Neural Network Structure: the designed architecture efficiently extracts patterns in genomic data to predict outcomes based on genome variants of each individual.

The addition of more than one convolutional layer allows the network to learn a hierarchical feature representation. Early layers focus on understanding simple patterns, such as edges or basic genomic motifs, while deeper layers capture more complex features, such as interactions between genetic variants [31]. To further mitigate underfitting, additional convolutional layers were incorporated, thus increasing the feature learnability by deepening the network structure [32].

Max-pooling layers, placed after each convolutional layer, are crucial in reducing the spatial dimensions of the feature maps, which aids in controlling overfitting by summarizing regions of the input [33]. These pooling layers enhance the network’s ability to generalize by down-sampling the feature maps, making the network less sensitive to minor variations in the input data, a typical issue in biological datasets [34].

Finally, after the convolution and pooling layers, a fully connected (FC) layer is included to aggregate the extracted features and pass them to the output layer for classification [35]. The FC layer transforms the multidimensional feature maps into a one-dimensional vector, which is subsequently used for classification tasks. The output layer uses a softmax activation function to predict the class probabilities, such as the presence or absence of a genetic condition, based on the input genome variants [36].

Once the architecture is implemented, the model undergoes training using a labeled dataset, and its performance is evaluated using metrics such as accuracy and loss. These metrics are indicative of the model’s predictive capability and efficiency in classifying genomic data [37]. Further refinements, such as hyperparameter tuning and regularization, are applied to optimize the network’s performance during training [38].

2.2 Dataset and preprocessing

In this study, we utilized data from the Sweden Schizophrenia Population-Based Case-Control Exome Sequencing dataset, which is stored in the dbGaP repository [39] under the accession ID code phs000473.v2.p2. The dataset includes information for 12,380 individuals, consisting of 6,245 controls and 6,135 cases. Of the cases, 4,969 were diagnosed with SCZ, while 1,166 were diagnosed with bipolar disorder. However, for the purposes of this study, only the SCZ cases were analyzed to maintain a binary classification framework. Participants in both the case and control groups were required to be at least 18 years old, with both parents born in Scandinavia [40], [41], [42]. Participants for the study were selected through data from the Swedish Hospital Discharge Register. To qualify, individuals need to have been hospitalized on at least two separate occasions, with each admission confirming a diagnosis of schizophrenia (SCZ). Any individuals diagnosed with other medical or psychiatric disorders that might compromise the precision or reliability of the SCZ diagnosis were not included in the research sample. Control cases were randomly selected from the general population and excluded if they had any recorded history of hospitalization for SCZ. To the best of our knowledge, this dataset represents one of the most considerable whole-exome and whole-genome sequencing resources available for studying schizophrenia, accessed through request.

The initial dataset comprised 1,811,204 variants, all of which had a call Phred-score (QUAL) above 30 [43]. For this study, BCFtools [44], along with its fill-tags plugin, was utilized to remove the bipolar samples from the dataset and to recalculate key metrics for each variant site, including mean read depth (DP), Allele count (AC), Number (AN), and frequency (AF). Further filtering was carried out using VCFTools (Version 0.1.15) [45], applying criteria such as allele count (AC = 0), depth (DP < 8), and a genotype call rate lower than 90 %. Multi-allelic variants were segregated with BCF tools. For newly identified variants, genotypes that were neither the reference nor the alternative allele were marked as missing. A second round of filtering was applied to the split variants, adhering to the same conditions: allele count (AC = 0), depth (DP < 8), and genotype call rate under 90 %. The Variant Quality Score Recalibration process from the Genome Analysis Toolkit Best Practices Workflow (Version 3.8) [46] was followed to guarantee variant call accuracy.

In order to explore potential relationships between genotypic variation and phenotypic expression, chi-squared tests were performed utilizing a 3 × 3 contingency framework, accommodating the three possible genotypic categories across both case and control groups. Variants located on sex chromosomes, as well as insertion-deletion polymorphisms (InDels), were deliberately excluded to maintain consistency and reliability in the statistical assessment. The resulting p-values were reported without correction for multiple testing. This decision was made to preserve a broader spectrum of genetic variants for subsequent evaluation, reducing the likelihood of overlooking variants that may exert meaningful biological influence. By adopting this more inclusive strategy, the analysis allowed for a more exhaustive examination of potential genotype-phenotype associations. After filtering, the dataset contained 18,970 variants, which were then annotated using the most up-to-date version of Annovar [47], aligned to the hg19 genome RefSeq [48], as of October 19, 2021. The annotated variants spanned 9,160 unique genes.

3.1 Initial model testing

In the first version, we allocated 99 % of the data for training and only 1 % for testing. While this model achieved an accuracy of 100 % with an impressively low loss value of 4.9557 × 10⁻⁵ %, these results were achieved within 6–9 training epochs. This raised concerns about overfitting, as such rapid convergence often indicates that the model may not generalize well to new data. To counter this, we explored parameter adjustments and model architecture modifications, aiming to improve its ability to generalize effectively.

3.2 Optimized data split and architectural refinements

Through a series of adjustments, we achieved more reliable results. We modified the data split to allocate 85 % for training and 15 % for testing, reduced the layer sizes, and increased the number of epochs from 15 to 20. Additionally, we expanded the number of neurons in the fully connected network (FCN) from 100 to 500. To improve clarity, we also simplified the way predictions were displayed, making performance analysis more intuitive. These refinements resulted in a more steady accuracy progression over epochs, with the loss value stabilizing at approximately 2 %. However, despite these improvements, around 50 % of the predictions remained incorrect, highlighting the need for further optimization.

3.3 Parameter optimization and additional refinements

To further enhance the model, we increased the number of epochs to 100 and introduced a new parameter: the percentage of data used for predictions. Additionally, we added an extra neuron layer to the FCN, incorporated dropout in both the max-pooling layers and the FCN, and implemented a function for graphical result visualization. These modifications led to more stable training performance, with accuracy and loss following expected trends. The testing accuracy improved to 68 %, while the prediction error rate decreased to 10 %. However, despite these enhancements, graphical analysis of the testing set still suggested the presence of overfitting, albeit at a reduced level.

3.4 Overfitting reduction and model optimization

To further minimize overfitting, we simplified the FCN by removing its second neuron layer and introduced an optimization technique. We tested both the “adam” and “SGD” optimizers and ultimately chose “adam” since the performance difference between the two was negligible. The modifications led to noticeable improvements in the testing set graphs, with smoother accuracy and loss curves showing fewer fluctuations. Accuracy saw a slight increase of approximately 1 %, while the prediction error dropped from 10 % to 8 %. Additionally, incorporating early stopping allowed the model to assess performance at each epoch and discard those that did not contribute to accuracy gains. This adjustment further stabilized the loss and accuracy curves, leading to more consistent and reliable performance.

3.5 K-Fold cross-validation and model adjustments

Following extensive testing, we introduced modifications to the model by incorporating a third dimension of size three, where each layer corresponded to a specific genotype (0, 1, and 2) [49]. Instead of using pre-trained weights, genotype-specific values were assigned, while the remaining layers were initialized with zeros. These adjustments produced accuracy and loss values comparable to previous tests, with slight improvements in learning curves. However, overfitting persisted, prompting the integration of K-fold cross-validation [50].

Even with the implementation of K-fold cross-validation, graphical analyses continued to indicate the presence of overfitting. Since the modifications did not make substantial improvements, we began exploring alternative approaches to enhance model generalization.

3.6 Final model

To further enhance performance, we experimented with various regularization methods, including Drop Blocks [51] and Spatial Dropouts. However, these approaches had minimal impact on the model’s performance. We also tested an attention mechanism [52] to improve both accuracy and interpretability, but the results did not show significant improvements. Ultimately, we opted to integrate Batch Normalization across most layers while maintaining dropout rates at 0.5.

One of the key breakthroughs came from implementing a Learning Rate Schedule, which dynamically adjusted the learning rate between epochs. This strategy helped address the early plateau observed with early stopping, leading to a notable improvement. As a result, model accuracy increased to 76 %, and the learning curves demonstrated better convergence, indicating enhanced stability during training.

In the final phase of model development, we introduced a Learning Rate Reduction function that dynamically decreased the learning rate when progress stagnated. To simplify the architecture, we used a single neuron layer with 1,000 neurons and reduced the number of filters. These refinements led to the highest recorded accuracy of 80 % (see Figure 2).

The learning curves for both the training and testing sets showed significant improvements, reflecting a more balanced ratio of correct to incorrect predictions and a noticeable reduction in overfitting. To evaluate the impact of different model configurations on performance, we tested several combinations of model dimension, filter sizes, and batch sizes. The results, presented in Table 1, highlight the most effective setups for both 2D and 3D CNNs. Additionally, we investigated the influence of varying neuron counts in different layers across both architectures. As shown in Table 2, increasing the number of neurons yielded only marginal improvements in accuracy, with no significant gains beyond a certain threshold.

Figure 2:

Performance evaluation of the final 2D model over 50 epochs, incorporating a single neuron layer with 1,000 neurons, a learning rate schedule, and a learning rate reduction function. The figure illustrates model accuracy, prediction error percentage, and learning curves for both training and testing sets.