<jupyter_start><jupyter_text>DeepVCF: A Variant Caller using deep learning  This tutorial on how to train and use DeepVCF to call variants and produce a Variant Call Format (VCF) file<jupyter_code>from pprint import pprint

import pandas as pd
import seaborn as sns

from DeepVCF import pathing
from DeepVCF import biopandas as bpd
from DeepVCF import alignment
from DeepVCF.core import DeepVCF

deepvcf = DeepVCF()<jupyter_output>=== No GPU Detected ===<jupyter_text>Pulling in Training & Testing Data We created these through the read simulator DWGSIM in the [tutorial](./creating-data-for-usage-demo.ipynb). The training and testing will be for both 2% platform base error (low) rate and 10% platform base error (high) rate with a .5% mutation rate for each. The goal is compare DeepVCF againset BCFtools regarding SNP variant calling when given low and high platform error rates.<jupyter_code>ref_train_be10percent = pathing('./example_files/staphylococcus/train/MRSA107.fasta')
bam_train_be10percent = pathing('./example_files/staphylococcus/train/train-base-error-10percent.mapped.bam')
vcf_train_be10percent = pathing('./example_files/staphylococcus/train/train-base-error-10percent.mutations.vcf')

ref_test_be10percent = pathing('./example_files/staphylococcus/test/GCF_000013425.1_ASM1342v1_genomic_reference.fasta')
bam_test_be10percent = pathing('./example_files/staphylococcus/test/test-base-error-10percent.mapped.bam')
vcf_test_be10percent = pathing('./example_files/staphylococcus/test/test-base-error-10percent.mutations.vcf')

ref_train_be2percent = pathing('./example_files/staphylococcus/train/MRSA107.fasta')
bam_train_be2percent = pathing('./example_files/staphylococcus/train/train-base-error-2percent.mapped.bam')
vcf_train_be2percent = pathing('./example_files/staphylococcus/train/train-base-error-2percent.mutations.vcf')

ref_test_be2percent = pathing('./example_files/staphylococcus/test/GCF_000013425.1_ASM1342v1_genomic_reference.fasta')
bam_test_be2percent = pathing('./example_files/staphylococcus/test/test-base-error-2percent.mapped.bam')
vcf_test_be2percent = pathing('./example_files/staphylococcus/test/test-base-error-2percent.mutations.vcf')

