<jupyter_start><jupyter_text>Create standard Kraken/Bracken databases The aim of this notebook is to create standard `Kraken 2` and `Bracken 2` databases of RefSeq *Archaea* and *Bacteria* genomes.   Init<jupyter_code>import os<jupyter_output><empty_output><jupyter_text>Var<jupyter_code># Conda env
Bracken_env = "Bracken2"
# Scripts dir
scripts_dir = "/ebio/abt3_projects/small_projects/jdelacuesta/DBs_benchmark/scripts"
SGE_dir = os.path.join(scripts_dir, "SGE_out")
if not os.path.exists(scripts_dir):
    os.makedirs(scripts_dir)
    os.makedirs(SGE_dir)
# Kraken dir
kraken_dbs = "/tmp/global/jdelacuesta/standard_DBs/kraken"
if not os.path.exists(kraken_dbs):
    os.makedirs(kraken_dbs)<jupyter_output><empty_output><jupyter_text>Build standard databases  Kraken<jupyter_code># Download taxonomy
taxonomy_cmd = "kraken2-build --download-taxonomy --db {k_db}"
taxonomy_job = taxonomy_cmd.format(k_db = kraken_dbs)
taxonomy_job = 'bash -c "source activate {}; {}"'.format(Bracken_env, taxonomy_job)
print(taxonomy_job)
# If there are errors downloading the taxonomy
# Use already dowloaded files
NCBI_taxonomy = "/ebio/abt3_projects/databases/Kraken/K2_Progenomes/Kraken/taxonomy"
rsync_cmd = "rsync -ah {0} {1}".format(NCBI_taxonomy, kraken_dbs)

# Only sync if there is no taxonomy folder
if not os.path.exists(os.path.join(kraken_dbs, "taxonomy")): 
    print(rsync_cmd)
    !$rsync_cmd
# Download Kraken 2 databases
download_cmd = "kraken2-build --download-library {domain} --db {k_db}"

# Download Bacteria DB
download_bacteria = download_cmd.format(domain = "bacteria", k_db = kraken_dbs)
download_bacteria = 'bash -c "source activate {}; {}"'.format(Bracken_env, download_bacteria)
print(download_bacteria)

# Download Archaea DB
download_archaea = download_cmd.format(domain = "archaea", k_db = kraken_dbs)
download_archaea = 'bash -c "source activate {}; {}"'.format(Bracken_env, download_archaea)
print(download_archaea)
# Build the database
kraken_build = """#!/bin/bash
#$ -N {name}
#$ -pe parallel {cpu}
#$ -l h_vmem=10G
#$ -l h_rt=200:0:0
#$ -o {SGE_dir}
#$ -j y
#$ -wd {workdir}
#$ -m ea
#$ -M <EMAIL>

export PATH='/ebio/abt3_projects/software/miniconda3_gt4.4/envs/Bracken2/bin':$PATH
/ebio/abt3_projects/software/miniconda3_gt4.4/envs/Bracken/bin/kmer2read_distr

kraken2-build --build --db {k_db} --threads {cpu}
"""
job_name = "kraken2_build.sh"
kraken_build_file = os.path.join(scripts_dir, job_name)
script_build = kraken_build.format(name = job_name,
                                         workdir = kraken_dbs,
                                         cpu = 30,
                                         k_db = kraken_dbs, 
                                         SGE_dir = SGE_dir)
   
print(script_build)

with open(kraken_build_file, "w") as f:
    f.write(script_build) 
!qsub $kraken_build_file<jupyter_output>Your job 4199461 ("kraken2_build.sh") has been submitted<jupyter_text>Bracken<jupyter_code># Create output directories
dir_100mer = os.path.join(kraken_dbs, "100mers")
dir_150mer = os.path.join(kraken_dbs, "150mers")
if not os.path.exists(dir_100mer):
    os.makedirs(dir_100mer)
    os.makedirs(dir_150mer)
# Create symlinks to kraken files
k2d_files = ["hash.k2d",  "opts.k2d",  "taxo.k2d"]
for k2d in k2d_files:
    file_path = os.path.join("..", k2d)
    link_100mer_path = os.path.join(dir_100mer, k2d)
    link_150mer_path = os.path.join(dir_150mer, k2d)
    symlink_cmd = 'ln -s {file} {link}'
    symlink_100mer_job = symlink_cmd.format(file = file_path, link = link_100mer_path)
    symlink_150mer_job = symlink_cmd.format(file = file_path, link = link_150mer_path)
    print(symlink_100mer_job, symlink_150mer_job)
    !$symlink_100mer_job; $symlink_150mer_job
