R-CRISPR: A Deep Learning Network to Predict Off-Target Activities with Mismatch, Insertion and Deletion in CRISPR-Cas9 System

The Clustered Regularly Interspaced Short Palindromic Repeats (CRISPR)—associated protein 9 (Cas9) system is a groundbreaking gene-editing tool, which has been widely adopted in biomedical research. However, the guide RNAs in CRISPR-Cas9 system may induce unwanted off-target activities and further affect the practical application of the technique. Most existing in silico prediction methods that focused on off-target activities possess limited predictive precision and remain to be improved. Hence, it is necessary to propose a new in silico prediction method to address this problem. In this work, a deep learning framework named R-CRISPR is presented, which devises an encoding scheme to encode gRNA-target sequences into binary matrices, a convolutional neural network as feature extractor, and a recurrent neural network to predict off-target activities with mismatch, insertion, or deletion. It is demonstrated that R-CRISPR surpasses six mainstream prediction methods with a significant improvement on mismatch-only datasets verified by GUIDE-seq. Compared with the state-of-art prediction methods, R-CRISPR also achieves competitive performance on datasets with mismatch, insertion, and deletion. Furthermore, experiments show that data concatenate could influence the quality of training data, and investigate the optimal combination of datasets.


Introduction
The CRISPR-Cas9 (Clustered Regularly Interspaced Short Palindromic Repeats (CRISPR)-CRISPR-associated protein 9 (Cas9)) system is a robust genome-editing tool with a broad range of applications in numerous research [1][2][3]. After the recognition of the 3-nucleotide protospacer adjacent motif (PAM), the endonuclease Cas9 uses a single guide RNA (gRNA) to form base pairs with any DNA target sequences of interest and introduce a site-specific double-strand break [1,4,5]. The high-efficiency and simplicity of CRISPR-Cas9 system enabled genome engineering has great potential in improving agriculture productivity and clinical application [6,7].
The CRISPR-Cas9 system is widely used to enable highly efficient genome editing in various species and cell types, but it may wrongly bind to the unwanted region and cause extra off-target activity. These off-target activities can confound research experiments and also affect the practical application of the technique [8]. The Cas9 can be programmed by altering the sequence of gRNA to target abundant sites in the genome, and the off-target effects of different gRNAs may vary greatly [9]. Therefore, it is crucial to design the offtarget prediction model to evaluate the on-and off-target activities of gRNA and choose gRNA with high on-target rate and low off-target effect [10].
From the perspective of gRNA binding to non-target regions, the off-target activities induced by CRISPR-Cas9 mechanism can be divided into three categories: (a) nucleic acid base mismatch with on-target sites; (b) nucleic acid base insertion from gRNA sequence; proposed CRISPR-Net [32], in which the Inception module that combined several kernels with different sizes were used as feature extractor in the convolutional layer, and the Long Short-Term Memory (LSTM) units were used to form a recurrent neural network in terms of its advantages of selective memory function. Although the method uses a data feature extractor to prevent information loss, it still needs to be further improved to preserve the original information. Meanwhile, since those existing prediction methods cannot satisfy enough precision for implementing CRISPR/Cas9 gene-editing techniques at the clinical level, it is pressing to propose a new method to address the problem.
In this work, we propose an off-target prediction model based on a recurrent convolutional network named R-CRISPR, predicting off-target activities of gRNA-target sequence with mismatch, insertion, and deletion. We first encode the target sequence pair into a binary matrix as the input of the prediction model and then use the preprocessing module on the basis of the RepVGG to extract data features. Finally, the bi-directional recurrent network constructed by Long Short Term Memory units is used for further training of data to improve learning efficiency.
This work provides the following contributions: 1. We developed R-CRISPR, a recurrent convolutional network to evaluate and predict off-target effects of gRNA-target sequence with mismatch, insertion, and deletion.
2. We compare the R-CRISPR with five mainstream prediction methods on datasets obtain from experimental methods to evaluate the model performance. Using the area under the curve of Receiver Operating Characteristic Curve (ROC) and Precision Recall Curve (PRC) as the measurement standard, the performance of R-CRISPR surpasses existing mainstream prediction models.
3. We compare the R-CRISPR with the state-of-art prediction model CRISPR-Net, the R-CRISPR model has an improvement of 0.2% and 1.9% on AUROC and AUPRC. 4. We make extended research to explore the performance difference on various combinations of training datasets, and improve the prediction accuracy by designing an ideal dataset combination.

