# Gene Position Index Mutation Detection Algorithm Based on Feedback Fast Learning Neural Network

Authors: Zhike Zuo; Chao Tang; Yu Xu; Ying Wang; Yongzhong Wu; Jun Qi; Xiaolong Shi
Journal: Computational Intelligence and Neuroscience
 (2021)

---

## Abstract

In the detection of genome variation, the research on the internal correlation of reference genome is deepening; the detection of variation in genome sequence has become the focus of research, and it has also become an effective path to find new genes and new functional proteins. The targeted sequencing sequence is used to sequence the exon region of a specific gene in cancer gene detection, and the sequencing depth is relatively large. Traditional alignment algorithms will lose some sequences, which will lead to inaccurate mutation detection. This paper proposes a mutation detection algorithm based on feedback fast learning neural network position index. By establishing a position index relationship for ACGT in the DNA sequence, the subsequence is decomposed into the position relationship of different subsequences corresponding to the main sequence. The positional relationship of the subsequence in the main sequence is determined by the positional relationship. Analyzing SNP and InDel mutations, even structural mutations, through the position correlation of sequences has the advantages of high precision and easy implementation by personal computers. The feedback fast learning neural network is used to verify whether there is a linear relationship between two or more positions. Experimental results show that the mutation points detected by position index are more than those detected by Bcftools, Freebye, Vanscan2, and Gatk.


---

## Body

## 1. Introduction

In recent years, the research scope of chemical genomics has gradually expanded, combining combinatorial chemistry, cell molecular biology, and genetics to form a fusion technology model, and using high-throughput screening technology to conduct data analysis from a more comprehensive perspective. It is precise because the chemical genome analyzes life sciences from a new perspective; it has a certain promotion value. During the operation of chemical genomics technology, the specificity between the small molecule compound probe and the target protein can be used as an entry point to complete gene transcription operations, gene processing operations, and translation procedures, thereby enabling in-depth regulation of specific life processes.Currently, an index is established for DNA sequences to improve the speed of searching and matching DNA sequences. DNA target regions can generally be divided into repetitive fragments and specific fragments from the composition of the sequence. Repetitive fragments refer to the sequence in the target area where there are more repeated sequences. Literature [1] proposes Hamming distance or PFD filter to find the repetitive area, but the algorithm efficiency is low, the memory is large, and the running time is long. Literature [2] proposes to search for repetitive sequences by indexing subsequent arrays, but the efficiency is still relatively low. In the method of searching for specific fragments, the specific region segment proposed in literature [3] is composed of four bases A, G, T, and C to form an irregular sequence combination. Literature [4] studies the hash index structure of the one-way hash function and the index retrieval method to search for specific fragments and similar sequences. Literature [5] also proposed a novel solution for searching for specific DNA sequences. For the construction of the hash index structure, in DNA sequence matching [6], the commonly used fixed sequences are stored in the DNA database, and the similarity is used to evaluate whether the sequences are matched successfully. Literature [7] is applied to the index values of pattern characters and subsequence characters and matches from left to right. It stores all the index values of all characters and checks the first character of the pattern, which character appears first in the pattern as the starting matching position. Literature [8] proposes to establish index structure in DNA and protein sequences, design a multithread matching model based on DNA sequence index, match sequences of multiple tasks simultaneously and has high efficiency in DNA matching accuracy. Literature [9] proposes a hash function based on Hash-q to eliminate conflicts, providing a perfect and efficient hash value generation method. Under the condition that q characters of pattern and text need not be compared, Hash-q has better performance in accuracy and time compared with Escherichia coli and human chromosome data sets. Literature [10] applies to multiple patterns matching of DNA sequences, uses index tables for strings and patterns, and uses the count variable to count the number of occurrences of each character. The character displayed with the smallest number in the pattern is considered the first choice for comparison, and the character with the largest number in the given string is matched first. In the second part of the article, the implementation method and comparison algorithm of feedback fast learning neural network are explained. The third part describes the realization process of DNA sequence positional relationship. In the fourth part, the advantages of the proposed method are verified by experiments.


## 2. Feedback Fast Learning Neural Network and DNA Comparison Algorithm

### 2.1. Fast Learning Network (FLN)

FLN [11, 12], as a new type of double parallel and feed-forward multilayer artificial neural network, has the advantages of compact network scale, short learning and training time, strong fitting ability, and so forth. FLN consists of three layers of neurons, including input layer, hidden layer, and output layer neurons; see Figure 1.Figure 1
Diagram of the fast learning network.The weight matrix from the input layer to the hidden layer isWm×n, the weight matrix from the hidden layer to the output layer is Vp×m, and the weight matrix from the input layer to the output layer directly without passing through the hidden layer is Up×n. FLN network can be described by (1)Lm×1=fθm×1+Wm×nXn×1,Yp×1=gδp×1+Vp×mLm×1+Up×nXn×1,where n, m, and p are the number of neurons in the input layer, hidden layer, and output layer, respectively; Xn×1, Lm×1, and Yp×1 are the input vector, hidden layer output vector, and output layer output vector of the network, respectively; θm×1 and δp×1are the hidden layer threshold vector and the output layer threshold vector, respectively; and f and g are the kernel functions of the hidden layer and the output layer, respectively.


### 2.2. Feedback Fast Learning Network

