<jupyter_start><jupyter_text>Prerequisites- dwgsim https://github.com/nh13/DWGSIM- samtools https://github.com/samtools/samtools You can install easily via anaconda ```bashconda install -y dwgsim samtools bcftools bwa```  This notebook shows how to create In Silico (computational) mutations to give training and testing examples for DeepVCF. We keep track of the mutations introduced via Variant Format Files (VCF). It's not important if the mutatations are found in the wild, but rather to demonstrate how a deep learning model will train on these mutations and determine variant calls. If you would like to see DeepVCF use a high-confidence provided from Genome in a Bottle [GIAB]() then please see the notebook [giab notebook](./giab)<jupyter_code>from DeepVCF import pathing
from DeepVCF import pandas as pd
from DeepVCF import alignment
from DeepVCF import dwgsim<jupyter_output>=== No GPU Detected ===<jupyter_text>Create Reads With In Silco Mutations<jupyter_code>dwgsim(
    ref_file='./example_files/staphylococcus/train/MRSA107.fasta',
    output_folder='./example_files/staphylococcus/train/',
    output_prefix='train-base-error-10percent',
    **{
        '-r': .005,  # mutation rate :: need to get about 20K mutations for usable training,
        '-e': '.1-.1',  # 10-10% error rate forward
        '-E': '.1-.1',  # 10-10% error rare backward
        '-C': 50,  # mean coverage :: to reduce runtime since 50 is plenty for the tensor
        '-1': 150, # read length
        '-2': 150, # read length
    }
)<jupyter_output>dwgsim -r 0.005 -e .1-.1 -E .1-.1 -C 50 -1 150 -2 150 /home/tmsincomb/Dropbox/git/DeepVCF/jupyter_nb/example_files/staphylococcus/train/MRSA107.fasta train-base-error-10percent<jupyter_text>Variant Format File (VCF) With In Silco Mutations<jupyter_code>vcf = pd.read_vcf('./example_files/staphylococcus/train/train-base-error-10percent.mutations.vcf')
print('Mutations Introduced:', vcf.shape[0])
vcf.head()<jupyter_output><empty_output><jupyter_text>Alignment Back to the Genome with the introduced mutations<jupyter_code>pysam_obj = alignment.bwa_mem(
    reference='./example_files/staphylococcus/train/MRSA107.fasta',
    reads=[
        './example_files/staphylococcus/train/train-base-error-10percent.bwa.read1.fastq', 
        './example_files/staphylococcus/train/train-base-error-10percent.bwa.read2.fastq'
    ],
    output_folder='./example_files/staphylococcus/train/',
    output_prefix='train-base-error-10percent',
)<jupyter_output>bwa index /home/tmsincomb/Dropbox/git/DeepVCF/jupyter_nb/example_files/staphylococcus/train/MRSA107.fasta
bwa mem -t 6 /home/tmsincomb/Dropbox/git/DeepVCF/jupyter_nb/example_files/staphylococcus/train/MRSA107.fasta /home/tmsincomb/Dropbox/git/DeepVCF/jupyter_nb/example_files/staphylococcus/train/train-base-error-10percent.bwa.read1.fastq /home/tmsincomb/Dropbox/git/DeepVCF/jupyter_nb/example_files/staphylococcus/train/train-base-error-10percent.bwa.read2.fastq > /home/tmsincomb/Dropbox/git/DeepVCF/jupyter_nb/example_files/staphylococcus/train/train-base-error-10percent.sam
samtools view -Shb /home/tmsincomb/Dropbox/git/DeepVCF/jupyter_nb/example_files/staphylococcus/train/train-base-error-10percent.sam > /home/tmsincomb/Dropbox/git/DeepVCF/jupyter_nb/example_files/staphylococcus/train/train-base-error-10percent.unsorted.bam
samtools sort -n /home/tmsincomb/Dropbox/git/DeepVCF/jupyter_nb/example_files/staphylococcus/train/train-base-error-10percent.unsorted.bam -o /home/tmsincomb/Dropbo[...]<jupyter_text>Sequence Alignment Mapping (SAM/BAM/CRAM) to be used for DeepVCFs tensor input<jupyter_code>bam = pd.read_seq('./example_files/staphylococcus/train/train-base-error-10percent.mapped.bam', format='bam')
bam.head()<jupyter_output><empty_output><jupyter_text>Create Testing Data :: Verbose = False to show silent option<jupyter_code># Quietly create testing
prefix='test-base-error-10percent'
folder_path=pathing('./example_files/staphylococcus/test')