Datasets
Seven off-target datasets that were validated by mainstream experimental methods were selected for model training and validation [31]. Those datasets were shifted into two categories: one category contains mismatch, insertion, and deletion off-target sites, while another just includes datasets with mismatch off-target sites.
As shown in Table 1, the total sites denote the total number of active off-target sites and inactive off-target sites obtained from Cas-Offinder [33], which search for potential binding targets for Cas9 RNA-guided endonucleases by given gRNA sequences. Dataset CIRCLE contains mismatch, insertion, and deletion, and was confirmed by in vitro method CIRCLEseq [17]. The highly sensitive unbiased method CIRCLE-seq is based on the principle of detecting new DNA cleavage events, acquired from purified circularized genomic DNA treated with Cas9:gRNA complex, by high-throughput analysis. A total of 7371 off-target sites were validated by CIRCLE-seq with 430 insertion and deletion sites, besides, from the 10 gRNA sequences, Cas-Offinder obtained 577,578 inactive off-target sites.  CIRCLE  -584,949  7371  430  10  PKD  -4853  2273  -65  PDH  -10,129  52  -19  SITE  -217,733  3767  -9  GUIDE_I  -294,534  354  -9  GUIDE_II  -95,829  54  -5  GUIDE_III  -383,463  56  -22 The datasets in the second category include mismatch sites only. Based on the protein knockout detection method, Dataset PKD was constructed by targeting on human coding sequence CD33 [9], which is composed of 4853 target sites and 2273 off-target sites. Confirmed by Polymerase Chain Reaction (PCR) amplification technology, Digenome-Seq and HTGTS [34], Haeussler et al. constructed Dataset PDH includes 10,129 target sites and 52 off-target sites. Dataset SITE contains 3767 positive off-target sites validated by SITE-seq [35] and 217,733 gRNA-target pairs in total. Tsai et al. constructed Datasets GUIDE_I, GUIDE_II and GUIDE_III [12], Kleinstiver et al. [8] and Listgarten et al. [23] based on the cellular method GUIDE-seq, which required individual transfections for each target. With 294,534 recognized sites in total, Dataset GUIDE_I contains 354 off-target sites with mismatch. And on the basis of the 5 and 22 gRNAs verified by GUIDE-seq, Datasets GUIDE_II and GUIDE_III include 54 and 56 active off-target sites, respectively. In the experimental section, Datasets CIRCLE, PKD, PDH, SITE, GUIDE_I were used as training datasets, while Datasets GUIDE_II and GUIDE_III served as test datasets.

R-CRISPR Model
The construction of R-CRISPR mainly includes three stages. In the first step, the input on-target and off-target sequences are encoded into binary matrices by an encoding scheme. The output of the encoding scheme is then transmitted into a convolutional layer composed of convolutional kernels and RepVGG blocks for data features extraction. Finally, the output of the feature extraction layer is passed to the bi-directional recurrent layers based on LSTM units to learn sequential patterns.

Encoding Matrix Scheme for gRNA-Target Pair
Suppose X on = [x on1 , x on2 , . . . , x onn ] could represents the on-target sequences, while x o f f n demonstrates the off-target sequence, where n denotes the length of sequences. Since the off-target activities could be divided into nucleic base mismatch, nucleic base insertion and nucleic base deletion as shown in Figure 1, x on and x o f f could be represented by x on , x o f f ∈ {A, C, G, T, _} where the symbol "_" denotes insertion or deletion. In terms of the above thoughts, each base in the sequence can be represented by a five-bit vector by one-hot encoding mechanism and the gRNA-target base pair could be encoded into a ten-bit vector as the suitable input of convolutional neural network (e.g., "0100000100" represents the mismatch "C → G"). However, since the off-target sites were analogous to the on-target sites with only difference on mismatch, insertion and deletion site, the encoding scheme could be further optimized. In the scheme, five-bit vectors (A, C, G, T, _) retained the nucleotides of each base pair, and two-bit vectors were used to represent the base pair type (i.e., match, mismatch, insertion, deletion). The seven-bit scheme not only reduces the input size of the neural network but preserves the various information of gRNA-target pairs. Mismatch "T → A" was encoded into "1001001", where "10010" represented the base pair included A and T, while "01" referred to A as the off-target site and T was the on-target site. Similarly, mismatch "A → T" could be encoded into "1001010", insertion "_ → A" could be presented by "1000101", and deletion "A → _" could be regarded as "1000110". As a result, every gRNA-target sequence could be represented by a 7 × 24 matrix E = [x 1 , x 2 , . . . , x T ] where 24 is the length of sequence that includes 3-bp PAM adjacent and 21 bases.

