<jupyter_start><jupyter_text>Execute profile pipeline on BGI250 data The aim of this notebook is to run the metagenome profile pipeline using the defaults and the custom databases on real samples from the BGI250 dataset  Init<jupyter_code>import os
import pandas as pd<jupyter_output>/ebio/abt3_projects/software/miniconda3_gt4.4/envs/py3/lib/python3.6/importlib/_bootstrap.py:219: RuntimeWarning: numpy.dtype size changed, may indicate binary incompatibility. Expected 96, got 88
  return f(*args, **kwds)<jupyter_text>Var<jupyter_code># Dirs
work_dir = "/ebio/abt3_projects/small_projects/jdelacuesta/DBs_benchmark"
pipeline_folder = os.path.join(work_dir, "bin/llmgp")<jupyter_output><empty_output><jupyter_text>Create folders<jupyter_code># SGE out folders
SGE_GTDB_dir = os.path.join(work_dir, "tmp/SGE_out/llmgp/GTDB")
SGE_proGenomes_dir = os.path.join(work_dir, "tmp/SGE_out/llmgp/progenomes")
SGE_Defaults_dir = os.path.join(work_dir, "tmp/SGE_out/llmgp/defaults")

if not os.path.exists(os.path.join(work_dir, "tmp/SGE_out/llmgp")):
    os.makedirs(SGE_GTDB_dir)
    os.makedirs(SGE_proGenomes_dir)
    os.makedirs(SGE_Defaults_dir)
# Profile out folders
out_GTDB_dir = os.path.join(work_dir, "data/profiles_BGI250/GTDB")
out_proGenomes_dir = os.path.join(work_dir, "data/profiles_BGI250/progenomes")
out_Defaults_dir = os.path.join(work_dir, "data/profiles_BGI250/defaults")

if not os.path.exists(os.path.join(work_dir, "data/profiles_BGI250")):
    os.makedirs(out_GTDB_dir)
    os.makedirs(out_proGenomes_dir)
    os.makedirs(out_Defaults_dir)<jupyter_output><empty_output><jupyter_text>Prepare config files<jupyter_code>samples_file = "/ebio/abt3_projects/databases_no-backup/TUK/metagenome/190100_TUK-all/samples_filt_BGI250.txt"