FLN is a feed-forward neural network, and its output is only related to the input of the network at the current moment [13, 14] but has nothing to do with the input and output of the network at the previous moment, that is, FLN ignores the connection between the output of the system at this moment and the previous output of the system. However, for systems with large inertia or delay, the output of the model is not only related to the input at the current moment but also related to the input at the previous moments, and the input at the current moment affects the output at the following moments of the model [13, 15, 16]. Based on this, a feedback fast learning network (B-FLN) is proposed to improve the performance of FLN by adding a delayed feedback channel from output to input on the basis of FLN [17, 18]. The structure diagram of B-FLN is shown in Figure 2.Figure 2
Diagram of feedback fast learning network.In the figure, Z−1 is a delayed feedback link and Bm×p is the weight of the network output Yt−1to the hidden layer neurons at the previous moment. The B-FLN mathematical model is described as follows, where T is the current time:(2)Lm×1t=fθm×1+Wm×nXn×1t+Bm×pYp×1t−1,Yp×1t=gδp×1+Vp×mLm×1t+Up×nXn×1t.For B-FLN, if the weightsW, B, V, and U of the network are determined, thresholdsθ, δand f, g of any input sample Xn×1 correspond to an output vector Yp×1. Generally speaking, the output layer of the three-layer network takes the linear output function, while the kernel function of the hidden layer takes the “Sigmoid” function. The weights from the input layer to the hidden layer of the B-FLN and any element of the threshold value W, B, and θ are initialized to random values within [0, 1].The weightsV and U and thresholds δ from hidden layer to output layer of B-FLN are solved by Moore–Penrose generalized inverse theory. Assuming Xt−i and Y^t−ii=0,…,N−1 that the input and output sample sequence collected sequentially from a certain system and N is the length of the sample sequence, the predicted output Ytof the B-FLN model is(3)Lt=fθI+WXt+BY^t−1,Yt=gδI+VLt+UXt,where I=1,1,…,11×N.Minimize predictedYtand true valuesY^t:(4)gU^,V^,δXtYtI−Y^t=minU,V,δYt−Y^t.B-FLN weight and threshold determined by Moore–Penrose generalize inverse:(5)U^,V^,δ=g−1Y^tXtYtI†,where the superscript symbol “†” represents the M-P generalized inverse of the matrix.The predicted outputYt of the B-FLN network can be described as a function of Xt and Yt−1, and Yt−1 can be described as a function of Xt−1 and Yt−2. Therefore, B-FLN actually establishes the mapping relationship from sequence Xt, Xt−1 to Yt. The learning training of the B-FLN network is carried out by using the collected data samples, and the weights and thresholds of the network are determined. The steps are as follows:(1)
Randomly initialize the weightsW and B and threshold vector θ(2)
The predicted valueYt is calculated by using equation (3)(3)
The weightsV and U and the threshold vector δ are calculated using equation (4)


### 2.3. DNA Comparison Algorithm

The dynamic programming algorithm was also used in the DNA sequence alignment algorithm in the early days, which is a global alignment algorithm. The Needleman–Wunsch algorithm proposed by Satra et al. [19] was first used in biological sequence alignment algorithms. It has been widely used in many fields [20–22]. The basic idea of the dynamic programming algorithm is to score a given sequence, taking the sequence S = “ACGTACAAAT,” T = “ACGGTAG” as an example, as shown in the following formula:(6)ϕS′i,T′j=+1,S′i=T′j,0,S′i≠T′j,−1,S′i=−orT′j=−.Step 1.
Initialize the sequence alignment matrix.
ConstructS+1×T+1 initialization matrix V, as shown in the following formula:(7)V0,0=0,Vi,0=Vi−1,0+ϕS′i,−,1≤i≤S,V0,j=V0,j−1+ϕ−,T′j,1≤i≤T.
To initialize the matrixV, the effect is shown in Figure 3.Figure 3
Initialization matrix of the Needleman–Wunsch algorithm.Step 2.
Fill the matrix.
The current matrix value is equal to the maximum value among diagonal, horizontal, and vertical positions, and the matrix is filled numerically by the following formula:(8)Vi,j=maxVi−1,j−1+ϕS′i,T′j,Vi−1,j+ϕS′i,−,Vi−1,j+ϕ−,T′j.
After formula (8) to fill the matrix Figure 1, the effect is shown in Figure 4.Figure 4
Effect picture of the Needleman–Wunsch algorithm after filling.Step 3.
Backtracking on the matrix
The optimal backtracking position for the global sequence alignment is the maximum value ofV (|S|, |T|) in the lower right corner of the matrix; from the diagonal, vertical, and horizontal directions of V (|S|, |T|) to V (0, 0), mark the optimal global path with the identifier “⟶,” and finally form the optimal global path. The effect is shown in Figure 5.
After optimization according to the path, the global optimal alignment sequence is obtained and the effect is shown in Figure6.Figure 5
Needleman–Wunsch backtracking path.Figure 6
Effect of global optimal ratio pair.


## 2.1. Fast Learning Network (FLN)

FLN [11, 12], as a new type of double parallel and feed-forward multilayer artificial neural network, has the advantages of compact network scale, short learning and training time, strong fitting ability, and so forth. FLN consists of three layers of neurons, including input layer, hidden layer, and output layer neurons; see Figure 1.Figure 1
Diagram of the fast learning network.The weight matrix from the input layer to the hidden layer isWm×n, the weight matrix from the hidden layer to the output layer is Vp×m, and the weight matrix from the input layer to the output layer directly without passing through the hidden layer is Up×n. FLN network can be described by (1)Lm×1=fθm×1+Wm×nXn×1,Yp×1=gδp×1+Vp×mLm×1+Up×nXn×1,where n, m, and p are the number of neurons in the input layer, hidden layer, and output layer, respectively; Xn×1, Lm×1, and Yp×1 are the input vector, hidden layer output vector, and output layer output vector of the network, respectively; θm×1 and δp×1are the hidden layer threshold vector and the output layer threshold vector, respectively; and f and g are the kernel functions of the hidden layer and the output layer, respectively.


## 2.2. Feedback Fast Learning Network