Preprocessing Module for Feature Extraction
The classic convolutional neural network VGG [36] achieved excellent performance in image recognition, which uses several 3 × 3 kernels to replace the larger ones and with a simple architecture composed of the convolutional kernel, ReLU activation, and pooling. To improve recognition accuracy, more complicated and well-designed architectures such as ResNet [37], Inception [38] were introduced into the area of computer vision. Though many complicated architectures deliver higher accuracy, there still exist significant drawbacks such as limited implementation and reduce memory utilization.
As shown in Figure 2, a RepVGG [39] block was used as the preprocessing module of R-CRISPR, which had the advantages of multi-branch designs and plain topology designs, to discover useful features and avoid biases introduced by hand-crafted rules. Inspired by ResNet, the structure of RepVGG block includes a 3 × 3 kernel, a 1 × 1 kernel and an identity branch, it becomes is a convolutional shortcut implemented by a 1 × 1 kernel, and is an identity branch.  [39]. A RepVGG module was comprised by a 3 × 3 kernel branch, a 1 × 1 kernel branch and an identity branch.

Long Short-Term Memory for Constructing RNN
LSTM [40] is a variant of RNN proposed to solve long-term dependencies problem (i.e., gradient explosion and gradient vanishing) while memorizing long-range information from sequence [41,42]. Meanwhile, LSTM layer is capable of automatically regulating self-connecting loops to memorize long-range information more effectively, since gene sequences could be regarded as the language of biology, such characteristic process significant advantage in learning sequences features.
LSTM composed of two states (i.e., c t and h t ), and three gates (i.e., input gate i t , forget gate f t , and output gate o t ). For each stage, the neuron of neural network provides input x t at time t, previous cell state c t−1 at time t − 1, and previous hidden state h t−1 at time t − 1.
The key equations of LSTM unit are as follows: and U ∈ R n×n refer to the weight matrix, while h t ∈ R n is the hidden state that uses n to represent number of the hidden states, y t is the output at time t. Initial value of c 0 and h 0 is 0 while the operator " * " denotes Hadamard product.

R-CRISPR Model Construction
The Long-term Recurrent Convolutional Neural Network (LRCN) that combines CNN and recurrent neural network achieved huge success in the areas of speech recognition and machine translation. Recently, LRCN architecture was also introduced into bioinformatics, and it is approved that the LRCN framework outperformed the CNN and RNN architectures on prediction of transcription factor binding site [33].
Off-target prediction model R-CRISPR was inspired by the LRCN framework and includes an encoding scheme [31] to convert the on-and off-target pair into suitable input for neural network, a convolutional layer, and a recurrent layer. The convolutional layer built on the architecture of CNN and RepVGG [39] module is used as a feature extractor, while the recurrent layer is composed of bi-directional LSTM RNN, and the output of the recurrent layer is passed to the subsequent dense layers. Figure 3 describes the network architecture of R-CRISPR.
On-and off-target sequence pair (X on , X o f f ) was regard as the input of R-CRISPR and passed to the encoding mechanism to be encoded into a binary matrix 1} 7 , and T referred to the length of the on-and off-target sequence pair. The matrix E was then transferred to the convolutional layer comprised forty convolutional filters with size of 1 × 1, learning a representation ∆(·) after convolution and batch normalization operation, and produce matrix C = ∆(E) = [c 1 , c 2 , . . . , c T ] where c t ∈ R 40 . was then proceed to the forty RepVGG module that comprised by 3 × 3 kernel, 1 × 1 kernel and identity branch, learning a special representation Φ(·) on C and output R i = Φ i (C) = [r 1 , r 2 , . . . , r T ] where r t ∈ R 40 . In view of the theory of structural re-parameterization, the input of bi-directional recurrent network G = R 1 + R 2 + R 3 and produce G = [g 1 , g 2 , . . . , g T ] where g t ∈ R 40 .
In order to obtain better analysis of the sequence features extracted by the preprocessing module, the recurrent layer was designed to combine two directional RNNs in which contains 15 LSTM units to learn forward patterns or backward patterns. For the forward direction, each LSTM unit maps the input g t and the previous hidden state h t−1 to produce the output y  Overview of R-CRISPR. The input on-and off-target sequence pairs is represented into a 7 × length (i.e., length represents the length of gRNA-target sequences) binary matrix as the input of the feature extraction layer. The feature extraction layer contains 40 convolutional kernels and 40 modified RepVGG modules composed of a 1 × 3 convolutional kernel, a 1 × 1 convolutional kernel, and an identity branch. The output of the feature extraction layer is then passed to the bi-directional recurrent layers, each direction is based on 15 LSTM units, to learn sequential patterns of the feature matrix. Followed by the recurrent layer, there are two dense layers with sigmoid as the activation function for final outputs.
The main task of R-CRISPR is to predict the on-and off-target effects, which could be seen as a binary classification task. The labels of the off-target sequence were labeled with label 1, while 0 could represent the other non-off-target sequences. And the Cross Entropy Loss Function was used as the loss function of this model.
the y refers to the distribution of true label, while the a refers to the distribution predicted after training. The Cross Entropy Loss Function could be used to measure the similarity between y and a, as well as weight update tardiness caused by quadratic loss function when sigmoid is used as activation function.