# Build the database
bracken_build = """#!/bin/bash
#$ -N {name}
#$ -pe parallel {cpu}
#$ -l h_vmem=10G
#$ -l h_rt=200:0:0
#$ -o {SGE_dir}
#$ -j y
#$ -wd {workdir}
#$ -m ea
#$ -M <EMAIL>

export PATH='/ebio/abt3_projects/software/miniconda3_gt4.4/envs/Bracken2/bin':$PATH

# Create database.kraken
kraken2 --db={k_db} --threads={cpu} <( find -L {k_db} \( -name "*.fna" -o -name "*.fa" -o -name "*.fasta" \) -exec cat {{}} + ) > {kraken_db_file}

# database 100mer
# kmer2read
{kmer2read} --seqid2taxid {k_db}/seqid2taxid.map \
    --taxonomy {k_db}/taxonomy \
    --kraken {kraken_db_file} \
    --output {output_100mer_2read} \
    -k 35 \
    -l 100 \
    -t {cpu}
    
# kmer distribution    
generate_kmer_distribution.py -i {output_100mer_2read} -o {output_100mer_distr}

# database 150mer
# kmer2read
{kmer2read} --seqid2taxid {k_db}/seqid2taxid.map \
    --taxonomy {k_db}/taxonomy \
    --kraken {kraken_db_file} \
    --output {output_150mer_2read} \
    -k 35 \
    -l 150 \
    -t {cpu}
    
# kmer distribution    
generate_kmer_distribution.py -i {output_150mer_2read} -o {output_150mer_distr}
"""
# Files and paths
kmer2read_path = "/ebio/abt3_projects/software/miniconda3_gt4.4/envs/Bracken/bin/kmer2read_distr"

kraken_database_file = os.path.join(kraken_dbs, "database.kraken")

# 100mers
k2read_100mer = os.path.join(dir_100mer, "database100mers.kraken")
k2distr_100mer = os.path.join(dir_100mer, "database100mers.kmer_distrib")

# 150mers
k2read_150mer = os.path.join(dir_150mer, "database150mers.kraken")
k2distr_150mer = os.path.join(dir_150mer, "database150mers.kmer_distrib")
job_name = "bracken2_build.sh"
bracken_build_file = os.path.join(scripts_dir, job_name)
script_build = bracken_build.format(name = job_name,
                                    workdir = kraken_dbs,
                                    cpu = 30,
                                    k_db = kraken_dbs, 
                                    SGE_dir = SGE_dir, 
                                    kraken_db_file = kraken_database_file, 
                                    kmer2read = kmer2read_path, 
                                    output_100mer_2read = k2read_100mer, 
                                    output_100mer_distr = k2distr_100mer,
                                    output_150mer_2read = k2read_150mer, 
                                    output_150mer_distr = k2distr_150mer)
print(script_build)
with open(bracken_build_file, "w") as f:
    f.write(script_build) 
!qsub $bracken_build_file<jupyter_output>Your job 4216709 ("bracken2_build.sh") has been submitted<jupyter_text>Move database to databases folder<jupyter_code>databases_project = os.path.join("/ebio/abt3_projects/databases/Kraken/K2_Standard")
if not os.path.exists(databases_project):
    os.makedirs(databases_project)

rsync_kraken = "rsync -ah --bwlimit 300m {0}/ {1}".format(kraken_dbs, databases_project)
print(rsync_kraken)
#!$rsync_kraken<jupyter_output>rsync -ah --bwlimit 300m /tmp/global/jdelacuesta/standard_DBs/kraken/ /ebio/abt3_projects/databases/Kraken/K2_Standard<jupyter_text>Session info<jupyter_code>!conda list -n Bracken2<jupyter_output># packages in environment at /ebio/abt3_projects/software/miniconda3_gt4.4/envs/Bracken2:
#
# Name                    Version                   Build  Channel
appdirs                   1.4.3                      py_1    conda-forge
asn1crypto                0.24.0                     py_1    conda-forge
attrs                     18.1.0                     py_1    conda-forge
automat                   0.7.0                    py36_0    conda-forge
backcall                  0.1.0                      py_0    conda-forge
blas                      1.1                    openblas    conda-forge
blast                     2.7.1                h4422958_6    bioconda
bleach                    2.1.3                      py_0    conda-forge
boost                     1.67.0           py36h3e44d54_0    conda-forge
boost-cpp                 1.67.0               h3a22d5f_0    conda-forge
bracken                   2.2              py36h2d50403_0    bioconda/label/broken
bzip2            [...]
{'fraction_non_alphanumeric': 0.11215666327568667, 'fraction_numerical': 0.04488809766022381, 'mean_word_length': 2.697696285848613, 'pattern_counts': {'":': 0, '<': 24, '<?xml version=': 0, '>': 24, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 40, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '6', 'dir': 'jupyter-structured-clean-dedup', 'id': '49411', 'license': 'MIT', 'n_tokens_mistral': 3175, 'n_tokens_neox': 2968, 'n_words': 508, 'path': 'benchmark/01_Standard-Kraken-Bracken-DBs.ipynb', 'repo_name': 'leylabmpi/Struo', 'size': '7906'}