# Note that I will use the cutsom HUMANn2 database, built using progenomes
config_file = os.path.join(pipeline_folder, 'config_custom-db.yaml')
#!cat $config_file
config_file = os.path.join(pipeline_folder, 'config.yaml')
!cat $config_file
config_default="""#-- I/O --#
# table with sample --> read_file information
samples_file: {samples_file}

# output location
output_dir: {output_dir}

# read file path
## use "None" if full file path is included in the samples_file
read_file_path: None

#-- DB --#
# NOTE: see the config_custom-db.yaml for using the progenomes db
# humann2-associated databases
genefamily_annotation_db: {genefamily_annotation_db}
humann2_nuc_db: {humann2_nuc_db}
humann2_prot_db: {humann2_prot_db}
metaphlan2_pkl_db: /ebio/abt3_projects/databases_no-backup/metaphlan2/mpa_v20_m200/mpa_v20_m200.pkl
metaphlan2_bt2_db: /ebio/abt3_projects/databases_no-backup/metaphlan2/mpa_v20_m200/mpa_v20_m200     # use the prefix for the *.bt2 files
utility_mapping_db: /ebio/abt3_projects/databases_no-backup/humann2/utility_mapping
## kraken/bracken (db selected automatically based on read length)
kraken_dbs:
  150bp: {kraken_dbs_150}
  100bp: {kraken_dbs_100}
# taxonomy
tax_dump: /ebio/abt3_projects/databases/Kraken/taxonkit/names.dmp

#-- subsample --#
# subsampling input reads 
## "Skip" skips subsampling; otherwise set the number of reads to subsample
subsample_depth: 1000000
subsample_seed: 18938

#-- include read2 (if paired-end) --#
# combine R1 & R2 or just use R1?
include_read2: True

#-- humann2 temporary files --#
# remove the large temporary files generated by humann2?
rm_humann2_tmp_files: True      

#-- humann2 groupings --#
# always have at least the "*_default" grouping
humann2_regroup:
  - uniref50_default
  - uniref50_go
  - uniref50_ko
  - uniref50_eggnog
  - uniref50_pfam
  - uniref50_level4ec
  - uniref50_infogo1000
  - uniref50_rxn

#-- software parameters --#
# Use "Skip" to skip steps.
# By skipping, you can run just humann2, kraken/bracken, and/or simka
params:
  # humann2
  metaphlan2: -t rel_ab
  humann2: --gap-fill on --diamond-2pass --search-mode uniref50 # --bypass-nucleotide-index 
  humann2_db_in_memory: Skip #True        # copy databases to memory; less I/O, more memory
  humann2_diamond: --sensitive --max-target-seqs 20 --block-size 3 --index-chunks 2
  humann2_diamond_evalue: 1
  reduce_taxonomic_profile: --function max --sort-by level
  humann2_renorm_table: --units relab
  # kraken/bracken (NOTE: dependent on read length)
  kraken: ""        
  bracken:  -t 100 -l S        # species level (S); `-r` parameter set automatically
  # simka
  simka: Skip #-kmer-size 31 -abundance-min 2 -simple-dist -max-reads 1000000
  simka_vis: Skip #-width 8 -height 8 -pca -heatmap
  # hulk
  hulk_histosketch: Skip #-k 21 -m 2
  hulk_distance: 
    - jaccard
    - braycurtis

#-- snakemake pipeline --#
## To use /tmp/global2/, see http://ilm.eb.local/user-guide/#Scratch-space-on-_002ftmp_002fglobal2
pipeline:
  snakemake_folder: ./
  script_folder: ./bin/scripts/
  temp_folder: /tmp/global2/      # your username will be added automatically to this path
"""
config_custom="""# DESCRIPTION:
## This is an example of running the pipeline with a custom humann2 db.
## This config is set up to just use the custom nucleotide db, but the protein db could be used also (or instead)

#-- I/O --#
# table with sample --> read_file information
samples_file: {samples_file}

# output location
output_dir: {output_dir}

# read file path
# use "None" if full file path is included in the samples_file
read_file_path: None

#-- DB --#
## humann2 
### custom humann2 databases
humann2_nuc_db: {humann2_nuc_db}
humann2_prot_db: {humann2_prot_db}
### required humann2 database files (no need to change this)
### To use UniRef90, change 
### genefamily_annotation_db: /ebio/abt3_projects/databases_no-backup/humann2/utility_mapping/map_uniref90_name.txt.bz2
genefamily_annotation_db: /ebio/abt3_projects/databases_no-backup/humann2/utility_mapping/map_uniref50_name.txt.bz2
metaphlan2_pkl_db: /ebio/abt3_projects/databases_no-backup/metaphlan2/mpa_v20_m200/mpa_v20_m200.pkl
metaphlan2_bt2_db: /ebio/abt3_projects/databases_no-backup/metaphlan2/mpa_v20_m200/mpa_v20_m200     
utility_mapping_db: /ebio/abt3_projects/databases_no-backup/humann2/utility_mapping
## kraken/bracken (db selected automatically based on read length)
kraken_dbs:
  150bp: {kraken_dbs_150}
  100bp: {kraken_dbs_100}
### NCBI taxonomy
tax_dump: /ebio/abt3_projects/databases_no-backup/GTDB/release86/LLMGP-DB/kraken2/taxonomy/names.dmp

#-- subsample --#
# subsampling input reads 
## "Skip" skips subsampling; otherwise set the number of reads to subsample
subsample_depth: 1000000
subsample_seed: 18938

#-- include read2 (if paired-end) --#
# combine R1 & R2 or just use R1?
include_read2: True

#-- humann2 temporary files --#
# remove the large temporary files generated by humann2?
rm_humann2_tmp_files: True      

#-- humann2 groupings --#
# always have at least the "*_default" grouping
humann2_regroup:
  - uniref50_default
  - uniref50_go
  - uniref50_ko
  - uniref50_eggnog
  - uniref50_pfam
  - uniref50_level4ec
  - uniref50_infogo1000
  - uniref50_rxn

#-- software parameters --#
# Use "Skip" to skip steps.
# By skipping, you can run just humann2, kraken/bracken, or simka
params:
  # humann2
  metaphlan2: -t rel_ab  
  humann2: --gap-fill on --bypass-nucleotide-index --diamond-2pass --search-mode uniref50
  humann2_db_in_memory: Skip #True        # copy databases to memory; less I/O, more memory
  humann2_diamond: --sensitive --max-target-seqs 20 --block-size 3 --index-chunks 2
  humann2_diamond_evalue: 1
  reduce_taxonomic_profile: --function max --sort-by level
  humann2_renorm_table: --units relab
  # kraken/bracken (NOTE: dependent on read length)
  kraken: ""
  bracken: -t 100 -l S        # species level (S); `-r` parameter set automatically
  # simka 
  simka: Skip # -kmer-size 31 -abundance-min 2 -simple-dist -max-reads 1000000
  simka_vis: Skip # -width 8 -height 8 -pca -heatmap
  # hulk
  hulk_histosketch: Skip # -k 21 -m 2
  hulk_distance: 
    - jaccard
    - braycurtis


#-- snakemake pipeline --#
pipeline:
  snakemake_folder: ./
  script_folder: ./bin/scripts/
  temp_folder: /tmp/global2/        # your username will be added automatically to this path
"""
# Config with GTDB databases
config_GTDB = config_custom.format(samples_file = samples_file, 
                                   output_dir = out_GTDB_dir, 
                                   humann2_nuc_db = "/ebio/abt3_projects/databases_no-backup/GTDB/release86/LLMGP-DB/humann2/all_genes_annot.fna.gz", 
                                   humann2_prot_db = "/ebio/abt3_projects/databases_no-backup/GTDB/release86/LLMGP-DB/humann2/all_genes.dmnd", 
                                   kraken_dbs_150 = "/ebio/abt3_projects/databases_no-backup/GTDB/release86/LLMGP-DB/kraken2/database150mers.kraken", 
                                   kraken_dbs_100 = "/ebio/abt3_projects/databases_no-backup/GTDB/release86/LLMGP-DB/kraken2/database100mers.kraken") 