Mainstream Prediction Methods
In the next experimental section, six mainstream in silico prediction methods will be selected to make a comparison with our model R-CRISPR. As the groundbreaking prediction method based on machine learning, CFD [9] constructed a Naïve Bayes to predict off-target activities and surpassed the hand-crafted rules models. The widely recognized regression model Elevation-score [23] contained two layers, the first layer using boosted regression tree to predict the off-target score for single mismatch of gRNAs, while the second layer constructed an L1-regularized linear regression combiner model to calculate the aggregate score of a single gRNA from multiple off-target activities related to it. With the same accuracy, the training speed of Ensemble SVM [25] was greatly improved, which made it comparatively more suitable for large datasets. In CNN_std [28], each sequence was encoded into a matrix as input firstly, and then multiple sizes of filters were used to capture features in different ranges, the feature matrix was passed several convolutional layers and a dense layer to learn sequential patterns. Also based on CNN architecture, AttnToMismatch_CNN [29] introduced an attention mechanism to select the information that was highly correlated to off-target activities as whole gRNAtarget sequence information. The state-of-art prediction model CRISPR-Net [32] combined the advantages of the inception module and LSTM units, which had achieved higher performance accuracy than previous models.

Results
In the training time, Adam optimizer dynamically optimized the learning rate to achieve both efficiency and effectiveness and the initial learning rate of weight was set as 0.0001. Besides, the batch size of each batch was set at 10,000 with the epoch number was set as 100. To systematically represent the performance of prediction models, ROC (Receiver Operating Characteristic curve) and PR (Precision-Recall) analysis were used as an evaluation criterion. As Table 2 shows, the hyperparameters are as follows: Besides, all components of R-CRISPR were implemented using Keras 2.2.4 with TensorFlow 2.3.0 backend.

Performance of R-CRISPR on Mismatch-Only gRNA-Target Prediction
In terms that base mismatch occupies a large proportion in three kinds of off-target types (i.e., nucleic acid base mismatch, nucleic acid base insertion, and nucleic acid base deletion), and the existed mainstream prediction methods were designed to predict offtarget sites with mismatch, we first verified the performance of R-CRISPR with six models (i.e., AttnToMismatch_CNN, Elevation-score, CFD, Ensemble SVM, CNN_std and CRISPR-Net) on mismatch-only datasets. Using the combination of datasets PKD, PDH, and GUIDE_I as training set and tested on dataset GUIDE_II. R-CRISPR achieved relatively highest performance both AUROC (Area under ROC curve) and AUPRC (Area under PR curve), with an accuracy of 0.991 on AUROC and 0.319 on AUPRC. As shown in Table 3, the difference is relatively slight on AUROC while AUPRC score appeared significant differences between diverse models, the R-CRISPR appeared maximum value of 0.319, and minimum value is 0.071 of AttnToMismatch_CNN. Though the AUROC score of R-CRISPR (0.991) was slightly lower than Ensemble SVM (0.993) and CRISPR-Net (0.993) on the GUIDE-seq dataset, R-CRISPR held an improvement of 18.8% and 2.7% on AUPRC.