samtools_variant_calls_be2percent = pathing('./example_files/staphylococcus/test-base-error-2percent.bcftools.vcf')
samtools_variant_calls_be10percent = pathing('./example_files/staphylococcus/test-base-error-10percent.bcftools.vcf')<jupyter_output><empty_output><jupyter_text>Train Model VCF file needed to train and assumed to be high-confidence or "truth"<jupyter_code>%%time
deepvcf.train(
    reference_file=ref_train_be2percent,
    alignment_file=bam_train_be2percent,
    vcf_file=vcf_train_be2percent,
    mimimum_alignment_coverage=.75,
    minimum_coverage=30,
    heterozygous_threshold=.25,
    minimum_variant_radius=15,
    # save_pileup_to_destination='./example_files/staphylococcus/train/train.pileup.npy',
    # use_saved_pileup='./example_files/staphylococcus/train/train.pileup.npy',  # Ignore the error, this will fail the first time if it doesn't exist
)<jupyter_output>=== Building Pileup From Scratch ===
=== Pilup Complete ===
=== Tensors Complete ===
Epoch 1/30
360/360 [==============================] - 3s 8ms/step - loss: 0.6732 - base_loss: 0.4467 - genotype_loss: 0.4529 - base_accuracy: 0.5299 - genotype_accuracy: 0.8521 - val_loss: 0.3158 - val_base_loss: 0.2760 - val_genotype_loss: 0.0796 - val_base_accuracy: 0.7369 - val_genotype_accuracy: 0.9817
Epoch 2/30
360/360 [==============================] - 3s 8ms/step - loss: 0.3369 - base_loss: 0.2880 - genotype_loss: 0.0978 - base_accuracy: 0.6474 - genotype_accuracy: 0.9752 - val_loss: 0.2960 - val_base_loss: 0.2638 - val_genotype_loss: 0.0644 - val_base_accuracy: 0.6795 - val_genotype_accuracy: 0.9867
Epoch 3/30
360/360 [==============================] - 3s 7ms/step - loss: 0.3143 - base_loss: 0.2724 - genotype_loss: 0.0837 - base_accuracy: 0.6479 - genotype_accuracy: 0.9805 - val_loss: 0.2955 - val_base_loss: 0.2658 - val_genotype_loss: 0.0593 - val_base_accuracy: 0.6889 - val_genotype_accuracy[...]<jupyter_text>Model Structure<jupyter_code>deepvcf.model.summary()<jupyter_output>Model: "model"
__________________________________________________________________________________________________
Layer (type)                    Output Shape         Param #     Connected to                     
==================================================================================================
input_1 (InputLayer)            [(None, 15, 4, 3)]   0                                            
__________________________________________________________________________________________________
conv2d (Conv2D)                 (None, 15, 4, 16)    208         input_1[0][0]                    
__________________________________________________________________________________________________
max_pooling2d (MaxPooling2D)    (None, 11, 4, 16)    0           conv2d[0][0]                     
__________________________________________________________________________________________________
conv2d_1 (Conv2D)               (None, 11, 4, 32)    4128        max_pooling2d[0][0]          [...]<jupyter_text>Example Pileup & Tensor<jupyter_code>position = 9
print(deepvcf.preprocess.ref_seq[position:position+15])
print('=== Tensor ===')
deepvcf.preprocess.pileup[4:19, :]
position = 9
print(deepvcf.preprocess.ref_seq[position:position+15])
print('=== Tensor ===')
deepvcf.preprocess.get_window(position+8).transpose()<jupyter_output>CATCTTTCTTAAATT
=== Tensor ===<jupyter_text>Model Internal Loss/Validation<jupyter_code>from matplotlib import rcParams, pyplot as plt
import dataframe_image as dfi
import seaborn as sns
# figure size in inches
rcParams['figure.figsize'] = 12,9
rcParams["legend.framealpha"] = 0.0
sns.set_theme()
sns.set(font_scale=1.0) 
snsplot = sns.lineplot(data=pd.DataFrame(deepvcf.model.history.history))
snsplot.set(xlabel='epoch', ylabel='score', title='Model History', xlim=0)
snsplot.figure.savefig('example_files/images/bacterial-model-history.png')<jupyter_output><empty_output><jupyter_text>Use Model to generate a Variant Format File (VCF)<jupyter_code>%%time
predicted_vcf_test = deepvcf.create_vcf(
    reference_file=ref_test_be2percent,
    alignment_file=bam_test_be2percent,
    # save_pileup_to_destination='./example_files/staphylococcus/test/test.pileup.npy',
    # use_saved_pileup='./example_files/staphylococcus/test/test.pileup.npy',  # Ignore the error, this will fail the first time if it doesn't exist
    mimimum_alignment_coverage=.75,
    minimum_coverage=30,
    heterozygous_threshold=.25,
    minimum_variant_radius=15,
    output_folder='./example_files/staphylococcus/test/',  # save as an actual .vcf file
    output_prefix='test' 
)<jupyter_output>=== Building Pileup From Scratch ===
=== Pilup Complete ===
=== Tensors Complete ===
CPU times: user 1min 43s, sys: 698 ms, total: 1min 43s
Wall time: 1min 41s<jupyter_text>Your predicted Variant Calls from the return<jupyter_code>predicted_vcf_test.head()<jupyter_output><empty_output><jupyter_text>Your predicted Variant Calls from the file saved<jupyter_code>predicted_vcf_test = bpd.read_vcf('./example_files/staphylococcus/test/test.deepvcf.vcf')
predicted_vcf_test<jupyter_output><empty_output><jupyter_text>Valitdation Metrics<jupyter_code>bpd.read_vcf(vcf_test_be2percent)
# compare test to predicted
metrics = deepvcf.validation(
    predicted_vcf=predicted_vcf_test,
    real_vcf=vcf_test_be2percent,
)
metrics<jupyter_output><empty_output><jupyter_text>Now we know the complete pipeline lets compare metrics for DeepVCF VS. BCFtools when base errors are 2% and 10%  2% trained - 2% tested<jupyter_code>deepvcf = DeepVCF()
pileup = './example_files/staphylococcus/train/train-base-error-2percent.pileup.npy'
deepvcf.train(
    reference_file=ref_train_be2percent,
    alignment_file=bam_train_be2percent,
    vcf_file=vcf_train_be2percent,
    mimimum_alignment_coverage=.75,
    minimum_coverage=30,
    heterozygous_threshold=.25,
    minimum_variant_radius=15,
    # save_pileup_to_destination=pileup,
    # use_saved_pileup=pileup,  # Ignore the error, this will fail the first time if it doesn't exist
    verbose=0,
)
dvf_be2 = deepvcf.create_vcf(
    reference_file=ref_test_be2percent,
    alignment_file=bam_test_be2percent,
    # save_pileup_to_destination='./example_files/staphylococcus/test/test-base-error-2percent.pileup.npy',
    # use_saved_pileup='./example_files/staphylococcus/test/test-base-error-2percent.pileup.npy',    
    mimimum_alignment_coverage=.75,
    minimum_coverage=30,
    heterozygous_threshold=.25,
    minimum_variant_radius=15,
    output_folder='./example_files',  # save as an actual .vcf file
    output_prefix='demo-be-2-2' 
)
dvf_be2_metrics = deepvcf.validation(real_vcf=vcf_test_be2percent, predicted_vcf=dvf_be2)
dvf_be2_metrics<jupyter_output><empty_output><jupyter_text>2% trained - 10% tested<jupyter_code>deepvcf = DeepVCF()
deepvcf.train(
    reference_file=ref_train_be2percent,
    alignment_file=bam_train_be2percent,
    vcf_file=vcf_train_be2percent,
    mimimum_alignment_coverage=.75,
    minimum_coverage=30,
    heterozygous_threshold=.25,
    minimum_variant_radius=15,
    # save_pileup_to_destination='./example_files/staphylococcus/train/train-base-error-2percent.pileup.npy',
    # use_saved_pileup='./example_files/staphylococcus/train/train-base-error-2percent.pileup.npy',  # Ignore the error, this will fail the first time if it doesn't exist
    verbose=0,
)
dvf_be_2_10 = deepvcf.create_vcf(
    reference_file=ref_test_be10percent,
    alignment_file=bam_test_be10percent,
    # save_pileup_to_destination='./example_files/staphylococcus/test/test-base-error-10percent.pileup.npy',
    # use_saved_pileup='./example_files/staphylococcus/test/test-base-error-10percent.pileup.npy',  # Ignore the error, this will fail the first time if it doesn't exist
    mimimum_alignment_coverage=.75,
    minimum_coverage=30,
    heterozygous_threshold=.25,
    minimum_variant_radius=15,
    output_folder='./example_files/',  # save as an actual .vcf file
    output_prefix='demo-be-2-10' 
)
dvf_be_2_10_metrics = deepvcf.validation(real_vcf=vcf_test_be10percent, predicted_vcf=dvf_be_2_10)
dvf_be_2_10_metrics<jupyter_output><empty_output><jupyter_text>10% trained - 10% tested<jupyter_code>deepvcf = DeepVCF()
deepvcf.train(
    reference_file=ref_train_be10percent,
    alignment_file=bam_train_be10percent,
    vcf_file=vcf_train_be10percent,
    mimimum_alignment_coverage=.75,
    minimum_coverage=30,
    heterozygous_threshold=.25,
    minimum_variant_radius=15,
    # save_pileup_to_destination='./example_files/staphylococcus/train/train-base-error-10percent.pileup.npy',
    # use_saved_pileup='./example_files/staphylococcus/train/train-base-error-10percent.pileup.npy',  # Ignore the error, this will fail the first time if it doesn't exist
    verbose=0,
)
dvf_be10 = deepvcf.create_vcf(
    reference_file=ref_test_be10percent,
    alignment_file=bam_test_be10percent,
    # save_pileup_to_destination='./example_files/staphylococcus/test/test-base-error-10percent.pileup.npy',
    # use_saved_pileup='./example_files/staphylococcus/test/test-base-error-10percent.pileup.npy',  # Ignore the error, this will fail the first time if it doesn't exist
    mimimum_alignment_coverage=.75,
    minimum_coverage=30,
    heterozygous_threshold=.25,
    minimum_variant_radius=15,
    output_folder='./example_files',  # save as an actual .vcf file
    output_prefix='demo-be-10-10' 
)
dvf_be10_metrics = deepvcf.validation(real_vcf=vcf_test_be10percent, predicted_vcf=dvf_be10)
dvf_be10_metrics
deepvcf = DeepVCF()
deepvcf.train(
    reference_file=ref_train_be10percent,
    alignment_file=bam_train_be10percent,
    vcf_file=vcf_train_be10percent,
    mimimum_alignment_coverage=.75,
    minimum_coverage=30,
    heterozygous_threshold=.25,
    minimum_variant_radius=15,
    # save_pileup_to_destination='./example_files/staphylococcus/train/train-base-error-10percent.pileup.npy',
    # use_saved_pileup='./example_files/staphylococcus/train/train-base-error-10percent.pileup.npy',  # Ignore the error, this will fail the first time if it doesn't exist
    verbose=0,
)
dvf_be_10_2 = deepvcf.create_vcf(
    reference_file=ref_test_be2percent,
    alignment_file=bam_test_be2percent,
    # save_pileup_to_destination='./example_files/staphylococcus/test/test-base-error-10percent.pileup.npy',
    # use_saved_pileup='./example_files/staphylococcus/test/test-base-error-2percent.pileup.npy',  # Ignore the error, this will fail the first time if it doesn't exist
    mimimum_alignment_coverage=.75,
    minimum_coverage=30,
    heterozygous_threshold=.25,
    minimum_variant_radius=15,
    output_folder='./example_files',  # save as an actual .vcf file
    output_prefix='demo-be-10-2' 
)
dvf_be_10_2_metrics = deepvcf.validation(real_vcf=vcf_test_be2percent, predicted_vcf=dvf_be_10_2)
dvf_be_10_2_metrics<jupyter_output><empty_output><jupyter_text>BCFtools Variant calling   2% base error<jupyter_code>bcf_be2_metrics = deepvcf.validation(
    predicted_vcf=samtools_variant_calls_be2percent,
    real_vcf=vcf_test_be2percent,
)
pprint(bcf_be2_metrics)<jupyter_output>{'het': {'Accuracy': 0.9932351172579675,
         'F1': 0.9946056101654279,
         'PPV': 0.9998794890335021,
         'Sensitivity': 0.9893870736942523},
 'hom_alt': {'Accuracy': 0.9977351653329307,
             'F1': 0.9964522232734153,
             'PPV': 0.9997626957759849,
             'Sensitivity': 0.9931636020744932}}<jupyter_text>10% base error<jupyter_code>bcf_be10_metrics = deepvcf.validation(
    predicted_vcf=samtools_variant_calls_be10percent,
    real_vcf=vcf_test_be10percent,
)
pprint(bcf_be10_metrics)<jupyter_output>{'het': {'Accuracy': 0.29926002409223884,
         'F1': 0.07222601959444065,
         'PPV': 0.9491017964071856,
         'Sensitivity': 0.03754144954997631},
 'hom_alt': {'Accuracy': 0.6540052651372696,
             'F1': 0.7745650575839255,
             'PPV': 0.8231770833333333,
             'Sensitivity': 0.7313743637204998}}<jupyter_text>Comparting DeepVCF with BCFtools when using noisy data<jupyter_code>metrics = [
    dvf_be2_metrics,
    dvf_be10_metrics,
    dvf_be_2_10_metrics,
    dvf_be_10_2_metrics,
    bcf_be2_metrics,
    bcf_be10_metrics,
]
names = [
    'DVCF-2',
    'DVCF-10',
    'DVCF-2-10',
    'DVCF-10-2',
    'BCF-2',
    'BCF-10'
]
hom_alts = [{**{'tool':name}, **m['hom_alt']} for name, m in zip(names, metrics)]
hets = [{**{'tool': name}, **m['het']} for name, m in zip(names, metrics)]
hom_alt = pd.DataFrame(hom_alts)
het = pd.DataFrame(hets)
hom_alt.style.background_gradient()
het.style.background_gradient()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11840856716822556, 'fraction_numerical': 0.05815372102480683, 'mean_word_length': 5.265392781316348, 'pattern_counts': {'":': 0, '<': 58, '<?xml version=': 0, '>': 58, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 1, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '16', 'dir': 'jupyter-structured-clean-dedup', 'id': '491875', 'license': 'MIT', 'n_tokens_mistral': 5892, 'n_tokens_neox': 5136, 'n_words': 801, 'path': 'jupyter_nb/tutorial.ipynb', 'repo_name': 'tmsincomb/DeepVCF', 'size': '14754'}