dwgsim(
    ref_file=folder_path / 'GCF_000013425.1_ASM1342v1_genomic_reference.fasta',
    output_folder=folder_path,
    output_prefix=prefix,
    verbose=False,
    **{
        '-r': .005,  # mutation rate :: need to get about 20K mutations for usable training,
        '-e': '.1-.1',  # 10-10% error rate forward :: for attempt of PacBio
        '-E': '.1-.1',  # 10-10% error rare backward
        '-C': 50,  # mean coverage
        '-1': 150, # read length
        '-2': 150, # read length
    }
)
_ = alignment.bwa_mem(
    reference=folder_path / 'GCF_000013425.1_ASM1342v1_genomic_reference.fasta',
    reads=[
        folder_path / f'{prefix}.bwa.read1.fastq', 
        folder_path / f'{prefix}.bwa.read2.fastq'
    ],
    output_folder=folder_path,
    output_prefix=prefix,
    verbose=False,
)<jupyter_output><empty_output><jupyter_text>Create a comparison test with default 2% platform base error rate<jupyter_code>### DEFAULT TRAIN ###
prefix='train-base-error-2percent'
folder_path=pathing('./example_files/staphylococcus/train')

dwgsim(
    ref_file=folder_path / 'MRSA107.fasta',
    output_folder=folder_path,
    output_prefix=prefix,
    verbose=False,
    **{
        '-r': .005,  # mutation rate :: need to get about 20K mutations for usable training,
        '-C': 50,  # mean coverage
        '-1': 150, # read length
        '-2': 150, # read length
    }
)
_ = alignment.bwa_mem(
    reference=folder_path / 'MRSA107.fasta',
    reads=[
        folder_path / f'{prefix}.bwa.read1.fastq', 
        folder_path / f'{prefix}.bwa.read2.fastq'
    ],
    output_folder=folder_path,
    output_prefix=prefix,
    verbose=False,
)

### DEFAULT TEST ###
prefix='test-base-error-2percent'
folder_path=pathing('./example_files/staphylococcus/test')

dwgsim(
    ref_file=folder_path / 'GCF_000013425.1_ASM1342v1_genomic_reference.fasta',
    output_folder=folder_path,
    output_prefix=prefix,
    verbose=False,
    **{
        '-r': .005,  # mutation rate :: need to get about 20K mutations for usable training,
        '-e': '.1-.1',  # 10-10% error rate forward :: for attempt of PacBio
        '-E': '.1-.1',  # 10-10% error rare backward
        '-C': 50,  # mean coverage
        '-1': 150, # read length
        '-2': 150, # read length
    }
)
_ = alignment.bwa_mem(
    reference=folder_path / 'GCF_000013425.1_ASM1342v1_genomic_reference.fasta',
    reads=[
        folder_path / f'{prefix}.bwa.read1.fastq', 
        folder_path / f'{prefix}.bwa.read2.fastq'
    ],
    output_folder=folder_path,
    output_prefix=prefix,
    verbose=False,
)<jupyter_output><empty_output><jupyter_text>We Now want to compare out variant calls to the Samtools/BCFTools variant caller.<jupyter_code>from DeepVCF import bcftools_vcf
samtools_variant_calls_be2percent = bcftools_vcf(
    ref_file='./example_files/staphylococcus/test/GCF_000013425.1_ASM1342v1_genomic_reference.fasta', 
    bam_file='./example_files/staphylococcus/test/test-base-error-2percent.mapped.bam',
    output_folder='./example_files/staphylococcus/',
    output_prefix='GCF_000013425-base-error-2percent'
)
samtools_variant_calls_be10percent = bcftools_vcf(
    ref_file='./example_files/staphylococcus/test/GCF_000013425.1_ASM1342v1_genomic_reference.fasta', 
    bam_file='./example_files/staphylococcus/test/test-base-error-10percent.mapped.bam',
    output_folder='./example_files/staphylococcus/',
    output_prefix='GCF_000013425-base-error-10percent'
)<jupyter_output>=== No GPU Detected ===<jupyter_text>Our Data is ready!<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/GCF_000013425-base-error-2percent.bcftools.vcf')
samtools_variant_calls_be10percent = pathing('./example_files/staphylococcus/GCF_000013425-base-error-10percent.bcftools.vcf')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1247929320817228, 'fraction_numerical': 0.04295969077857537, 'mean_word_length': 4.603960396039604, 'pattern_counts': {'":': 0, '<': 33, '<?xml version=': 0, '>': 35, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '9', 'dir': 'jupyter-structured-clean-dedup', 'id': '333750', 'license': 'MIT', 'n_tokens_mistral': 3557, 'n_tokens_neox': 3127, 'n_words': 507, 'path': 'jupyter_nb/creating-data-for-usage-demo.ipynb', 'repo_name': 'tmsincomb/DeepVCF', 'size': '9055'}