<jupyter_start><jupyter_text>Fix reads<jupyter_code>import os

# Fix warning about locale unset
os.environ['LANG'] = 'en_US.UTF-8'

!pushd output/reads; prename 's/initial_//' *.fq.gz; popd<jupyter_output>~/workspace/thesis-data-simulation/jackalope/salmonella/output/reads ~/workspace/thesis-data-simulation/jackalope/salmonella
~/workspace/thesis-data-simulation/jackalope/salmonella<jupyter_text>Jackalope produces reads with non-standard identifiers where pairs of reads don't have matching identifiers. For example:* Pair 1: `@SH08-001-NC_011083-3048632-R/1`* Pair 2: `@SH08-001-NC_011083-3048396-F/2`In order to run snippy, these paired identifiers need to match (except for the `/1` and `/2` suffix).So, I have to replace them all with something unique, but which matches in each pair of files. I do this by replacing the position (I think) with the read number (as it appears in the file). So the above identifiers become:* Pair 1: `@SH08-001-NC_011083-1/1`* Pair 2: `@SH08-001-NC_011083-1/2`<jupyter_code>import glob
import os

files = [os.path.basename(f) for f in glob.glob('output/reads/*.fq.gz')]
!parallel -j 24 -I% 'gzip -d --stdout output/reads/% | perl ../scripts/replace-fastq-header.pl | gzip > output/%' \
    ::: {' '.join(files)}<jupyter_output><empty_output><jupyter_text>Create input file for snippy<jupyter_code>import os
import glob

reference_file = 'input/S_HeidelbergSL476.fasta.gz'

# snippy only runs with uncompressed reference
!gunzip -f -k {reference_file}

reference_file_abs = os.path.abspath('input/S_HeidelbergSL476.fasta')

snippy_out = os.path.abspath('phylogeny')

if not os.path.exists(snippy_out):
    os.mkdir(snippy_out)

with open(f'{snippy_out}/snippy.fofn', 'w') as snippy_fofn:
    directory = 'output'
    for file in glob.glob(f'{directory}/*_R1.fq.gz'):
        sample = os.path.basename(file).rsplit('_R1.fq.gz')[0]
        
        files = [f'{directory}/{sample}_R1.fq.gz', f'{directory}/{sample}_R2.fq.gz']
        files = [os.path.abspath(f) for f in files]
        values = [sample]
        values.extend(files)
        snippy_fofn.write('\t'.join(values)+'\n')
!head -n 1 {snippy_out}/snippy.fofn<jupyter_output>SH14-013	/home/CSCScience.ca/apetkau/workspace/thesis-data-simulation/jackalope/salmonella/output/SH14-013_R1.fq.gz	/home/CSCScience.ca/apetkau/workspace/thesis-data-simulation/jackalope/salmonella/output/SH14-013_R2.fq.gz<jupyter_text>Run snippy<jupyter_code>!conda run --name snippy snippy-multi {snippy_out}/snippy.fofn \
    --reference {reference_file_abs} --cpus 6 > {snippy_out}/snippy-commands-all.sh
!head -n-2 {snippy_out}/snippy-commands-all.sh > {snippy_out}/snippy-commands-variant.sh
!tail -n 2 {snippy_out}/snippy-commands-all.sh > {snippy_out}/snippy-commands-core.sh
!tail -n 2 {snippy_out}/snippy-commands-variant.sh
!echo '****'
!tail {snippy_out}/snippy-commands-core.sh
# Run variant calling in parallel
!(pushd {snippy_out} && conda run --name snippy \
  parallel -j 12 -a {snippy_out}/snippy-commands-variant.sh && popd) > {snippy_out}/snippy-variant.log 2>&1
# Run core in serial
!(pushd {snippy_out} && conda run --name snippy \
  bash {snippy_out}/snippy-commands-core.sh && popd) > {snippy_out}/snippy-core.log 2>&1
!column -s$'\t' -t phylogeny/core.txt<jupyter_output>ID         LENGTH   ALIGNED  UNALIGNED  VARIANT  HET   MASKED  LOWCOV
SH08-001   4888768  4829501  48095      18425    2385  0       8787
SH09-29    4888768  4829610  48255      15585    2378  0       8525
SH10-001   4888768  4829370  47886      16770    2459  0       9053
SH10-002   4888768  4829251  48106      17282    2425  0       8986
SH10-014   4888768  4829822  48272      20812    2467  0       8207
SH10-015   4888768  4829414  48035      20341    2486  0       8833
SH10-30    4888768  4829590  48155      16013    2495  0       8528
SH11-001   4888768  4829479  47904      16234    2446  0       8939
SH11-002   4888768  4829426  48204      20843    2539  0       8599
SH12-001   4888768  4828792  48042      17178    2429  0       9505
SH12-002   4888768  4829176  48203      17406    2448  0       8941
SH12-003   4888768  4829802  48159      17171    2397  0       8410
SH12-004   4888768  4829943  48005      17169    2493  0       8327
SH12-005   4888768  4829544  48008      17165 [...]<jupyter_text>Tree<jupyter_code>!iqtree -redo -s phylogeny/core.aln -T 24 | tail -n 30
!sed -i.bak -e 's/Reference/reference/' phylogeny/core.aln.treefile<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12028725314183124, 'fraction_numerical': 0.14295332136445243, 'mean_word_length': 3.335603112840467, 'pattern_counts': {'":': 0, '<': 18, '<?xml version=': 0, '>': 26, 'https://': 0, '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': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '578085', 'license': 'Apache-2.0', 'n_tokens_mistral': 2156, 'n_tokens_neox': 1726, 'n_words': 379, 'path': 'jackalope/salmonella/run-snippy.ipynb', 'repo_name': 'apetkau/thesis-data-simulation', 'size': '4456'}