# Write config file
config_GTDB_file = os.path.join(pipeline_folder, 'BGI250_GTDB.yaml')
with open(config_GTDB_file, 'w') as outF:
    outF.write(config_GTDB)
#!cat $config_GTDB_file
# Config with proGenomes databases
config_proGenomes = config_custom.format(samples_file = samples_file, 
                                   output_dir = out_proGenomes_dir, 
                                   humann2_nuc_db = "/ebio/abt3_projects/databases/humann2_progenomes/progenomes_HUMANn2_UniRef.fna", 
                                   humann2_prot_db = "/ebio/abt3_projects/databases/humann2_progenomes/progenomes_HUMANn2.dmnd", 
                                   kraken_dbs_150 = "/ebio/abt3_projects/databases/Kraken/K2_Progenomes/Kraken/150mers/database150mers.kraken", 
                                   kraken_dbs_100 = "/ebio/abt3_projects/databases/Kraken/K2_Progenomes/Kraken/100mers/database100mers.kraken")

# Write config file
config_proGenomes_file = os.path.join(pipeline_folder, 'BGI250_proGenomes.yaml')
with open(config_proGenomes_file, 'w') as outF:
    outF.write(config_proGenomes)
#!cat $config_proGenomes_file
# Config with default databases
# "/ebio/abt3_projects/databases_no-backup/humann2/uniref90_ec_filt/uniref90.ec_filtered.1.1.dmnd"
config_Defaults = config_default.format(samples_file = samples_file, 
                                        output_dir = out_Defaults_dir, 
                                        humann2_nuc_db = "/ebio/abt3_projects/databases_no-backup/humann2/chocophlan/README.md", 
                                        humann2_prot_db = "/ebio/abt3_projects/databases_no-backup/humann2/uniref50/uniref50_annotated.1.1.dmnd",
                                        genefamily_annotation_db = "/ebio/abt3_projects/databases_no-backup/humann2/utility_mapping/map_uniref50_name.txt.bz2",
                                        kraken_dbs_150 = "/ebio/abt3_projects/databases/Kraken/K2_Standard/150mers/database150mers.kraken",
                                        kraken_dbs_100 = "/ebio/abt3_projects/databases/Kraken/K2_Standard/100mers/database100mers.kraken")


# Write config file
config_Defaults_file = os.path.join(pipeline_folder, 'BGI250_Defaults.yaml')
with open(config_Defaults_file, 'w') as outF:
    outF.write(config_Defaults)
