Detection and classification of genetic mutations in BRCA1 gene in human chromosome 17
Abstract
Detection and classification of genetic mutations in BRCA1 gene in human chromosome 17
Abstract
One of the major tasks in clinical genomics is the identification of mutations associated with human genetic diseases. Typically, genome-wide genetic studies identify a large number of variants that have potential association with a disease. However, to narrow this search and to pinpoint the variants that most likely cause a disease, a number of methods have been developed . These methods use the evolutionary conservation of nucleotide positions and/or the functional consequences of mutations to distinguish disease-associated variants from neutral and benign variants.
Breast cancer is a common disease. Each year, approximately 200,000 women in the United States are diagnosed with breast cancer, and one in nine American women will develop breast cancer in her lifetime. In 1994, the first gene associated with breast cancer — BRCA1 (for BReast CAncer1) was identified on chromosome 17. When individuals carry a mutated form of BRCA1, they have an increased risk of developing breast or ovarian cancer at some point in their lives. Children of parents with a BRCA1 mutation have a 50 percent chance of inheriting the gene mutation.
“The simplest denomination of breast cancer is based upon inherited susceptibility to breast cancer vs sporadic occurrences of breast cancer. Heightened breast cancer risk may be due to a genetic alteration that increases susceptibility based upon an inherited heterozygous gene defect in for example BRCA1, TP53, PTEN or other tumor suppressors”
The quote is taken from DNA damage and breast cancer written by Jennifer D Davis and Shiaw-Yih Lin.
Overview
The current project is focused on applying deep learning approaches for detecting and classifying potential point mutations in patient’s BRCA1 gene.Triple-negative breast cancers often contain inactivation of the DNA repair gene BRCA1. In fact, as much as 30% of breast cancers are thought to have some degree of BRCA1 inactivation. A key feature of the structure of a genes like BRCA1 is that their transcripts are typically subdivided into exon and intron regions. Exon regions are retained in the final mature mRNA molecule, while intron regions are cut out during post-transcriptional process. Indeed my focus was on exon’s mutations.
Point mutations
Genetic mutations can be classified based on their effect on the protein structure:
- Missense: these mutations change the coded amino acid, hence they influence the final protein structure. Their effect can be uncertain or pathogenic.
- Nonsense: these mutations cause a shift in the reading frame or the formation of a premature stop codon (truncated protein). Very often these mutations are pathogenic.
But there’s also another class of genetic mutations, rather infrequent, but possible:
- Start-loss: these mutations affect the initiation codon (the very first amino acid of the protein — a Methionine), and their effect on the final protein structure is very often pathogenic.
- Stop-loss: even rarer than start-loss, these mutations affect the last protein amino acid. Also for these mutations, the effect isn’t easy to understand
I want to try to clarify the strand issue. Consider the following stretch of double stranded DNA which encodes a short peptide:
The actual biological transcription process works from the template strand of the DNA. This is the reason why the reference and the patient sequences are indeed taken from reverse strand (aka Watson strand, strand −1).
The Data
There is used various types of data:
- Reference and patient DNA of BRCA1 gene and exon’s boundaries were taken from National Library of Medicine.
- All the DNA variants (mutation) were taken from Health University of Utah
Patient data
For the current research there is a patient DNA sequence with anomalies in exons: 2, 3, 4, 7, 13 and 14. Most of them are non pathogenic. Only one (in exon 2) is a pathogenic variant c.1A>G - mutation type:start loss.
“Start-loss: these mutations affect the initiation codon, i.e. the very first amino acid of the protein (which is a Methionine), and their effect on the final protein structure (and therefore on the individual’s clinical picture), is anything but easily deducible.” The quote was taken from Breda Genetics
Implementation
After wide research I decided to use a convolutional autoencoder for anomaly detection in patient DNA. For the anomaly detection part I transformed the DNA sequences into arrays of numbers (1–5).Each one corresponding to each nucleotide and 5 is for letter ’N’ (unknown) . After many experiments I decided to use different approach for the classification part. I transformed the DNA sequences into arrays of one hot encoded nucleotides. The latter seemed more robust and reliable.
Because of the relatively small size of described DNA variants (about 240), I had to find something different than regular classifier. I decide to use a siamese network with contrastive loss function for classifying the anomalies.
def vectorize(seq):
"""
Vectorize the the DNA sequence
"""
str_sequence = ''
arr_sequence = []
#Concatenate lines of strings
for line in seq[0]:
str_sequence += line
#Populate the vector with the nucleotides
for n in str_sequence:
arr_sequence.append(n)
return np.array(arr_sequence)# Find the length of the longest exon
EXON_MAX_LENGTH = 0
for i in range(len(exon_boundaries)):
start = exon_boundaries.loc[i].start
end = exon_boundaries.loc[i].end
diff = end - start + 1
if diff > EXON_MAX_LENGTH:
EXON_MAX_LENGTH = diff#Vectorize the referent and patient data
ref_vec = vectorize(reference_data)
pat_vec = vectorize(patient_data)def normalize(seq):
"""
Normalize the sequence
"""
return seq / np.max(seq)def nucleotide_to_num(sequence):
"""
Transforms DNA sequence to a sequence of numbers coresponding to certain nucleotide
"""
ref_sequence_to_num = []
nucleotide_to_num = {
'A': 1.,
'C': 2.,
'G': 3.,
'T': 4.,
'N': 5.
}
for nucleodtide in sequence:
ref_sequence_to_num.append(nucleotide_to_num[nucleodtide])
return np.array(normalize(ref_sequence_to_num))Exons boundaries
The exon boundaries were taken from National Center for Biotechnology Information. They are basically the start and end position of the exon according the order of Human Genome version called Assembly GRCh38.p13 - Cr17.
For the current project only the exon sequences will be used. The following function is extracting only exons and making their length even number. I needed this even lengths because of the nature of the autoencoder.
Convolutional reconstruction autoencoder
The dataset for the anomaly detection part of the project is actually the entire DNA sequence of the BRCA1 gene. The idea is to train the model with the referent (healthy) DNA sequence. After that the model should be able to detect potential changes in the patient DNA. Since the input data is a 3d array i decided to use convolutional reconstruction autoencoder model The model will take input of shape (batch_size, sequence_length, num_features) and return output of the same shape. In this case, sequence_length is 81070 (the length of BRCA1 gene) and num_features is 4 (which corresponds to one hot encoded nucleotide).
An autoencoder is a special type of neural network that is trained to copy its input to its output.It will first encode the sequence into a lower dimensional latent representation, then decodes the latent representation back to an sequence.
The autoencoder was trained to minimize reconstruction error with the referent DNA (normal) only, then I used it to reconstruct patient data. My hypothesis was that the mutated sequence will have higher reconstruction error.
For the encoder part were used two Conv1D layers. For the "bottleneck" I used a Flatten layout. And for decoder part were used two Conv1DTranspose layers and one Conv1D with output size 1. The Conv1DTranspose purpose is to increase the volume size back to the original array spatial dimensions.
A special class MutationDetector was created to capsulate the autoencoder. This way it was easier to use it multiple times for each exon separately.
KRNL_SIZE = 7
latent_dim = 20
class MutationDetector(Model):
"""
Capsulate an autoencoder in a class.
This way it can be used multiple times as separate object.
"""
def __init__(self, input_a):
super(MutationDetector, self).__init__()
#Encoder part ----------------------------------------------
inputs = Input(shape=(input_a, 1))
conv = Conv1D(
filters=8, kernel_size=KRNL_SIZE, strides=2, padding="same", activation = 'relu')(inputs)
conv = Conv1D(
filters=4, kernel_size=KRNL_SIZE, strides=1, padding="same", activation="relu")(conv)
"""
Storing the shape of the last convolutional layer
to use it in the decoder part
"""
vol_size = K.int_shape(conv)
"""
Flatten the convolutional output to a 1d lattent space
Pass it to a dense layer
"""
flat = Flatten()(conv)
latent = Dense(latent_dim)(flat)
#Decoder part ----------------------------------------------
dense = Dense(np.prod(vol_size[1:]))(latent)
reshape = Reshape((vol_size[1], vol_size[2]))(dense)
conv_trans = Conv1DTranspose(
filters=4, kernel_size=KRNL_SIZE, strides=1, padding="same", activation = 'relu')(reshape)
conv_trans = Conv1DTranspose(
filters=8, kernel_size=KRNL_SIZE, strides=2,padding="same", activation = 'relu')(conv_trans)
outputs = Conv1DTranspose(
filters=1, kernel_size=KRNL_SIZE, padding="same", activation='sigmoid')(conv_trans)
self.autoencoder = Model(inputs, outputs)
def call(self, x):
return self.autoencoder(x)
def show_summary(self):
#Returns the summary of the current object
return self.autoencoder.summary()The reference_set was used as both the input and the target since this is a reconstruction model.
Anomaly detection
Here , for every exon was created a train set ,a test set and an autoencoder. Every model was trained and it made it’s predictions separately from the others.
for exon_n in range(len(referent_exons)):
#Creating the exon names
exon_name = 'exon_0' + str(exon_n+1) if exon_n+1 < 10 else 'exon_' + str(exon_n+1)
#Converting the sequences into required shape
reference_set = create_dataset(referent_exons[exon_n])
reference_set = reference_set.reshape(COPY_COUNT, len(referent_exons[exon_n]), 1)
#Calling the autoencoder
autoencoder = MutationDetector(reference_set.shape[1])
#Compile the model using optimizer adam and loss: mean absolute error
autoencoder.compile(optimizer=Adam(learning_rate=0.001), loss="mae")
tr_title = 'Training for {}'.format(exon_name)
print( '==========='+ tr_title +'============')
#Fitting the model
history = autoencoder.fit(reference_set, reference_set,
epochs=30,
batch_size=16,
validation_split=0.2,
verbose=2)
print(autoencoder.show_summary())
#Plot training result
training_plot(history, tr_title)
"""
Predict the referent data to calculate the training loss
it is needed to determine the reconstruction loss
"""
predict = autoencoder.predict(reference_set)
# Get reconstruction loss threshold.
threshold = calculate_threshold(predict[0], reference_set[0])
print("Reconstruction error threshold:", threshold)
# Loading patient exon sequence
patient_set = patient_exons[exon_n].reshape(1, len(referent_exons[exon_n]), 1)
test_predict = autoencoder.predict(patient_set)
#Searching for anomalies
detect_anomalies(test_predict[0] , patient_set[0], threshold, exon_name)
print("======== End of processing {} ======== \n".format(exon_name))===========Training for exon_01============
_________________________________________________________________
Layer (type) Output Shape Param #
=================================================================
input_1 (InputLayer) [(None, 94, 1)] 0
_________________________________________________________________
conv1d (Conv1D) (None, 47, 8) 64
_________________________________________________________________
conv1d_1 (Conv1D) (None, 47, 4) 228
_________________________________________________________________
flatten (Flatten) (None, 188) 0
_________________________________________________________________
dense (Dense) (None, 20) 3780
_________________________________________________________________
dense_1 (Dense) (None, 188) 3948
_________________________________________________________________
reshape (Reshape) (None, 47, 4) 0
_________________________________________________________________
conv1d_transpose (Conv1DTran (None, 47, 4) 116
_________________________________________________________________
conv1d_transpose_1 (Conv1DTr (None, 94, 8) 232
_________________________________________________________________
conv1d_transpose_2 (Conv1DTr (None, 94, 1) 57
=================================================================
Total params: 8,425
Trainable params: 8,425
Non-trainable params: 0
_________________________________________________________________
None
Reconstruction error threshold: 0.04515981674194336
Anomaly detection - exon_01
Count of anomalous nucleotides: 0
======== End of processing exon_01 ========
===========Training for exon_02============
Model: "model_1"
_________________________________________________________________
Layer (type) Output Shape Param #
=================================================================
input_2 (InputLayer) [(None, 80, 1)] 0
_________________________________________________________________
conv1d_2 (Conv1D) (None, 40, 8) 64
_________________________________________________________________
conv1d_3 (Conv1D) (None, 40, 4) 228
_________________________________________________________________
flatten_1 (Flatten) (None, 160) 0
_________________________________________________________________
dense_2 (Dense) (None, 20) 3220
_________________________________________________________________
dense_3 (Dense) (None, 160) 3360
_________________________________________________________________
reshape_1 (Reshape) (None, 40, 4) 0
_________________________________________________________________
conv1d_transpose_3 (Conv1DTr (None, 40, 4) 116
_________________________________________________________________
conv1d_transpose_4 (Conv1DTr (None, 80, 8) 232
_________________________________________________________________
conv1d_transpose_5 (Conv1DTr (None, 80, 1) 57
=================================================================
Total params: 7,277
Trainable params: 7,277
Non-trainable params: 0
_________________________________________________________________
None
Reconstruction error threshold: 0.22216740436851978
Anomaly detection - exon_02
Count of anomalous nucleotides: 2
Position of anomalous nucleotide: (array([ 0, 29]),)
Test losses [0.50054242 0.2221675 ]
======== End of processing exon_02 ========
===========Training for exon_03============
Model: "model_2"
_________________________________________________________________
Layer (type) Output Shape Param #
=================================================================
input_3 (InputLayer) [(None, 54, 1)] 0
_________________________________________________________________
conv1d_4 (Conv1D) (None, 27, 8) 64
_________________________________________________________________
conv1d_5 (Conv1D) (None, 27, 4) 228
_________________________________________________________________
flatten_2 (Flatten) (None, 108) 0
_________________________________________________________________
dense_4 (Dense) (None, 20) 2180
_________________________________________________________________
dense_5 (Dense) (None, 108) 2268
_________________________________________________________________
reshape_2 (Reshape) (None, 27, 4) 0
_________________________________________________________________
conv1d_transpose_6 (Conv1DTr (None, 27, 4) 116
_________________________________________________________________
conv1d_transpose_7 (Conv1DTr (None, 54, 8) 232
_________________________________________________________________
conv1d_transpose_8 (Conv1DTr (None, 54, 1) 57
=================================================================
Total params: 5,145
Trainable params: 5,145
Non-trainable params: 0
_________________________________________________________________
None
Reconstruction error threshold: 0.24999821186065674
Anomaly detection - exon_03
Count of anomalous nucleotides: 1
Position of anomalous nucleotide: (array([52]),)
Test losses [0.25322703]
======== End of processing exon_03 ========
===========Training for exon_04============
Model: "model_3"
_________________________________________________________________
Layer (type) Output Shape Param #
=================================================================
input_4 (InputLayer) [(None, 78, 1)] 0
_________________________________________________________________
conv1d_6 (Conv1D) (None, 39, 8) 64
_________________________________________________________________
conv1d_7 (Conv1D) (None, 39, 4) 228
_________________________________________________________________
flatten_3 (Flatten) (None, 156) 0
_________________________________________________________________
dense_6 (Dense) (None, 20) 3140
_________________________________________________________________
dense_7 (Dense) (None, 156) 3276
_________________________________________________________________
reshape_3 (Reshape) (None, 39, 4) 0
_________________________________________________________________
conv1d_transpose_9 (Conv1DTr (None, 39, 4) 116
_________________________________________________________________
conv1d_transpose_10 (Conv1DT (None, 78, 8) 232
_________________________________________________________________
conv1d_transpose_11 (Conv1DT (None, 78, 1) 57
=================================================================
Total params: 7,113
Trainable params: 7,113
Non-trainable params: 0
_________________________________________________________________
None
Reconstruction error threshold: 0.49754321575164795
Anomaly detection - exon_04
Count of anomalous nucleotides: 1
Position of anomalous nucleotide: (array([19]),)
Test losses [0.7474227]
======== End of processing exon_04 ========
===========Training for exon_05============
Model: "model_4"
_________________________________________________________________
Layer (type) Output Shape Param #
=================================================================
input_5 (InputLayer) [(None, 90, 1)] 0
_________________________________________________________________
conv1d_8 (Conv1D) (None, 45, 8) 64
_________________________________________________________________
conv1d_9 (Conv1D) (None, 45, 4) 228
_________________________________________________________________
flatten_4 (Flatten) (None, 180) 0
_________________________________________________________________
dense_8 (Dense) (None, 20) 3620
_________________________________________________________________
dense_9 (Dense) (None, 180) 3780
_________________________________________________________________
reshape_4 (Reshape) (None, 45, 4) 0
_________________________________________________________________
conv1d_transpose_12 (Conv1DT (None, 45, 4) 116
_________________________________________________________________
conv1d_transpose_13 (Conv1DT (None, 90, 8) 232
_________________________________________________________________
conv1d_transpose_14 (Conv1DT (None, 90, 1) 57
=================================================================
Total params: 8,097
Trainable params: 8,097
Non-trainable params: 0
_________________________________________________________________
None
Reconstruction error threshold: 0.14869280457496642
Anomaly detection - exon_05
Count of anomalous nucleotides: 1
Position of anomalous nucleotide: (array([11]),)
Test losses [0.14869289]
======== End of processing exon_05 ========
===========Training for exon_06============
Model: "model_5"
_________________________________________________________________
Layer (type) Output Shape Param #
=================================================================
input_6 (InputLayer) [(None, 140, 1)] 0
_________________________________________________________________
conv1d_10 (Conv1D) (None, 70, 8) 64
_________________________________________________________________
conv1d_11 (Conv1D) (None, 70, 4) 228
_________________________________________________________________
flatten_5 (Flatten) (None, 280) 0
_________________________________________________________________
dense_10 (Dense) (None, 20) 5620
_________________________________________________________________
dense_11 (Dense) (None, 280) 5880
_________________________________________________________________
reshape_5 (Reshape) (None, 70, 4) 0
_________________________________________________________________
conv1d_transpose_15 (Conv1DT (None, 70, 4) 116
_________________________________________________________________
conv1d_transpose_16 (Conv1DT (None, 140, 8) 232
_________________________________________________________________
conv1d_transpose_17 (Conv1DT (None, 140, 1) 57
=================================================================
Total params: 12,197
Trainable params: 12,197
Non-trainable params: 0
_________________________________________________________________
None
Reconstruction error threshold: 0.2492716908454895
Anomaly detection - exon_06
Count of anomalous nucleotides: 0
======== End of processing exon_06 ========
===========Training for exon_07============
Model: "model_6"
_________________________________________________________________
Layer (type) Output Shape Param #
=================================================================
input_7 (InputLayer) [(None, 106, 1)] 0
_________________________________________________________________
conv1d_12 (Conv1D) (None, 53, 8) 64
_________________________________________________________________
conv1d_13 (Conv1D) (None, 53, 4) 228
_________________________________________________________________
flatten_6 (Flatten) (None, 212) 0
_________________________________________________________________
dense_12 (Dense) (None, 20) 4260
_________________________________________________________________
dense_13 (Dense) (None, 212) 4452
_________________________________________________________________
reshape_6 (Reshape) (None, 53, 4) 0
_________________________________________________________________
conv1d_transpose_18 (Conv1DT (None, 53, 4) 116
_________________________________________________________________
conv1d_transpose_19 (Conv1DT (None, 106, 8) 232
_________________________________________________________________
conv1d_transpose_20 (Conv1DT (None, 106, 1) 57
=================================================================
Total params: 9,409
Trainable params: 9,409
Non-trainable params: 0
_________________________________________________________________
None
Reconstruction error threshold: 0.2598261833190918
Anomaly detection - exon_07
Count of anomalous nucleotides: 1
Position of anomalous nucleotide: (array([27]),)
Test losses [0.49831474]
======== End of processing exon_07 ========
===========Training for exon_08============
Model: "model_7"
_________________________________________________________________
Layer (type) Output Shape Param #
=================================================================
input_8 (InputLayer) [(None, 78, 1)] 0
_________________________________________________________________
conv1d_14 (Conv1D) (None, 39, 8) 64
_________________________________________________________________
conv1d_15 (Conv1D) (None, 39, 4) 228
_________________________________________________________________
flatten_7 (Flatten) (None, 156) 0
_________________________________________________________________
dense_14 (Dense) (None, 20) 3140
_________________________________________________________________
dense_15 (Dense) (None, 156) 3276
_________________________________________________________________
reshape_7 (Reshape) (None, 39, 4) 0
_________________________________________________________________
conv1d_transpose_21 (Conv1DT (None, 39, 4) 116
_________________________________________________________________
conv1d_transpose_22 (Conv1DT (None, 78, 8) 232
_________________________________________________________________
conv1d_transpose_23 (Conv1DT (None, 78, 1) 57
=================================================================
Total params: 7,113
Trainable params: 7,113
Non-trainable params: 0
_________________________________________________________________
None
Reconstruction error threshold: 0.22729121446609502
Anomaly detection - exon_08
Count of anomalous nucleotides: 0
======== End of processing exon_08 ========
===========Training for exon_09============
Model: "model_8"
_________________________________________________________________
Layer (type) Output Shape Param #
=================================================================
input_9 (InputLayer) [(None, 90, 1)] 0
_________________________________________________________________
conv1d_16 (Conv1D) (None, 45, 8) 64
_________________________________________________________________
conv1d_17 (Conv1D) (None, 45, 4) 228
_________________________________________________________________
flatten_8 (Flatten) (None, 180) 0
_________________________________________________________________
dense_16 (Dense) (None, 20) 3620
_________________________________________________________________
dense_17 (Dense) (None, 180) 3780
_________________________________________________________________
reshape_8 (Reshape) (None, 45, 4) 0
_________________________________________________________________
conv1d_transpose_24 (Conv1DT (None, 45, 4) 116
_________________________________________________________________
conv1d_transpose_25 (Conv1DT (None, 90, 8) 232
_________________________________________________________________
conv1d_transpose_26 (Conv1DT (None, 90, 1) 57
=================================================================
Total params: 8,097
Trainable params: 8,097
Non-trainable params: 0
_________________________________________________________________
None
Reconstruction error threshold: 0.02076256275177002
Anomaly detection - exon_09
Count of anomalous nucleotides: 0
======== End of processing exon_09 ========
===========Training for exon_10============
Model: "model_9"
_________________________________________________________________
Layer (type) Output Shape Param #
=================================================================
input_10 (InputLayer) [(None, 172, 1)] 0
_________________________________________________________________
conv1d_18 (Conv1D) (None, 86, 8) 64
_________________________________________________________________
conv1d_19 (Conv1D) (None, 86, 4) 228
_________________________________________________________________
flatten_9 (Flatten) (None, 344) 0
_________________________________________________________________
dense_18 (Dense) (None, 20) 6900
_________________________________________________________________
dense_19 (Dense) (None, 344) 7224
_________________________________________________________________
reshape_9 (Reshape) (None, 86, 4) 0
_________________________________________________________________
conv1d_transpose_27 (Conv1DT (None, 86, 4) 116
_________________________________________________________________
conv1d_transpose_28 (Conv1DT (None, 172, 8) 232
_________________________________________________________________
conv1d_transpose_29 (Conv1DT (None, 172, 1) 57
=================================================================
Total params: 14,821
Trainable params: 14,821
Non-trainable params: 0
_________________________________________________________________
None
Reconstruction error threshold: 0.4937223792076111
Anomaly detection - exon_10
Count of anomalous nucleotides: 0
======== End of processing exon_10 ========
===========Training for exon_11============
Model: "model_10"
_________________________________________________________________
Layer (type) Output Shape Param #
=================================================================
input_11 (InputLayer) [(None, 128, 1)] 0
_________________________________________________________________
conv1d_20 (Conv1D) (None, 64, 8) 64
_________________________________________________________________
conv1d_21 (Conv1D) (None, 64, 4) 228
_________________________________________________________________
flatten_10 (Flatten) (None, 256) 0
_________________________________________________________________
dense_20 (Dense) (None, 20) 5140
_________________________________________________________________
dense_21 (Dense) (None, 256) 5376
_________________________________________________________________
reshape_10 (Reshape) (None, 64, 4) 0
_________________________________________________________________
conv1d_transpose_30 (Conv1DT (None, 64, 4) 116
_________________________________________________________________
conv1d_transpose_31 (Conv1DT (None, 128, 8) 232
_________________________________________________________________
conv1d_transpose_32 (Conv1DT (None, 128, 1) 57
=================================================================
Total params: 11,213
Trainable params: 11,213
Non-trainable params: 0
_________________________________________________________________
None
Reconstruction error threshold: 0.30231815576553345
Anomaly detection - exon_11
Count of anomalous nucleotides: 0
======== End of processing exon_11 ========
===========Training for exon_12============
Model: "model_11"
_________________________________________________________________
Layer (type) Output Shape Param #
=================================================================
input_12 (InputLayer) [(None, 192, 1)] 0
_________________________________________________________________
conv1d_22 (Conv1D) (None, 96, 8) 64
_________________________________________________________________
conv1d_23 (Conv1D) (None, 96, 4) 228
_________________________________________________________________
flatten_11 (Flatten) (None, 384) 0
_________________________________________________________________
dense_22 (Dense) (None, 20) 7700
_________________________________________________________________
dense_23 (Dense) (None, 384) 8064
_________________________________________________________________
reshape_11 (Reshape) (None, 96, 4) 0
_________________________________________________________________
conv1d_transpose_33 (Conv1DT (None, 96, 4) 116
_________________________________________________________________
conv1d_transpose_34 (Conv1DT (None, 192, 8) 232
_________________________________________________________________
conv1d_transpose_35 (Conv1DT (None, 192, 1) 57
=================================================================
Total params: 16,461
Trainable params: 16,461
Non-trainable params: 0
_________________________________________________________________
None
Reconstruction error threshold: 0.40868947505950926
Anomaly detection - exon_12
Count of anomalous nucleotides: 0
======== End of processing exon_12 ========
===========Training for exon_13============
Model: "model_12"
_________________________________________________________________
Layer (type) Output Shape Param #
=================================================================
input_13 (InputLayer) [(None, 312, 1)] 0
_________________________________________________________________
conv1d_24 (Conv1D) (None, 156, 8) 64
_________________________________________________________________
conv1d_25 (Conv1D) (None, 156, 4) 228
_________________________________________________________________
flatten_12 (Flatten) (None, 624) 0
_________________________________________________________________
dense_24 (Dense) (None, 20) 12500
_________________________________________________________________
dense_25 (Dense) (None, 624) 13104
_________________________________________________________________
reshape_12 (Reshape) (None, 156, 4) 0
_________________________________________________________________
conv1d_transpose_36 (Conv1DT (None, 156, 4) 116
_________________________________________________________________
conv1d_transpose_37 (Conv1DT (None, 312, 8) 232
_________________________________________________________________
conv1d_transpose_38 (Conv1DT (None, 312, 1) 57
=================================================================
Total params: 26,301
Trainable params: 26,301
Non-trainable params: 0
_________________________________________________________________
None
Reconstruction error threshold: 0.36978744268417363
Anomaly detection - exon_13
Count of anomalous nucleotides: 2
Position of anomalous nucleotide: (array([125, 289]),)
Test losses [0.74480762 0.36979764]
======== End of processing exon_13 ========
===========Training for exon_14============
Model: "model_13"
_________________________________________________________________
Layer (type) Output Shape Param #
=================================================================
input_14 (InputLayer) [(None, 88, 1)] 0
_________________________________________________________________
conv1d_26 (Conv1D) (None, 44, 8) 64
_________________________________________________________________
conv1d_27 (Conv1D) (None, 44, 4) 228
_________________________________________________________________
flatten_13 (Flatten) (None, 176) 0
_________________________________________________________________
dense_26 (Dense) (None, 20) 3540
_________________________________________________________________
dense_27 (Dense) (None, 176) 3696
_________________________________________________________________
reshape_13 (Reshape) (None, 44, 4) 0
_________________________________________________________________
conv1d_transpose_39 (Conv1DT (None, 44, 4) 116
_________________________________________________________________
conv1d_transpose_40 (Conv1DT (None, 88, 8) 232
_________________________________________________________________
conv1d_transpose_41 (Conv1DT (None, 88, 1) 57
=================================================================
Total params: 7,933
Trainable params: 7,933
Non-trainable params: 0
_________________________________________________________________
None
Reconstruction error threshold: 0.2496626377105713
Anomaly detection - exon_14
Count of anomalous nucleotides: 1
Position of anomalous nucleotide: (array([4]),)
Test losses [0.49999022]
======== End of processing exon_14 ========
===========Training for exon_15============
Model: "model_14"
_________________________________________________________________
Layer (type) Output Shape Param #
=================================================================
input_15 (InputLayer) [(None, 78, 1)] 0
_________________________________________________________________
conv1d_28 (Conv1D) (None, 39, 8) 64
_________________________________________________________________
conv1d_29 (Conv1D) (None, 39, 4) 228
_________________________________________________________________
flatten_14 (Flatten) (None, 156) 0
_________________________________________________________________
dense_28 (Dense) (None, 20) 3140
_________________________________________________________________
dense_29 (Dense) (None, 156) 3276
_________________________________________________________________
reshape_14 (Reshape) (None, 39, 4) 0
_________________________________________________________________
conv1d_transpose_42 (Conv1DT (None, 39, 4) 116
_________________________________________________________________
conv1d_transpose_43 (Conv1DT (None, 78, 8) 232
_________________________________________________________________
conv1d_transpose_44 (Conv1DT (None, 78, 1) 57
=================================================================
Total params: 7,113
Trainable params: 7,113
Non-trainable params: 0
_________________________________________________________________
None
Reconstruction error threshold: 0.49994945526123047
Anomaly detection - exon_15
Count of anomalous nucleotides: 0
======== End of processing exon_15 ========
===========Training for exon_16============
Model: "model_15"
_________________________________________________________________
Layer (type) Output Shape Param #
=================================================================
input_16 (InputLayer) [(None, 42, 1)] 0
_________________________________________________________________
conv1d_30 (Conv1D) (None, 21, 8) 64
_________________________________________________________________
conv1d_31 (Conv1D) (None, 21, 4) 228
_________________________________________________________________
flatten_15 (Flatten) (None, 84) 0
_________________________________________________________________
dense_30 (Dense) (None, 20) 1700
_________________________________________________________________
dense_31 (Dense) (None, 84) 1764
_________________________________________________________________
reshape_15 (Reshape) (None, 21, 4) 0
_________________________________________________________________
conv1d_transpose_45 (Conv1DT (None, 21, 4) 116
_________________________________________________________________
conv1d_transpose_46 (Conv1DT (None, 42, 8) 232
_________________________________________________________________
conv1d_transpose_47 (Conv1DT (None, 42, 1) 57
=================================================================
Total params: 4,161
Trainable params: 4,161
Non-trainable params: 0
_________________________________________________________________
None
Reconstruction error threshold: 0.10628277063369751
Anomaly detection - exon_16
Count of anomalous nucleotides: 0
======== End of processing exon_16 ========
===========Training for exon_17============
Model: "model_16"
_________________________________________________________________
Layer (type) Output Shape Param #
=================================================================
input_17 (InputLayer) [(None, 84, 1)] 0
_________________________________________________________________
conv1d_32 (Conv1D) (None, 42, 8) 64
_________________________________________________________________
conv1d_33 (Conv1D) (None, 42, 4) 228
_________________________________________________________________
flatten_16 (Flatten) (None, 168) 0
_________________________________________________________________
dense_32 (Dense) (None, 20) 3380
_________________________________________________________________
dense_33 (Dense) (None, 168) 3528
_________________________________________________________________
reshape_16 (Reshape) (None, 42, 4) 0
_________________________________________________________________
conv1d_transpose_48 (Conv1DT (None, 42, 4) 116
_________________________________________________________________
conv1d_transpose_49 (Conv1DT (None, 84, 8) 232
_________________________________________________________________
conv1d_transpose_50 (Conv1DT (None, 84, 1) 57
=================================================================
Total params: 7,605
Trainable params: 7,605
Non-trainable params: 0
_________________________________________________________________
None
Reconstruction error threshold: 0.7493552491650917
Anomaly detection - exon_17
Count of anomalous nucleotides: 0
======== End of processing exon_17 ========"Mutations of BRCA1 gene were found in: {}".format(anomalous_exon_names)"Mutations of BRCA1 gene were found in: ['exon_02', 'exon_03', 'exon_04', 'exon_05', 'exon_07', 'exon_13', 'exon_14']"All the anomalous exons were found by the autoencoder. The problem is that there is one additional — exon_5. It is not in the list of known anomalous exons of the patient DNA. The exons of interest were found and now i had to classify them.
Classifying the the anomalous sequences
Siamese network
During my research I realized that there aren’t so much DNA variants for different exons in BRCA1 gene. The overall count of the variants is 241.
So I came to the conclusion that I need something different than common classifier to determine which anomalous sequence belongs to which class (pathogenic , not pathogenic) , mutation type or nucleotide change.
Then I decide to try siamese network with contrastive loss function.
Practical, real-world use cases of siamese networks include face recognition, signature verification, prescription pill identification, and many more!
Furthermore, siamese networks can be trained with very little data, which actually I was aiming for.
The concept of a siamese network:
Two convolutional network to be merged and to output the similarity of two entries in a pair.
DNA Variants
The DNA variants were taken from BRCA1 database of Health University of Utah. The data includes all the recorded mutations classified as Definitely pathogenic , likely not pathogenic and not pathogenic.All the nucleotide changes have mutation type: Start loss,Splice site, Missense or Nonsence
variants_data.head()
Only the variants of the anomalous exons were taken for the classification. Next steps were to vectorize the DNA sequences of the variants. Because the network requires to the sequences to be with equal lengths, the shorter sequences were made as long as the longest sequence by adding sufficient count of ‘N’s at their ends.
Missing information for variants
Unfortunately, during my research I couldn’t find any variants for exons: 1, 8 and 10. So this is the reason why I didn’t considered them for the classification.
# Vectorize the patient BRCA1 gene sequence
vec_patient_sequence = vectorize(patient_data)
patient_exons = {}
exon_boundaries_df = exon_boundaries.set_index('name')
missing_data = ['exon_01','exon_08','exon_10']
# Create dictionary with detetected by the autencoder anomalies
for name in anomalous_exon_names:
if name in missing_data:
continue
address = exon_boundaries_df.loc[name]
patient_exons[name] = vec_patient_sequence[address.start : address.end]def nucleotide_to_vector(sequence):
"""
Transforms DNA sequence to a array of one-hot encoded nucleotides
"""
ref_sequence_to_vec = []
nucleotide_to_vec = {
'A': [1., 0., 0., 0.],
'C': [0., 1., 0., 0.],
'G': [0. ,0. ,1., 0.],
'T': [0., 0., 0., 1.],
'N': [0., 0., 0., 0.]
}
for nucleodtide in sequence:
ref_sequence_to_vec.append(nucleotide_to_vec[nucleodtide])
return np.array(ref_sequence_to_vec)Creating the variants data frame with vectorized sequences, so that all the variants have the same length. The shorter sequences received subsequence with letter ’N’ at the end.
EXON_MAX_LENGTH = len(max(variants_data.variant, key=len))
vec_variants_dict = {'variant':[],'location':[],'classification':[],'nucleotide change':[],
'mutation type':[]}
for i in range(len(variants_data)):
row = variants_data.loc[i]
sequence = row.variant
# Make the exons as long as the longest one
sequence = sequence + 'N'* (EXON_MAX_LENGTH-len(sequence))
#Transform the nucleotides sequence to one hot vector
sequence = nucleotide_to_vector(sequence)
#Populate the variant dictionary
vec_variants_dict['variant'].append(sequence)
vec_variants_dict['location'].append(row.location.strip())
vec_variants_dict['classification'].append(row.classification.strip())
vec_variants_dict['nucleotide change'].append(row[' nucleotide change'].strip())
vec_variants_dict['mutation type'].append(row['mutation type'])
#Populate a dataframe with required data
vec_variants_df = pd.DataFrame(vec_variants_dict).sample(frac=1)
vec_variants_df = vec_variants_df.reset_index()
vec_variants_df = vec_variants_df.drop(columns=['index'])vec_variants_df.head()
"""
Create and populate a dictionary with data, so that after the predictions,
the data of predicted variant should be acessible
"""
pairs_dict={'current':[], 'candidate':[], 'current nuc change':[],
'current class':[],'paired class':[], 'paired nuc change':[], 'class':[]}
def populate_pairs_data(df, row, is_similar):
"""
Populate a dictionary with pairs of sequences
"""
current_seq = row.variant
current_nuc_change = row['nucleotide change']
for j in range(len(df)):
c_row = df.loc[j]
pairs_dict['current'].append(current_seq)
pairs_dict['candidate'].append(c_row.variant)
sim_class = np.array([1.]) if is_similar else np.array([0.])
pairs_dict['class'].append(sim_class)
pairs_dict['current nuc change'].append(current_nuc_change)
pairs_dict['paired nuc change'].append(c_row['nucleotide change'])
pairs_dict['paired class'].append(c_row.classification)
pairs_dict['current class'].append(row.classification)Since there are two subnetworks, there must be two inputs to the model. When training siamese networks I need to have positive pairs and negative pairs:
- Positive pairs: Two sequences that belong to the same class (pathogenic — pathogenic).
- Negative pairs: Two sequences that belong to different classes (non pathogenic — pathogenic).
Next steps were to create such dataset that can fit the requirements of the siamese network inputs.
TRAIN_SIZE = 21000
pairs_df = pd.DataFrame()
def generate_pairs(df):
"""
initialize two empty lists to hold the (sequence, sequence) pairs and
labels to indicate if a pair is positive or negative
"""
global pairs_df
pair_sequences = []
pair_classes = []
non_pair_columns = ['class','current nuc change', 'paired nuc change',
'paired class','current class']
# loop over all dataset
for i in range(len(df)):
row = df.loc[i]
# take the current sequence and it's class
current_seq = row.variant
current_class = row.classification
current_location = row.location
#Create a Dataframe with variants with same classes of the current exon
current_similars_df=df[(df.location == current_location)
& (df.classification == current_class)].reset_index()
populate_pairs_data(current_similars_df, row, True)
#Create a Dataframe with variants with different classes of the current exon
current_non_similar_df=df[(df.classification != current_class)].reset_index()
populate_pairs_data(current_non_similar_df, row, False)
#Populate a dataframe and shuffle
pairs_df = pd.DataFrame(pairs_dict).sample(frac=1)
pairs_df = pairs_df.reset_index()
pairs_df = pairs_df.drop(columns=['index'])
#Construct and return the pairs
pair_sequences = np.array(pairs_df.drop(columns=non_pair_columns, axis=1).values.tolist()).astype(np.float)
pair_classes = np.array(pairs_df['class'].values.tolist()).astype(np.float)
return (pair_sequences[:TRAIN_SIZE],pair_sequences[TRAIN_SIZE:],
pair_classes[:TRAIN_SIZE],pair_classes[TRAIN_SIZE:])#Defining train and test sets for the siamese network
train_pairs, test_pairs, train_classes, test_classes = generate_pairs(vec_variants_df)print("Train pairs shape:{} - Train labels shape:{}".format(train_pairs.shape,train_classes.shape))Train pairs shape:(21000, 2, 312, 4) - Train labels shape:(21000, 1)print("Test pairs shape:{} - Test labels shape:{}".format(test_pairs.shape,test_classes.shape))Test pairs shape:(7525, 2, 312, 4) - Test labels shape:(7525, 1)So the train data has 21000 pairs with shape 312 nucleotides and last dimension is representing the encoded nucleotide itself.
seq_input_shape = [train_pairs.shape[2],train_pairs.shape[3]]The twin prototype
I constructed the prototype of the twins model, defining three sets of Conv1D layer with relu activation. Each convolutional layer has a total of 64 filters with size 7 following by MaxPooling1D with size = 2.
GlobalAveragePooling1D - This layer performs exactly the same operation as the 1D Average pooling layer, except that the pool size is the size of the entire input of the layer,it computes a single average value for each of the input channels (the second dimension). A fully-connected layer was defined with the specified size = 48.
Finally I normalized the features using L2 normalization before using Contrastive Loss.
def build_twin_model(input_shape):
inputs = Input(input_shape)
layer = Conv1D(filters=64, kernel_size=7, padding="same", activation="relu")(inputs)
layer = MaxPooling1D(pool_size=2)(layer)
layer = Conv1D(filters=64, kernel_size=7, padding="same", activation="relu")(layer)
layer = MaxPooling1D(pool_size=2)(layer)
layer = Conv1D(filters=64, kernel_size=7, padding="same", activation="relu")(layer)
layer = MaxPooling1D(pool_size=2)(layer)
pooled = GlobalAveragePooling1D()(layer)
dense = Dense(48)(pooled)
#Normalizing the output distances
outputs = Lambda(lambda x: K.l2_normalize(x,axis=1))(dense)
twin = Model(inputs, outputs)
return twin#Create the separate inputs for the twin nets
pos_input = Input(seq_input_shape)
neg_input = Input(seq_input_shape)
#Construct the twins
prototype_twin = build_twin_model(seq_input_shape)
pos_twin = prototype_twin(pos_input)
neg_twin = prototype_twin(neg_input)Here a Lambda layer was used to compute the euclidean distance between the outputs of the twin networks. Eventually they became the output of the sieamese network through a Dence layer.
def euclidean_distance(vectors):
"""Calculating the euclidian distance between twin ouputs"""
# unpack the vectors into separate lists
(vec_A, vec_B) = vectors
# compute the sum of squared distances between the vectors
squared = K.sum(K.square(vec_A - vec_B), axis=1, keepdims=True)
# Return the distances
return K.sqrt(K.maximum(squared, K.epsilon()))#construct the siamese network
distance = Lambda(euclidean_distance)([pos_twin, neg_twin])
outputs = Dense(1, activation="sigmoid")(distance)
siamese = Model(inputs=[pos_input, neg_input], outputs=outputs)Contrastive loss
For the current task binary cross entropy can be valid loss function. Тhe goal of a siamese network isn’t to classify a set of pairs but instead to differentiate between them. Essentially, contrastive loss is evaluating how good the siamese network is distinguishing between the pairs. There is distance based loss function called contrastive loss.
def contrastive_loss(true_labels, dist, margin=1):
"""
Calculate the contrastive loss between the true labels and
the predicted distances
"""
squared_dist = K.square(dist)
squared_margin = K.square(K.maximum(margin - dist, 0))
return K.mean(1/2 * true_labels * squared_dist + 1/2 * (1 - true_labels) * squared_margin)siamese.compile(loss=contrastive_loss, optimizer=Adam(learning_rate=0.005))
siamese.summary()Model: "model_18"
____________________________________________________________________
Layer (type) Output Shape Param # Connected to
====================================================================
input_18 (InputLayer) [(None, 312, 4)] 0
____________________________________________________________________
input_19 (InputLayer) [(None, 312, 4)] 0
____________________________________________________________________
model_17 (Functional) (None, 48) 62448 input_18[0][0]
input_19[0][0]
____________________________________________________________________
lambda_1 (Lambda) (None, 1) 0 model_17[0][0]
model_17[1][0]
____________________________________________________________________
dense_35 (Dense) (None, 1) 2 lambda_1[0][0]
====================================================================
Total params: 62,450
Trainable params: 62,450
Non-trainable params: 0
____________________________________________________________________# Fit the siamese network
history = siamese.fit(
[train_pairs[:,0], train_pairs[:,1]], [train_classes],
validation_data=([test_pairs[:,0], test_pairs[:,1]], [test_classes]),
batch_size=128,
epochs=25)Epoch 1/25
165/165 [==============================] - 20s 81ms/step - loss: 0.0740 - val_loss: 0.0338
Epoch 2/25
165/165 [==============================] - 4s 22ms/step - loss: 0.0314 - val_loss: 0.0279
Epoch 3/25
165/165 [==============================] - 4s 22ms/step - loss: 0.0276 - val_loss: 0.0254
Epoch 4/25
165/165 [==============================] - 4s 22ms/step - loss: 0.0245 - val_loss: 0.0233
Epoch 5/25
165/165 [==============================] - 4s 21ms/step - loss: 0.0225 - val_loss: 0.0203
Epoch 6/25
165/165 [==============================] - 3s 21ms/step - loss: 0.0196 - val_loss: 0.0191
Epoch 7/25
165/165 [==============================] - 3s 21ms/step - loss: 0.0178 - val_loss: 0.0180
Epoch 8/25
165/165 [==============================] - 4s 21ms/step - loss: 0.0157 - val_loss: 0.0121
Epoch 9/25
165/165 [==============================] - 3s 21ms/step - loss: 0.0124 - val_loss: 0.0102
Epoch 10/25
165/165 [==============================] - 3s 21ms/step - loss: 0.0094 - val_loss: 0.0097
Epoch 11/25
165/165 [==============================] - 3s 21ms/step - loss: 0.0086 - val_loss: 0.0095
Epoch 12/25
165/165 [==============================] - 4s 22ms/step - loss: 0.0067 - val_loss: 0.0061
Epoch 13/25
165/165 [==============================] - 4s 22ms/step - loss: 0.0060 - val_loss: 0.0046
Epoch 14/25
165/165 [==============================] - 4s 21ms/step - loss: 0.0049 - val_loss: 0.0039
Epoch 15/25
165/165 [==============================] - 4s 22ms/step - loss: 0.0037 - val_loss: 0.0042
Epoch 16/25
165/165 [==============================] - 4s 22ms/step - loss: 0.0034 - val_loss: 0.0030
Epoch 17/25
165/165 [==============================] - 4s 22ms/step - loss: 0.0026 - val_loss: 0.0021
Epoch 18/25
165/165 [==============================] - 4s 22ms/step - loss: 0.0021 - val_loss: 0.0017
Epoch 19/25
165/165 [==============================] - 4s 22ms/step - loss: 0.0019 - val_loss: 0.0016
Epoch 20/25
165/165 [==============================] - 4s 22ms/step - loss: 0.0016 - val_loss: 0.0014
Epoch 21/25
165/165 [==============================] - 4s 22ms/step - loss: 0.0015 - val_loss: 0.0013
Epoch 22/25
165/165 [==============================] - 4s 22ms/step - loss: 0.0015 - val_loss: 0.0012
Epoch 23/25
165/165 [==============================] - 4s 23ms/step - loss: 0.0012 - val_loss: 0.0011
Epoch 24/25
165/165 [==============================] - 4s 21ms/step - loss: 0.0014 - val_loss: 0.0011
Epoch 25/25
165/165 [==============================] - 4s 22ms/step - loss: 0.0019 - val_loss: 0.0021#Plot siamese loss functions
plt.figure()
plt.plot(history.history['loss'])
plt.plot(history.history['val_loss'])
plt.title('model loss')
plt.ylabel('loss')
plt.xlabel('epoch')
plt.legend(['train', 'validation'])
plt.show()
predict = siamese.predict([test_pairs[:, 0], test_pairs[:, 1]])Here I have created a data frame with the ground truth labels and the predicted distances between the pairs. We have to keep in mind that the smallest the distance is — the higher is the probability the sequences to belong of the same class. So in fact the prediction of 0.9 should correspond to label 0 (opposite classes) and vice versa.
#Create a dataframe with predicted distances and ground truth labels
pr_data = {'predicted':predict.reshape( -1).round(2), 'actual':test_classes.reshape(-1)}
comparisson_df = pd.DataFrame(data = pr_data).astype('float32')comparisson_df.head()
#Sorting the predicted distances
sorted_distances = np.sort(comparisson_df.predicted.unique())# Getting the train part of the combined_df dataframe
combined_df = pairs_df[TRAIN_SIZE:]
combined_df = combined_df.reset_index()
#Join the columns with predictions
combined_df = combined_df.join(comparisson_df.predicted)
#Drop not needed columns
combined_df = combined_df.drop(columns=['current', 'candidate'])Here I have constructed a data frame which helped me to find the related data to the predicted variants.
combined_df.head()
The margin
Checking the point from which we can say that a pair similar or not. Here I looped through the distances and calculated the count of false positive and false negative pairs. Then I found the distance where false negative and false positive counts were minimal. In this case the distance was 0.58.
counts = {'dist':[],'false positive':[] ,'false negative':[]}
for d in sorted_distances:
# Find count of false positive and false negative counts for the current distance
false_positive_count = len(combined_df[(combined_df['predicted'] > d) & (combined_df['class'] == 1)])
false_negaitive_count = (len(combined_df[(combined_df['predicted'] < d) & (combined_df['class'] == 0)]))
#Populate a dictionary
counts['dist'].append(d)
counts['false positive'].append(false_positive_count)
counts['false negative'].append(false_negaitive_count)
print('For distance: {} -> False positive count: {} , False negative count: {} '.format(d,false_positive_count,
false_negaitive_count))
print('False positive median: {} - False negative median: {}'.format(np.median(counts['false positive']),
np.median(counts['false negative'])))For distance: 0.07999999821186066 -> False positive count: 556 , False negative count: 0
For distance: 0.09000000357627869 -> False positive count: 397 , False negative count: 0
For distance: 0.10000000149011612 -> False positive count: 259 , False negative count: 0
For distance: 0.10999999940395355 -> False positive count: 184 , False negative count: 0
For distance: 0.11999999731779099 -> False positive count: 146 , False negative count: 0
For distance: 0.12999999523162842 -> False positive count: 118 , False negative count: 0
For distance: 0.14000000059604645 -> False positive count: 97 , False negative count: 0
For distance: 0.15000000596046448 -> False positive count: 77 , False negative count: 0
For distance: 0.1599999964237213 -> False positive count: 60 , False negative count: 2
For distance: 0.17000000178813934 -> False positive count: 51 , False negative count: 3
For distance: 0.18000000715255737 -> False positive count: 45 , False negative count: 3
For distance: 0.1899999976158142 -> False positive count: 41 , False negative count: 5
For distance: 0.20000000298023224 -> False positive count: 39 , False negative count: 5
For distance: 0.20999999344348907 -> False positive count: 37 , False negative count: 6
For distance: 0.2199999988079071 -> False positive count: 35 , False negative count: 7
For distance: 0.23000000417232513 -> False positive count: 34 , False negative count: 7
For distance: 0.23999999463558197 -> False positive count: 30 , False negative count: 7
For distance: 0.25 -> False positive count: 27 , False negative count: 7
For distance: 0.28999999165534973 -> False positive count: 25 , False negative count: 8
For distance: 0.3100000023841858 -> False positive count: 24 , False negative count: 8
For distance: 0.3199999928474426 -> False positive count: 22 , False negative count: 8
For distance: 0.3400000035762787 -> False positive count: 22 , False negative count: 8
For distance: 0.3700000047683716 -> False positive count: 22 , False negative count: 9
For distance: 0.3799999952316284 -> False positive count: 21 , False negative count: 10
For distance: 0.38999998569488525 -> False positive count: 20 , False negative count: 12
For distance: 0.4099999964237213 -> False positive count: 20 , False negative count: 12
For distance: 0.41999998688697815 -> False positive count: 19 , False negative count: 14
For distance: 0.4300000071525574 -> False positive count: 19 , False negative count: 14
For distance: 0.44999998807907104 -> False positive count: 17 , False negative count: 16
For distance: 0.46000000834465027 -> False positive count: 15 , False negative count: 18
For distance: 0.47999998927116394 -> False positive count: 14 , False negative count: 18
For distance: 0.5 -> False positive count: 14 , False negative count: 18
For distance: 0.5099999904632568 -> False positive count: 13 , False negative count: 19
For distance: 0.5299999713897705 -> False positive count: 12 , False negative count: 19
For distance: 0.5400000214576721 -> False positive count: 10 , False negative count: 20
For distance: 0.550000011920929 -> False positive count: 7 , False negative count: 23
For distance: 0.5600000023841858 -> False positive count: 6 , False negative count: 25
For distance: 0.5799999833106995 -> False positive count: 5 , False negative count: 25
For distance: 0.6000000238418579 -> False positive count: 5 , False negative count: 25
For distance: 0.6100000143051147 -> False positive count: 5 , False negative count: 27
For distance: 0.6200000047683716 -> False positive count: 5 , False negative count: 29
For distance: 0.6399999856948853 -> False positive count: 5 , False negative count: 30
For distance: 0.6600000262260437 -> False positive count: 5 , False negative count: 31
For distance: 0.6700000166893005 -> False positive count: 4 , False negative count: 32
For distance: 0.6800000071525574 -> False positive count: 4 , False negative count: 32
For distance: 0.6899999976158142 -> False positive count: 2 , False negative count: 33
For distance: 0.699999988079071 -> False positive count: 1 , False negative count: 33
For distance: 0.7099999785423279 -> False positive count: 1 , False negative count: 35
For distance: 0.7200000286102295 -> False positive count: 1 , False negative count: 37
For distance: 0.7300000190734863 -> False positive count: 1 , False negative count: 38
For distance: 0.7400000095367432 -> False positive count: 1 , False negative count: 39
For distance: 0.75 -> False positive count: 1 , False negative count: 40
For distance: 0.7599999904632568 -> False positive count: 1 , False negative count: 42
For distance: 0.7699999809265137 -> False positive count: 1 , False negative count: 43
For distance: 0.7799999713897705 -> False positive count: 0 , False negative count: 44
For distance: 0.7900000214576721 -> False positive count: 0 , False negative count: 47
For distance: 0.8100000023841858 -> False positive count: 0 , False negative count: 49
For distance: 0.8199999928474426 -> False positive count: 0 , False negative count: 56
For distance: 0.8299999833106995 -> False positive count: 0 , False negative count: 58
For distance: 0.8399999737739563 -> False positive count: 0 , False negative count: 62
For distance: 0.8500000238418579 -> False positive count: 0 , False negative count: 64
For distance: 0.8700000047683716 -> False positive count: 0 , False negative count: 67
For distance: 0.8799999952316284 -> False positive count: 0 , False negative count: 69
For distance: 0.8899999856948853 -> False positive count: 0 , False negative count: 72
For distance: 0.8999999761581421 -> False positive count: 0 , False negative count: 80
For distance: 0.9100000262260437 -> False positive count: 0 , False negative count: 84
For distance: 0.9200000166893005 -> False positive count: 0 , False negative count: 101
For distance: 0.9300000071525574 -> False positive count: 0 , False negative count: 118
For distance: 0.9399999976158142 -> False positive count: 0 , False negative count: 145
For distance: 0.949999988079071 -> False positive count: 0 , False negative count: 254
For distance: 0.9599999785423279 -> False positive count: 0 , False negative count: 545
For distance: 0.9700000286102295 -> False positive count: 0 , False negative count: 867
For distance: 0.9800000190734863 -> False positive count: 0 , False negative count: 1668
For distance: 0.9900000095367432 -> False positive count: 0 , False negative count: 2682
For distance: 1.0 -> False positive count: 0 , False negative count: 5895
False positive median: 5.0 - False negative median: 25.0counts_df = pd.DataFrame(counts)
best_dist = counts_df.loc[(counts_df.dist > 0.57) & (counts_df.dist < 0.59)]
best_dist
# compute final accuracy on test sets
faults_count = best_dist['false positive'] + best_dist['false negative']
'* Accuracy : %0.2f%%' % (100 - faults_count / (len(test_classes) / 100))'* Accuracy : 99.60%'The results
It is time to get the final results. Pairs between the variants and anomalous sequences were created and passes to the model.
def vec_to_sequence(seq):
"""
Converting the vector back to a string sequence
"""
sequence = ''
for c in range(len(seq)):
sequence+=seq[c]
return sequenceConstructing an array with pairs of anomalous patient data and variants data and make predictions
for exon_name in patient_exons.keys():
x = None
# Converting the vector back to a string sequence
sequence = vec_to_sequence(patient_exons[exon_name])
# Making the patient exon with the same length as the reference
anchor = sequence + 'N'* (EXON_MAX_LENGTH-len(sequence))
# Converting the sequence to one hot vector
anchor = nucleotide_to_vector(anchor)
# Construct a dataframe with variants of the current exon
anchor_variants = vec_variants_df[vec_variants_df.location == exon_name].reset_index()
#Construct pairs for the predictions
for v in range(len(anchor_variants)):
variant = anchor_variants.loc[v].variant
if x is None:
x = [[anchor], [variant]]
else:
x[0].append(anchor)
x[1].append(variant)
x = np.array(x)
"""
Make a prediction for the current annomaly
Get the index of the prediction with the highest probabiliry
Get the result from the variants dataframe
"""
index = np.argmin(siamese.predict([x[0], x[1]]), axis=0)
predicted_variant = anchor_variants.loc[index].values
#Printing the result for curen exon
print("Result for {} - mutation type: {}, nucleotide change: {} class: {}"
.format(exon_name,predicted_variant[0][5],
predicted_variant[0][4],
predicted_variant[0][3]))Result for exon_02 - mutation type: start loss, nucleotide change: c.1A>G class: pathogenic
Result for exon_03 - mutation type: missense, nucleotide change: c.133A>C class: non pathogenic
Result for exon_04 - mutation type: missense, nucleotide change: c.154C>A class: non pathogenic
Result for exon_05 - mutation type: -, nucleotide change: - class: non pathogenic
Result for exon_07 - mutation type: misssense, nucleotide change: c.469T>C class: non pathogenic
Result for exon_13 - mutation type: misssense, nucleotide change: c.4402A>C class: non pathogenic
Result for exon_14 - mutation type: misssense, nucleotide change: c.4520G>C class: non pathogenicThe results were quite accurate. Despite the anomaly in exon_5 all the other variants were exact.
Conclusion and Future work
DNA sequencing is a complex, multi-step process that is prone to errors and anomalies. The goal of this project was to identify such anomalies (not every anomaly is pathogenic) and classify them according to existing data base. During this vast research I found that the deep neural networks (particular autoencoders and siamese networks) performed quite well despite a massive data imbalance.
This project was focused only on point mutations. This means that the length of the patient DNA sequence is as long as the referent one. With more time and computational resources, there can be develop more robust project, which can detect and classify the other type of mutation of the DNA (deletions, additions , etc.)
Reference
메타데이터
- post_id
- 84c82b8b400a
- slug
- dna-cnn-autoencoder-and-siamese-network-84c82b8b400a
- url
- https://medium.com/@boyko-zhelev/dna-cnn-autoencoder-and-siamese-network-84c82b8b400a
- canonical_url
- https://medium.com/@boyko-zhelev/dna-cnn-autoencoder-and-siamese-network-84c82b8b400a
- author_url
- https://medium.com/@boyko-zhelev
- status
- ok
- fetched_at
- 2026-07-28 02:01:45