FLN is a feed-forward neural network, and its output is only related to the input of the network at the current moment [13, 14] but has nothing to do with the input and output of the network at the previous moment, that is, FLN ignores the connection between the output of the system at this moment and the previous output of the system. However, for systems with large inertia or delay, the output of the model is not only related to the input at the current moment but also related to the input at the previous moments, and the input at the current moment affects the output at the following moments of the model [13, 15, 16]. Based on this, a feedback fast learning network (B-FLN) is proposed to improve the performance of FLN by adding a delayed feedback channel from output to input on the basis of FLN [17, 18]. The structure diagram of B-FLN is shown in Figure 2.Figure 2
Diagram of feedback fast learning network.In the figure, Z−1 is a delayed feedback link and Bm×p is the weight of the network output Yt−1to the hidden layer neurons at the previous moment. The B-FLN mathematical model is described as follows, where T is the current time:(2)Lm×1t=fθm×1+Wm×nXn×1t+Bm×pYp×1t−1,Yp×1t=gδp×1+Vp×mLm×1t+Up×nXn×1t.For B-FLN, if the weightsW, B, V, and U of the network are determined, thresholdsθ, δand f, g of any input sample Xn×1 correspond to an output vector Yp×1. Generally speaking, the output layer of the three-layer network takes the linear output function, while the kernel function of the hidden layer takes the “Sigmoid” function. The weights from the input layer to the hidden layer of the B-FLN and any element of the threshold value W, B, and θ are initialized to random values within [0, 1].The weightsV and U and thresholds δ from hidden layer to output layer of B-FLN are solved by Moore–Penrose generalized inverse theory. Assuming Xt−i and Y^t−ii=0,…,N−1 that the input and output sample sequence collected sequentially from a certain system and N is the length of the sample sequence, the predicted output Ytof the B-FLN model is(3)Lt=fθI+WXt+BY^t−1,Yt=gδI+VLt+UXt,where I=1,1,…,11×N.Minimize predictedYtand true valuesY^t:(4)gU^,V^,δXtYtI−Y^t=minU,V,δYt−Y^t.B-FLN weight and threshold determined by Moore–Penrose generalize inverse:(5)U^,V^,δ=g−1Y^tXtYtI†,where the superscript symbol “†” represents the M-P generalized inverse of the matrix.The predicted outputYt of the B-FLN network can be described as a function of Xt and Yt−1, and Yt−1 can be described as a function of Xt−1 and Yt−2. Therefore, B-FLN actually establishes the mapping relationship from sequence Xt, Xt−1 to Yt. The learning training of the B-FLN network is carried out by using the collected data samples, and the weights and thresholds of the network are determined. The steps are as follows:(1)
Randomly initialize the weightsW and B and threshold vector θ(2)
The predicted valueYt is calculated by using equation (3)(3)
The weightsV and U and the threshold vector δ are calculated using equation (4)


## 2.3. DNA Comparison Algorithm

The dynamic programming algorithm was also used in the DNA sequence alignment algorithm in the early days, which is a global alignment algorithm. The Needleman–Wunsch algorithm proposed by Satra et al. [19] was first used in biological sequence alignment algorithms. It has been widely used in many fields [20–22]. The basic idea of the dynamic programming algorithm is to score a given sequence, taking the sequence S = “ACGTACAAAT,” T = “ACGGTAG” as an example, as shown in the following formula:(6)ϕS′i,T′j=+1,S′i=T′j,0,S′i≠T′j,−1,S′i=−orT′j=−.Step 1.
Initialize the sequence alignment matrix.
ConstructS+1×T+1 initialization matrix V, as shown in the following formula:(7)V0,0=0,Vi,0=Vi−1,0+ϕS′i,−,1≤i≤S,V0,j=V0,j−1+ϕ−,T′j,1≤i≤T.
To initialize the matrixV, the effect is shown in Figure 3.Figure 3
Initialization matrix of the Needleman–Wunsch algorithm.Step 2.
Fill the matrix.
The current matrix value is equal to the maximum value among diagonal, horizontal, and vertical positions, and the matrix is filled numerically by the following formula:(8)Vi,j=maxVi−1,j−1+ϕS′i,T′j,Vi−1,j+ϕS′i,−,Vi−1,j+ϕ−,T′j.
After formula (8) to fill the matrix Figure 1, the effect is shown in Figure 4.Figure 4
Effect picture of the Needleman–Wunsch algorithm after filling.Step 3.
Backtracking on the matrix
The optimal backtracking position for the global sequence alignment is the maximum value ofV (|S|, |T|) in the lower right corner of the matrix; from the diagonal, vertical, and horizontal directions of V (|S|, |T|) to V (0, 0), mark the optimal global path with the identifier “⟶,” and finally form the optimal global path. The effect is shown in Figure 5.
After optimization according to the path, the global optimal alignment sequence is obtained and the effect is shown in Figure6.Figure 5
Needleman–Wunsch backtracking path.Figure 6
Effect of global optimal ratio pair.


## 3. Construction Method Based on the Location Index Topology Map Path

### 3.1. DNA Sequence Position Index

For a given target DNA sequenceS and subsequence T, use |S| to denote the length of S and |T| to denote the length of T.In the DNA sequenceS, which is all composed of A, C, G, and T, the positions of the four bases are solved, {“AAAA,” “AAAC,” “AAAG,”...,“TTTT”} for a total of 256 kinds of combination. Find all the positions of each combination in DNA, similar to the index in BWT [23, 24], and its structure is shown in Figure 7.Figure 7
Path structure diagram based on location index topology diagram.In Figure7, P11 represents the first position of the sequence “AAAA” in the DNA, and Pn11represents the n1 position of the sequence “AAAA” in the DNA. n1,n2,...,n256 are not equal, and the number of positions of each sequence in DNA is not the same. P11 points to P12, which means that, in the DNA position structure relationship, P11 is in front of P12 and the position relationship is close together.Theorem 1.
Each subsequence in the DNA sequenceS can be described by each position relationship, and the number of all positions is equal to |S|:(9)S=P11⟶P22⟶P3255⟶Pnij⋯⟶Pnkl,i,k,l=1,2,3,…,256,n1+n2+n3+⋯+n256=|S|−3.Proof.
SupposeS = “a1a2a3a4…an,” |S| = n, ai∈A,C,G,T. Divide S into a sequence of four characters, namely, “a1a2a3a4,” “a2a3a4a5,” “a3a4a5a6,”…,“an − 3an − 2an − 1an” form. The possible positional relationship of P (“a1a2a3a4”) is described as {1, 8, 51, …, n1}, the position of P (“a2a3a4a5”) is {3, 9, 10, …, n2}, and the positional relationship of P (“a3a4a5a6”) is described as {4, 12, 98, …, n3}; then, P (“an − 3an − 2an − 1an”) represents the position {43, 98, …, nn − 3}. Let P11 = {1, 8, 51, …, n1}, P22 = {2, 9, 10, …, n2}, and Pnkl = {43, 98, …, nn − 3}.
IfP (“a1a2a3a4”) ∩ {P (“a2a3a4a5”) − 1} = {1, 8, 51, …, n1} ∩ {{2, 9, 10, …, n2} − 1} ≠ Φ, then P11 ⟶ P22.
IfP (“a2a3a4a5”) ∩ {P (“a3a4a5a6”) − 1} = {2, 9, 10, …, n2} ∩ {{4, 12, 98, …, n3}-1} ≠ Φ, then P11 ⟶ P3255.
Similarly,P (“an − 4an − 3an − 2an − 1”) ∩ {P (“an − 3an − 2an − 1an”) − 1} ≠ Φ means Pik−1l−1⟶Pikl.
TheS positional relationship is described as follows.
For eachS, the decomposed sequence into a group of 4 can describe P {“AAAA”} = {3, 4, …, n1} for the “AAAA” sequence; there are a total of n1 positions in S; and P {“AAAC”} = {5, 2, …, n2} has a total of n2 positions in S for the “AAAC” sequence. Similarly, P {“TTTT”} = {14, 67, …, n256} has a total of n256 positions in S for the “TTTT” sequence. It can be seen that P {“AAAA”} + P {“AAAC”} + ⋯ + P {“TTTT”} = {3, 4, …, n1}+{5, 2, …, n2} + ⋯ + {14, 67, …, n256} = {1, 2, 3, 4, …,|S| − 3}; then, it can be described as n1 + n2 + n3 + ⋯ + n256 = |S| − 3.Theorem 2.
The DNA subsequenceT can be decomposed into a group of 4 subsequences, and the subsequences describe the sequence T through the positional relationship:(10)Pi=P1i,P2i,P3i,...,Pni,i=1,2,3,...,256.Among them,P (i) is the position information of the fourth quaternary combination, which refers to the position information corresponding to {P11, P21, P31, …, Pn11}. Whether there is a correlation between any two P (i) and P (j), that is, Pi⟶Pj, can be calculated as follows:(11)Pj=Pi−j−i.Regarding whether the sequenceT is a subsequence in S, T can be decomposed into the form of T = P (n1)P (n2) ... P (ni), if P (n1) has links with all P (ni). There are subsequences, namely,(12)Pn1=Pni−4∗ni,i=1,2,3,…,256.For example,T = “ACGAACCCCTAGAGACTAGCTAACCGGAATCAGCTA” is decomposed into T = P (7)P (6)P (93)P (115)P (113)P (91)P (14)P (46), and finally “A” is not considered. Among them, if the link of P (7)⟶, P (6)⟶, P (93)⟶, P (115)⟶, P (113)⟶, P (91)⟶, P (14)⟶, P (46) exists in S, then compare whether there is a link at the end of “A.” If the above process is established, it means that the sequence T is in the subsequence of S and the starting position of P (7) is the specific position of the subsequence T in S.


### 3.2. Mutation Detection Based on Position Index

The method of detecting mutation based on position index uses the correlation between positions to analyze whether there are mutation points in the sequence.Theorem 3.
If there are mutations such as SNP and InDel in the subsequence [25], the subsequence is decomposed into Pi=P1i,P2i,…,Pnii, i = 1, 2, 3, …, 256. If the subsequence exists in the decomposition, P (1) ⟶ P (2) ⟶ ⋯ ⟶ P (i) ↦ P (m) ↦ P (j) ⟶ ⋯ ⟶ P (n) means that the two cannot be directly connected but there is a correlation. The variation judgment formula (14) is as follows:(13)SNP,ifPj−Pi=8,Insert,ifPj−Pi<8,Delete,ifPj−Pi>8.IfP (1) ⟶ P (2) ⟶ P (i) ↦ ⋯ ↦ P (j) ⟶ P (j + 1) ⟶ P (n) forms two segments, when Pj−Pi≫8, the positional relationship between the two is far greater than 8, and the two segments may belong to different genes. This variation is called structural variation. The local comparison algorithm is used for matching to determine the mutation position, and see the following:(14)partial matchPm,SubstrT,Pi+4,Pj−1.Proof.
The following examples are used to illustrate the entire matching process, targeting the reference sequenceT = “CATCCTCACTACCT,” decomposing T [P (1)] = “CATC,” T[P (2)] = “CTCA,” T [P (3)] = “CTAC,” normal matching P (1) ⟶ P (2) ⟶ P (3) is successful, indicating that there is a match between the sequences. If P (2) mutates, it becomes P (1) ↦ P (2) ↦ P (3).
IfP (3) − P (1) = 8, SNP mutation occurs. Suppose the sequencing sequence becomes S = “CATCGTCACTACCT,” T [P (1)] = “CATC,” T [P (2)] = “GTCA” and SubstrT,P1+4,P3−1, “GTCA” and “CATC” found the fifth position is caused by “C” mutation into “G.”
IfP (3) − P (1) < 8, insert mutation occurs. Suppose the sequencing sequence becomes T = “CATCGCTCACTACCT,” T [P (1)] = “CATC,” T [P (2)] = “GCTC,” T [P (3)] = “ACTA.” Use T [P (2)] and SubstrT,P1+4,P3−1, “GCTC” and “CTC” to find that “G” is inserted at the 5th position.
IfP (3) − P (1) > 8, delete mutation occurs. Suppose the sequencing sequence becomes S = “CATCTCACTACCT,” T [P (1)] = “CATC,” T [P (2)] = “TCAC,” T [P (3)] = “TACC.” Use T [P (2)] and SubstrT,P1+4,P3−1, “TCAC” and “CTCAC” to find the fifth position and delete “C.”
Targeted sequencing is a method that refers to the sequencing of specific gene exon regions, with low cost and a sequencing depth of up to 1000. The exon regions of different cancer-targeted sequencing are different, and the raw data of targeted sequencing of different target regions are shown in Figure8.
Figure8 shows the sequencing data in two directions of the sequencing data. There are two files R1.fq and R2.fq, respectively. When the exon region is relatively short, the two sequences will overlap.Figure 8
Original data structure of DNA targeted sequencing.


### 3.3. Position Optimization of Fast Learning Neural Network Based on Position Feedback

When a position index mismatch occurs in a position index relationship, there is no linear relationship between different positions, but most positions in the position index show a position similarity relationship, so it is necessary to introduce the position information into the feedback learning neural network for learning and determine whether there is a linear relationship between two or more position index relationships by inputting the position relationship, as shown in Figure9.Figure 9
Feedback learning neural network index position relationship model.


## 3.1. DNA Sequence Position Index

For a given target DNA sequenceS and subsequence T, use |S| to denote the length of S and |T| to denote the length of T.In the DNA sequenceS, which is all composed of A, C, G, and T, the positions of the four bases are solved, {“AAAA,” “AAAC,” “AAAG,”...,“TTTT”} for a total of 256 kinds of combination. Find all the positions of each combination in DNA, similar to the index in BWT [23, 24], and its structure is shown in Figure 7.Figure 7
Path structure diagram based on location index topology diagram.In Figure7, P11 represents the first position of the sequence “AAAA” in the DNA, and Pn11represents the n1 position of the sequence “AAAA” in the DNA. n1,n2,...,n256 are not equal, and the number of positions of each sequence in DNA is not the same. P11 points to P12, which means that, in the DNA position structure relationship, P11 is in front of P12 and the position relationship is close together.Theorem 1.
Each subsequence in the DNA sequenceS can be described by each position relationship, and the number of all positions is equal to |S|:(9)S=P11⟶P22⟶P3255⟶Pnij⋯⟶Pnkl,i,k,l=1,2,3,…,256,n1+n2+n3+⋯+n256=|S|−3.Proof.
SupposeS = “a1a2a3a4…an,” |S| = n, ai∈A,C,G,T. Divide S into a sequence of four characters, namely, “a1a2a3a4,” “a2a3a4a5,” “a3a4a5a6,”…,“an − 3an − 2an − 1an” form. The possible positional relationship of P (“a1a2a3a4”) is described as {1, 8, 51, …, n1}, the position of P (“a2a3a4a5”) is {3, 9, 10, …, n2}, and the positional relationship of P (“a3a4a5a6”) is described as {4, 12, 98, …, n3}; then, P (“an − 3an − 2an − 1an”) represents the position {43, 98, …, nn − 3}. Let P11 = {1, 8, 51, …, n1}, P22 = {2, 9, 10, …, n2}, and Pnkl = {43, 98, …, nn − 3}.
IfP (“a1a2a3a4”) ∩ {P (“a2a3a4a5”) − 1} = {1, 8, 51, …, n1} ∩ {{2, 9, 10, …, n2} − 1} ≠ Φ, then P11 ⟶ P22.
IfP (“a2a3a4a5”) ∩ {P (“a3a4a5a6”) − 1} = {2, 9, 10, …, n2} ∩ {{4, 12, 98, …, n3}-1} ≠ Φ, then P11 ⟶ P3255.
Similarly,P (“an − 4an − 3an − 2an − 1”) ∩ {P (“an − 3an − 2an − 1an”) − 1} ≠ Φ means Pik−1l−1⟶Pikl.
TheS positional relationship is described as follows.
For eachS, the decomposed sequence into a group of 4 can describe P {“AAAA”} = {3, 4, …, n1} for the “AAAA” sequence; there are a total of n1 positions in S; and P {“AAAC”} = {5, 2, …, n2} has a total of n2 positions in S for the “AAAC” sequence. Similarly, P {“TTTT”} = {14, 67, …, n256} has a total of n256 positions in S for the “TTTT” sequence. It can be seen that P {“AAAA”} + P {“AAAC”} + ⋯ + P {“TTTT”} = {3, 4, …, n1}+{5, 2, …, n2} + ⋯ + {14, 67, …, n256} = {1, 2, 3, 4, …,|S| − 3}; then, it can be described as n1 + n2 + n3 + ⋯ + n256 = |S| − 3.Theorem 2.
The DNA subsequenceT can be decomposed into a group of 4 subsequences, and the subsequences describe the sequence T through the positional relationship:(10)Pi=P1i,P2i,P3i,...,Pni,i=1,2,3,...,256.Among them,P (i) is the position information of the fourth quaternary combination, which refers to the position information corresponding to {P11, P21, P31, …, Pn11}. Whether there is a correlation between any two P (i) and P (j), that is, Pi⟶Pj, can be calculated as follows:(11)Pj=Pi−j−i.Regarding whether the sequenceT is a subsequence in S, T can be decomposed into the form of T = P (n1)P (n2) ... P (ni), if P (n1) has links with all P (ni). There are subsequences, namely,(12)Pn1=Pni−4∗ni,i=1,2,3,…,256.For example,T = “ACGAACCCCTAGAGACTAGCTAACCGGAATCAGCTA” is decomposed into T = P (7)P (6)P (93)P (115)P (113)P (91)P (14)P (46), and finally “A” is not considered. Among them, if the link of P (7)⟶, P (6)⟶, P (93)⟶, P (115)⟶, P (113)⟶, P (91)⟶, P (14)⟶, P (46) exists in S, then compare whether there is a link at the end of “A.” If the above process is established, it means that the sequence T is in the subsequence of S and the starting position of P (7) is the specific position of the subsequence T in S.


## 3.2. Mutation Detection Based on Position Index

The method of detecting mutation based on position index uses the correlation between positions to analyze whether there are mutation points in the sequence.Theorem 3.
If there are mutations such as SNP and InDel in the subsequence [25], the subsequence is decomposed into Pi=P1i,P2i,…,Pnii, i = 1, 2, 3, …, 256. If the subsequence exists in the decomposition, P (1) ⟶ P (2) ⟶ ⋯ ⟶ P (i) ↦ P (m) ↦ P (j) ⟶ ⋯ ⟶ P (n) means that the two cannot be directly connected but there is a correlation. The variation judgment formula (14) is as follows:(13)SNP,ifPj−Pi=8,Insert,ifPj−Pi<8,Delete,ifPj−Pi>8.IfP (1) ⟶ P (2) ⟶ P (i) ↦ ⋯ ↦ P (j) ⟶ P (j + 1) ⟶ P (n) forms two segments, when Pj−Pi≫8, the positional relationship between the two is far greater than 8, and the two segments may belong to different genes. This variation is called structural variation. The local comparison algorithm is used for matching to determine the mutation position, and see the following:(14)partial matchPm,SubstrT,Pi+4,Pj−1.Proof.
The following examples are used to illustrate the entire matching process, targeting the reference sequenceT = “CATCCTCACTACCT,” decomposing T [P (1)] = “CATC,” T[P (2)] = “CTCA,” T [P (3)] = “CTAC,” normal matching P (1) ⟶ P (2) ⟶ P (3) is successful, indicating that there is a match between the sequences. If P (2) mutates, it becomes P (1) ↦ P (2) ↦ P (3).
IfP (3) − P (1) = 8, SNP mutation occurs. Suppose the sequencing sequence becomes S = “CATCGTCACTACCT,” T [P (1)] = “CATC,” T [P (2)] = “GTCA” and SubstrT,P1+4,P3−1, “GTCA” and “CATC” found the fifth position is caused by “C” mutation into “G.”
IfP (3) − P (1) < 8, insert mutation occurs. Suppose the sequencing sequence becomes T = “CATCGCTCACTACCT,” T [P (1)] = “CATC,” T [P (2)] = “GCTC,” T [P (3)] = “ACTA.” Use T [P (2)] and SubstrT,P1+4,P3−1, “GCTC” and “CTC” to find that “G” is inserted at the 5th position.
IfP (3) − P (1) > 8, delete mutation occurs. Suppose the sequencing sequence becomes S = “CATCTCACTACCT,” T [P (1)] = “CATC,” T [P (2)] = “TCAC,” T [P (3)] = “TACC.” Use T [P (2)] and SubstrT,P1+4,P3−1, “TCAC” and “CTCAC” to find the fifth position and delete “C.”
Targeted sequencing is a method that refers to the sequencing of specific gene exon regions, with low cost and a sequencing depth of up to 1000. The exon regions of different cancer-targeted sequencing are different, and the raw data of targeted sequencing of different target regions are shown in Figure8.
Figure8 shows the sequencing data in two directions of the sequencing data. There are two files R1.fq and R2.fq, respectively. When the exon region is relatively short, the two sequences will overlap.Figure 8
Original data structure of DNA targeted sequencing.


## 3.3. Position Optimization of Fast Learning Neural Network Based on Position Feedback

When a position index mismatch occurs in a position index relationship, there is no linear relationship between different positions, but most positions in the position index show a position similarity relationship, so it is necessary to introduce the position information into the feedback learning neural network for learning and determine whether there is a linear relationship between two or more position index relationships by inputting the position relationship, as shown in Figure9.Figure 9
Feedback learning neural network index position relationship model.


## 4. Experimental Analysis of Targeted DNA Sequencing SNP Discovery Algorithm

### 4.1. Preexperiment Processing Flow

Since targeted sequencing is performed for specific gene exons, the sequenced sequence is shorter and the sequencing depth is deeper (the test depth is 1000). The range of the target sequence is shown in Table1.Table 1
Target sequence range.
NumberChrStartEndWidthExon1chr2214728618214728776159BARD12chr2214728732214728909178BARD13chr221472887028734619177BARD1………………702chr222873443229126542188CHEK2703chr222873458028734761182CHEK2In Table1, width is the sequencing width of the sequencer, which already includes the range of exons and some introns. The length of the illumina sequencing sequence is 150. Because it is paired-end sequencing, when the exons in the width are very short, paired-end sequencing will cause overlap. In Table 1, it is shown that the sequencing herein is flux-targeted DNA sequencing for 20 genes. The sequencing target region covers all coding regions of 20 genes, exon-intron junction (20–50 bp), and part of intron region of BRCA1/2 gene, with a total of 703 exon regions.


### 4.2. Sequence Alignment Software Read Quantity Comparison

In the early stage, the performance of Bowtie2 [26], BWA [27], Hisat2 [28], and Subread [29] was compared, and the number of comparisons was compared with the number of reads based on the location relationship matching algorithm. Count the number of reads of 703 exons, as shown in Table 2.Table 2
Exon read number statistics.
ChrStartEndBWASubreadBowtie2Hisat2Posindexchr2214728618214728776260249257129243chr22147287322147289091064105110555791020chr2214728870214729046501495489489561chr221473036521473055515515515580155……………………chr22287303922873055413868696867chr22287343122873447432473146321516153207chr22287344322873461940533908256820243446chr22287345802873476111265595605601969It can be seen from Table2 that Posindex, based on the positional relationship indexing algorithm, has certain advantages in most exon regions. The BWA algorithm also performs quite well in the exon regions, and the worst effect is Hisat2. From Table 2, compare the average sequencing numbers of the exon regions of 20 genes, as shown in Figure 10.Figure 10
Average number of gene exons.It can be seen from Figure10 that the Posindex algorithm has a clear advantage in the average number of exons in terms of statistics. In terms of gene CDH1 statistics, BWA, Subread, and Bowtie2 algorithms are better than Posindex algorithm; in terms of CHEK2, MAP3K1, and TLR4 statistics, Posindex algorithm is slightly inferior to BWA algorithm, better than Subread, Hisat2, and Bowtie2; in other gene penetrances, on the other hand, the Posindex algorithm has obvious advantages. Next, analyze the matching effect of the five algorithms from the overall matching rate of the number of reads. The matching rate is equal to the ratio of the number of successful matches in the target exon region to the number of matches in the Hg38 genome, as shown in Figure 11.Figure 11
Exon region matching rate.


### 4.3. Comparison of SNP and InDel Variation Quantity

In the above comparison of commonly used algorithms, BWA software has the highest overall matching rate. Taking the Bam file generated by BWA as the research object, the mutation detection software packages Varscan2, GATK, Bcftools, and Freebayes were used to detect SNP and InDel and then compared with the location index detection algorithm to evaluate the performance. The common visualization software IGV [30] can view the variation of exon regions in bam files. The following table shows the number of SNP and InDel of exons, as shown in Tables 3 and 4.Table 3
SNP number in the exon region.
ChrStartEndBcftoolsVarscan2FreebyesGatkPosindexchr111883898188408401101chr1110827940710827954610101chr1110827949710827965610101chr5568818685688206711111chr5568819925688218611111........................Table 4
InDel number in exon region.
ChrStartEndBcftoolVarscan2FreebyesGatkPosindexchr111886659188681611112chr111881380188155511111chr111881456188163911112chr13323325723233274011121chr13323325763233274011121........................In Table3, the Posindex algorithm has certain advantages in terms of the number of SNPs, and the effect of Gatk is also ideal. Many SNP points are due to the lack of advantages of other software in terms of number, and the depth of sequencing will not meet the filtering requirements.In Table4, Gatk has an advantage in terms of quantity, but the variation points in Gatk do not exist statistically in other software, and false positives are relatively high. The Posindex algorithm is statistically reasonable.Then, compare the number of SNP and InDel in the exon regions of 20 genes, as shown in Figures12 and 13.Figure 12
SNP number in the gene region.Figure 13
InDel number in the gene region.In Figure12, the Posindex algorithm has a great advantage in the statistics of SNPs in exon regions.In genes PTEN, CCND1, PALB2, and TOX3, other algorithms did not find SNP mutation points; on LSP1, BRCA2, BRIP1, STK11, BARD1, and MAP3K1 genes, the five types of algorithms have more SNP mutation points, and the Gatk effect is also ideal. Judging from the statistical number of SNP points in the overall gene region, the Gatk and Posindex algorithms are ideal in terms of detection effects, while the other three types of algorithms have average detection effects and similar numbers.In Figure13, the InDel detection quantity is similar to the SNP statistical quantity, and the Gatk and Posindex algorithm detection is ideal.Compare the positions of the overall SNP and InDel with those of other software, and the effect is shown in Figures14 and 15.Figure 14
Comparison of SNP position differences.Figure 15
Comparison of InDel position differences.By analyzing the analysis process of SNP and InDel, the proposed method Posindex can find more mutation points, and the Gatk detection effect is also ideal, while Bcftools, Varscan2, and Freebyes detect fewer SNP and InDel mutation points. A high number of detection results does not mean a good detection effect, because there will be many false positives in the detection process, so the number cannot explain the detection advantage. The difference between Gatk and position index is that the false positive rate of Gatk is high, while the correct rate of position index is high.The specific reasons are as follows:(1)
In the original sequencing data, there are a large number of reverse sequences, containing a large amount of public data and index data, which cannot be completely deleted when cleaning up.(2)
The problem of false negatives: there are many similar sequences in the DNA sequence, which may lead to relative sequencing data pointing to other locations. However, these sequences are not target sequences, resulting in a large number of false positives or false negatives, and many SNP and InDel mutation points will be lost.


## 4.1. Preexperiment Processing Flow

Since targeted sequencing is performed for specific gene exons, the sequenced sequence is shorter and the sequencing depth is deeper (the test depth is 1000). The range of the target sequence is shown in Table1.Table 1
Target sequence range.
NumberChrStartEndWidthExon1chr2214728618214728776159BARD12chr2214728732214728909178BARD13chr221472887028734619177BARD1………………702chr222873443229126542188CHEK2703chr222873458028734761182CHEK2In Table1, width is the sequencing width of the sequencer, which already includes the range of exons and some introns. The length of the illumina sequencing sequence is 150. Because it is paired-end sequencing, when the exons in the width are very short, paired-end sequencing will cause overlap. In Table 1, it is shown that the sequencing herein is flux-targeted DNA sequencing for 20 genes. The sequencing target region covers all coding regions of 20 genes, exon-intron junction (20–50 bp), and part of intron region of BRCA1/2 gene, with a total of 703 exon regions.


## 4.2. Sequence Alignment Software Read Quantity Comparison

In the early stage, the performance of Bowtie2 [26], BWA [27], Hisat2 [28], and Subread [29] was compared, and the number of comparisons was compared with the number of reads based on the location relationship matching algorithm. Count the number of reads of 703 exons, as shown in Table 2.Table 2
Exon read number statistics.
ChrStartEndBWASubreadBowtie2Hisat2Posindexchr2214728618214728776260249257129243chr22147287322147289091064105110555791020chr2214728870214729046501495489489561chr221473036521473055515515515580155……………………chr22287303922873055413868696867chr22287343122873447432473146321516153207chr22287344322873461940533908256820243446chr22287345802873476111265595605601969It can be seen from Table2 that Posindex, based on the positional relationship indexing algorithm, has certain advantages in most exon regions. The BWA algorithm also performs quite well in the exon regions, and the worst effect is Hisat2. From Table 2, compare the average sequencing numbers of the exon regions of 20 genes, as shown in Figure 10.Figure 10
Average number of gene exons.It can be seen from Figure10 that the Posindex algorithm has a clear advantage in the average number of exons in terms of statistics. In terms of gene CDH1 statistics, BWA, Subread, and Bowtie2 algorithms are better than Posindex algorithm; in terms of CHEK2, MAP3K1, and TLR4 statistics, Posindex algorithm is slightly inferior to BWA algorithm, better than Subread, Hisat2, and Bowtie2; in other gene penetrances, on the other hand, the Posindex algorithm has obvious advantages. Next, analyze the matching effect of the five algorithms from the overall matching rate of the number of reads. The matching rate is equal to the ratio of the number of successful matches in the target exon region to the number of matches in the Hg38 genome, as shown in Figure 11.Figure 11
Exon region matching rate.


## 4.3. Comparison of SNP and InDel Variation Quantity

In the above comparison of commonly used algorithms, BWA software has the highest overall matching rate. Taking the Bam file generated by BWA as the research object, the mutation detection software packages Varscan2, GATK, Bcftools, and Freebayes were used to detect SNP and InDel and then compared with the location index detection algorithm to evaluate the performance. The common visualization software IGV [30] can view the variation of exon regions in bam files. The following table shows the number of SNP and InDel of exons, as shown in Tables 3 and 4.Table 3
SNP number in the exon region.
ChrStartEndBcftoolsVarscan2FreebyesGatkPosindexchr111883898188408401101chr1110827940710827954610101chr1110827949710827965610101chr5568818685688206711111chr5568819925688218611111........................Table 4
InDel number in exon region.
ChrStartEndBcftoolVarscan2FreebyesGatkPosindexchr111886659188681611112chr111881380188155511111chr111881456188163911112chr13323325723233274011121chr13323325763233274011121........................In Table3, the Posindex algorithm has certain advantages in terms of the number of SNPs, and the effect of Gatk is also ideal. Many SNP points are due to the lack of advantages of other software in terms of number, and the depth of sequencing will not meet the filtering requirements.In Table4, Gatk has an advantage in terms of quantity, but the variation points in Gatk do not exist statistically in other software, and false positives are relatively high. The Posindex algorithm is statistically reasonable.Then, compare the number of SNP and InDel in the exon regions of 20 genes, as shown in Figures12 and 13.Figure 12
SNP number in the gene region.Figure 13
InDel number in the gene region.In Figure12, the Posindex algorithm has a great advantage in the statistics of SNPs in exon regions.In genes PTEN, CCND1, PALB2, and TOX3, other algorithms did not find SNP mutation points; on LSP1, BRCA2, BRIP1, STK11, BARD1, and MAP3K1 genes, the five types of algorithms have more SNP mutation points, and the Gatk effect is also ideal. Judging from the statistical number of SNP points in the overall gene region, the Gatk and Posindex algorithms are ideal in terms of detection effects, while the other three types of algorithms have average detection effects and similar numbers.In Figure13, the InDel detection quantity is similar to the SNP statistical quantity, and the Gatk and Posindex algorithm detection is ideal.Compare the positions of the overall SNP and InDel with those of other software, and the effect is shown in Figures14 and 15.Figure 14
Comparison of SNP position differences.Figure 15
Comparison of InDel position differences.By analyzing the analysis process of SNP and InDel, the proposed method Posindex can find more mutation points, and the Gatk detection effect is also ideal, while Bcftools, Varscan2, and Freebyes detect fewer SNP and InDel mutation points. A high number of detection results does not mean a good detection effect, because there will be many false positives in the detection process, so the number cannot explain the detection advantage. The difference between Gatk and position index is that the false positive rate of Gatk is high, while the correct rate of position index is high.The specific reasons are as follows:(1)
In the original sequencing data, there are a large number of reverse sequences, containing a large amount of public data and index data, which cannot be completely deleted when cleaning up.(2)
The problem of false negatives: there are many similar sequences in the DNA sequence, which may lead to relative sequencing data pointing to other locations. However, these sequences are not target sequences, resulting in a large number of false positives or false negatives, and many SNP and InDel mutation points will be lost.


## 5. Conclusion

In the detection of chemical genomic mutations, a DNA sequence matching algorithm based on the position index relationship is proposed to solve the problems of low accuracy and large differences in the detection of SNP and InDel mutations in the targeted sequencing sequence, which aims to establish the DNA sequence Position index relationship analysis SNP and InDel variation. First, divide the subsequence intok fixed sequences and establish links; secondly, analyze the position difference in the optimal link and establish a judgment model of position variation; finally, target the sequencing target area to cover the BRCA1/2 gene: the entire coding region, the exon-intron junction region (20–50 bp), and part of the intron region, a total of 703 exon regions. The actual data captured in the 101.3k area is used as an example to verify. The experimental results show that the location-based indexing method detects more mutation points than Bcftools, Freebyes, Vanscan2, and Gatk. After detecting SNP and InDel in this paper, it is found that the location-based index method has the best detection performance, but whether other data sets have the same effect needs a one-step test. Sequencing in other cancers will be applied and further analysis will be carried out later.



---

*Source: 1716182-2021-07-07.xml*