#!cat $config_Defaults_file<jupyter_output><empty_output><jupyter_text>Prepare snakemake command<jupyter_code># GTDB
#conda_env = 'conda activate snakemake'
conda_env = 'conda activate snakemake_dev'
P_cmd = "cd {llmgqc} && {conda_env} && screen -L -S BGI_llmgp {exe} {config_file} \
    cluster.json {SGE_out} {jobs} \
    --keep-going --rerun-incomplete --dryrun"
GTDB_cmd = P_cmd.format(conda_env = conda_env, 
                        llmgqc = pipeline_folder,  
                        exe = './snakemake_sge.sh', 
                        config_file = config_GTDB_file, 
                        SGE_out = SGE_GTDB_dir, 
                        jobs = 5)
print(GTDB_cmd)
proGenomes_cmd = P_cmd.format(conda_env = conda_env, 
                        llmgqc = pipeline_folder,  
                        exe = './snakemake_sge.sh', 
                        config_file = config_proGenomes_file, 
                        SGE_out = SGE_proGenomes_dir, 
                        jobs = 10)
print(proGenomes_cmd)
Defaults_cmd = P_cmd.format(conda_env = conda_env, 
                        llmgqc = pipeline_folder,  
                        exe = './snakemake_sge.sh', 
                        config_file = config_Defaults_file, 
                        SGE_out = SGE_Defaults_dir, 
                        jobs = 20)
print(Defaults_cmd)<jupyter_output>cd /ebio/abt3_projects/small_projects/jdelacuesta/DBs_benchmark/bin/llmgp && conda activate snakemake_dev && screen -L -S BGI_llmgp ./snakemake_sge.sh /ebio/abt3_projects/small_projects/jdelacuesta/DBs_benchmark/bin/llmgp/BGI250_Defaults.yaml     cluster.json /ebio/abt3_projects/small_projects/jdelacuesta/DBs_benchmark/tmp/SGE_out/llmgp/defaults 20     --keep-going --rerun-incomplete --dryrun<jupyter_text>Session Info<jupyter_code>sessionInfo = "find {0} -name '*.yaml' | xargs head -n 1000".format(os.path.join(pipeline_folder, 'bin', 'envs'))
!$sessionInfo<jupyter_output>==> /ebio/abt3_projects/small_projects/jdelacuesta/DBs_benchmark/bin/llmgp/bin/envs/bowtie2.yaml <==
channels:
- conda-forge
- bioconda
dependencies:
- pigz
- bioconda::samtools
- bioconda::bedtools
- bioconda::bowtie2

==> /ebio/abt3_projects/small_projects/jdelacuesta/DBs_benchmark/bin/llmgp/bin/envs/compress.yaml <==
channels:
- conda-forge
- bioconda
dependencies:
- pigz
- bioconda::dsrc

==> /ebio/abt3_projects/small_projects/jdelacuesta/DBs_benchmark/bin/llmgp/bin/envs/fastqc.yaml <==
channels:
- conda-forge
- bioconda
dependencies:
- bioconda::fastqc

==> /ebio/abt3_projects/small_projects/jdelacuesta/DBs_benchmark/bin/llmgp/bin/envs/hadley.yaml <==
channels:
- conda-forge
dependencies:
- conda-forge::r-ape
- conda-forge::r-dplyr
- conda-forge::r-tidyr
- conda-forge::r-ggplot2

==> /ebio/abt3_projects/small_projects/jdelacuesta/DBs_benchmark/bin/llmgp/bin/envs/hulk.yaml <==
channels:
- conda-forge
- bioconda
dependencies:
- pigz
- bioconda[...]
{'fraction_non_alphanumeric': 0.10394137734775186, 'fraction_numerical': 0.032939669891861124, 'mean_word_length': 3.4287964713295525, 'pattern_counts': {'":': 0, '<': 27, '<?xml version=': 0, '>': 29, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 16, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '6', 'dir': 'jupyter-structured-clean-dedup', 'id': '330448', 'license': 'MIT', 'n_tokens_mistral': 5571, 'n_tokens_neox': 5187, 'n_words': 918, 'path': 'benchmark/05_BGI_profiles.ipynb', 'repo_name': 'leylabmpi/Struo', 'size': '14056'}