Performance of R-CRISPR on Multiple gRNA-Target Prediction
In the previous section, we had evaluated the performance of R-CRISPR on mismatchonly datasets and proved that R-CRISPR outperformed the six existing models in the previous experiment, during the second stage, we explored how nucleic acid base insertion and deletion affect prediction accuracy and made comparison with the state-of-art offtarget prediction method CRISPR-Net. CRISPR-Net is built upon a long-term recurrent convolutional neural network and could recognize off-target activities with base mismatch, insertion and deletion. Moreover, Elevation-score was served as a benchmark to better evaluate model performance.
Since Dataset CIRCLE was the only dataset that contained three categories of off-target activities, three models were evaluated with 5-fold cross-validation. For each validation, one subset was used as the test dataset and the other four subsets were served as the training dataset. Figure 4 shows that compared with CRISPR-Net, though CRISPR-Net represented a tiny higher accuracy on AUROC (0.1%), R-CRISPR achieved an improvement of 4.1% on AUPRC. Combination of datasets CIRCLE, PKD, PDH, SITE, and GUIDE_I as training datasets by concatenating, to preserve the biological information of insertion and deletion while adding more mismatch sites. As shown in the Figure 5, drafting ROC curve and PR curve based on the prediction result on dataset GUIDE_II, R-CRISPR (AUROC = 0.991, AUPRC = 0.312) outperformed than CRISPR-Net (AUROC = 0.993, AUPRC = 0.297) on AUPRC with an improvement of 2.2%, and also surpassed the Benchmark (AUROC = 0.993, AUPRC = 0.131) on AUPRC with an improvement of 18.1%. Furthermore, we believe that the data concatenate may affect the quality of training data and improve model performance. Thus, we investigated various combinations of datasets to improve the performance of R-CRISPR in the next section.

Performance of R-CRISPR with Different Training Datasets
In previous study, we figured out that the model predictive performance could be influenced by the quality of training datasets. Thus, we generated seven training datasets (see Table 4) from five experimental datasets (i.e., datasets CIRCLE, PKD, PDH, SITE and GUIDE_I) in which the active off-target sites were validated by CIRCLE-seq, Digenome-seq, SITE-seq and GUIDE-seq. Testing on dataset GUIDE_II, seven R-CRISPR models represented competitive performance on the ROC curve as shown in Figure 6, with an average AUROC of 0.992. However, the test results of seven training datasets were numerous on the PR curve, dataset B (Combination of datasets PKD, PDH, SITE and GUIDE_I) achieved the highest AUPRC of 0.319, while dataset C (Only includes dataset CIRCLE appeared the lowest AUPRC score of 0.173. The result indicated that the designing of training dataset could improve predictive performance significantly, which may be because those datasets achieved higher accuracy also contains more abundant gRNAs and off-target sites. The R-CRISPR trained on combined datasets surpassed those trained on a single dataset among all seven models on dataset GUIDE_II. The model trained on dataset B reached the highest AUPRC of 0.319, which possessed an advantage of 0.7% on AUPRC compared to the second best model (AUROC = 0.991, AUPRC = 0.312) trained on dataset D (Combination of dataset CIRCLE, PKD, PDH, SITE and GUIDE_I), and had an improvement of 5.4% on AUPRC compared to the third best model (AUROC = 0.992, AUPRC = 0.265) trained on dataset F (Combination of dataset PKD, PDH, and GUIDE_I). As shown in Figure 7, For further exploring the efficiency of models trained on datasets B, D and dataset F, we tested those models on dataset GUIDE_III in which concludes 56 off-target sites and 22 diverse gRNAs. Table 5 shows it is obvious that the model trained on dataset B could achieve better performance (AUROC = 0.998, AUPRC = 0.184), and appeared an improvement of 0.4% and 3.4% on both AUROC and AUPRC than the second best model (AUROC = 0.994, AUPRC = 0.150) trained on dataset F.

Hyperparameters Optimization
The optimization process of large-scale machine learning usually contained a large number of hyper parameters that needed to be fixed by users according to a certain application, and the design of hyper parameters could directly influence the model performance. In this optimization section, we would like to explore hyperparameters combination that could achieve higher performance based on five kinds of hyper parameters (i.e., dropout_rate, learning_rate, batch_size, and epochs).
Given a set of hyperparameters and the potential assignments from its parameter space, the fundamental method Grid Search was used as the search practice to select the combination that outperformed others. Furthermore, we selected the Dataset CIRCLE as test data since it contains most off-target activities as well as various off-target categories, and AUROC was used to evaluate the certain performance of hyperparameters combinations. As Figure 8 shows, the best combination achieved 0.98877 on AUROC, in which dropout_rate = 0.5, learning_rate = 0.001, batch_size = 10,000 and epochs = 50. Significantly, learning_rate was inappropriate to set too high, while dropout_rate and epochs were not suitable to be too low.

Discussion
The accurate evaluation of off-target activities in the CRISPR-Cas9 system is a severe issue when applying machine learning. Since the early prediction models remained handcrafted rules and limited predictive accuracy. In this study, we first used an encoding scheme to encode each gRNA-target sequence into a 7 × 24 matrix as the input of an improved convolutional neural network for data feature extraction. Then, given the above strategies, we proposed R-CRISPR, an off-target prediction model based on a recurrent convolutional network with a Cross Entropy Loss Function to solve the problem. Since the mainstream in silicon off-target activities prediction methods lacked further research on gRNA-target pairs insertion and deletion problems, we optimized R-CRISPR to satisfy the demands of insertion and deletion detection. We first explored the prediction accuracy of mismatch problems in terms that nucleic acid base mismatch occupies the main proportion of off-target sites and most existing predictive methods were designed for mismatch-only problems. On mismatch-only off-target dataset GUIDE_II verified by GUIDE-seq, experiments show that R-CRISPR outperformed six existing mainstream predictive methods on both ROC and RC analysis with an average accuracy of 0.991 on AUROC and 0.319 on AUPRC. In addition, we set a 5-fold cross-validation test based on the off-target dataset confirmed by CIRCLE-seq (with nucleic acid base insertion and deletion) to investigate how insertion and deletion problems affect the off-target prediction. We trained and compared R-CRISPR with the state-of-art prediction method CRISPR-Net, which could also measure off-target sites with insertion and deletion, on different combinations of datasets. R-CRISPR achieved a higher accuracy of 0.976 on AUROC and 0.460 on AUPRC with an improvement of 0.1% and 4.1% than CRISPR-Net. Furthermore, we also explored how the quality of training data is influenced by data concatenation and designed seven combinations of datasets to test the performance of R-CRISPR. Seven R-CRISPR models expressed competitive performance on ROC analysis with an average accuracy of 0.992 on AUROC, while the test results were numerous on PR analysis with the highest accuracy achieved 0.319 and lowest one appeared 0.173. The experiments indicated that the designing of training datasets could affect predictive results significantly, and the R-CRISPR trained on combined datasets surpassed those trained on a single dataset. We believed that the combination of multiple datasets could obtain multifarious information of off-target activities, and produce a more comprehensive dataset, hence improving the model performance. Meanwhile, we speculated that the sample imbalance caused by fewer positive samples was also a crucial point for model performance. Since the off-target activities only occupied a minority number in the whole biological process, the datasets obtained from most experiments were unbalanced, which required further optimization.

Conclusions
In our work, we developed R-CRISPR to contribute to the quantification of off-target activities with nucleic acid base mismatch and deletion problems. The architecture of R-CRISPR demonstrated the practicality of convolutional recurrent neural network in predicting off-target sites between gRNA sequence and target DNA sequence. Since convolutional network could be used to do preliminary information extraction, we applied the RepVGG module in the convolutional layer to capture features for the target sequence with unclear biological structure, and combined a bi-directional recurrent network based on LSTM units for further training. Furthermore, as complementary off-target sequences and related datasets become available, the efficiency and predictive accuracy are expected to be improved. We will also carefully investigate more superior model architecture based on deep learning and an optimized combination of training datasets to improve model performance. In a nutshell, the experimental results in our work fully demonstrated that R-CRISPR is an effective off-target prediction method and can contribute to the gRNA design in the CRISPR-Cas9 system.  Data Availability Statement: Data used in this article was obtained from Jiecong Lin, and they are available at https://codeocean.com/capsule/9553651/tree/v1 (accessed on 24 September 2021).

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: