<jupyter_start><jupyter_text>GoalJacobo de la Cuesta-Zuluaga.The goal of this notebook is to generate phylogenomic trees of *Methanomassiliicoccales* using `phylophlan`. Note that the version used here is still in development, so there are several workarounds to make it run. I will try to document them along the notebook and perform them here, however, there are a couple that should be run directly from the terminal.This assumes that the `phylophlan` repo has been cloned and the conda environment with the requirements is ready.  Download phylophlanI downloaded the development branch of `phylophlan` according to the instructions [here](https://bitbucket.org/nsegata/phylophlan/wiki/phylophlan2markdown-header-conda-environment-medium) as follows:```bashhg clone https://bitbucket.org/nsegata/phylophlancd phylophlanhg up dev```also see [this issue](https://bitbucket.org/nsegata/phylophlan/issues/27/phylophlan2-implementation) in the bitbucket repo  Init<jupyter_code>import os<jupyter_output><empty_output><jupyter_text>Var<jupyter_code># Phylophlan env
phylophlan_env = "py3_phylophlan"

# Work dir
work_dir = "/ebio/abt3_projects/vadinCA11/data/V11"

# phylophlan paths
phylophlan_dir = "/ebio/abt3_projects/vadinCA11/bin/phylophlan"
phylophlan_exe = os.path.join(phylophlan_dir, "phylophlan2.py")

# Download databases
phylophlan_db = os.path.join(phylophlan_dir, "phylophlan_databases")
if not os.path.exists(phylophlan_db):
    os.makedirs(phylophlan_db)
    
wget_cmd = "wget {0} -O {1}"

phylophlan_tar = os.path.join(phylophlan_db, "PhyloPhlAn.tar")
phylophlan_md5 = os.path.join(phylophlan_db, "PhyloPhlAn.md5")

wget_tar = wget_cmd.format("https://bitbucket.org/nsegata/phylophlan/downloads/phylophlan.tar", 
                    phylophlan_tar)
wget_md5 = wget_cmd.format("https://bitbucket.org/nsegata/phylophlan/downloads/phylophlan.md5", 
                    phylophlan_md5)

print(wget_tar)
!$wget_tar
print(wget_md5)
!$wget_md5<jupyter_output>wget https://bitbucket.org/nsegata/phylophlan/downloads/phylophlan.tar -O /ebio/abt3_projects/vadinCA11/bin/phylophlan/phylophlan_databases/PhyloPhlAn.tar
--2019-01-24 15:36:34--  https://bitbucket.org/nsegata/phylophlan/downloads/phylophlan.tar
Resolving bitbucket.org (bitbucket.org)... 192.168.127.12, 192.168.127.12, 192.168.127.12, ...
Connecting to bitbucket.org (bitbucket.org)|192.168.127.12|:443... connected.
HTTP request sent, awaiting response... 302 Found
Location: https://bbuseruploads.s3.amazonaws.com/571bc798-2b2d-49f2-81c9-18a331937270/downloads/e675ba50-0946-4bc6-9daf-5da24e8eca08/phylophlan.tar?Signature=cxpvqKniAScM1Av43t5U%2BxqFbwI%3D&Expires=1548342395&AWSAccessKeyId=AKIAIQWXW6WLXMB5QZAQ&versionId=<KEY>9F&response-content-disposition=attachment%3B%20filename%3D%22phylophlan.tar%22 [following]
--2019-01-24 15:36:35--  https://bbuseruploads.s3.amazonaws.com/571bc798-2b2d-49f2-81c9-18a331937270/downloads/e675ba50-0946-4bc6-9daf-5da24e8eca08/phylophlan.tar?Si[...]<jupyter_text>Prepare input files<jupyter_code># IO dirs
input_dir = os.path.join(work_dir, "phylophlan")
phylophlan_out = os.path.join(input_dir, "trees")
if not os.path.exists(phylophlan_out):
    os.makedirs(phylophlan_out)
# Move faa files to input folder
# Methanomassilii proteins
included_genomes = os.path.join(work_dir, "prokka_output")
cp_cmd = "cp {0}/*.faa {1}"
cp_job = cp_cmd.format(included_genomes, input_dir)
print(cp_job)
!$cp_job

# Methanobrevibacter proteins
methanobrevibacter_genomes = os.path.join(work_dir, "genomes/Thermoplasmata")
cp_cmd = "cp {0}/*.faa {1}"
cp_job = cp_cmd.format(methanobrevibacter_genomes, input_dir)
print(cp_job)
!$cp_job<jupyter_output>cp /ebio/abt3_projects/vadinCA11/data/V11/prokka_output/*.faa /ebio/abt3_projects/vadinCA11/data/V11/phylophlan
cp /ebio/abt3_projects/vadinCA11/data/V11/genomes/Thermoplasmata/*.faa /ebio/abt3_projects/vadinCA11/data/V11/phylophlan<jupyter_text>Execute phylophlan<jupyter_code># Print version used
phylophlan_cmd = "python {0} --version".format(phylophlan_exe)
phylophlan_job = 'bash -c "source activate {0}; {1}"'.format(phylophlan_env, phylophlan_cmd)
print(phylophlan_job)
!$phylophlan_job
# Create template config file
config_base = os.path.join(phylophlan_dir, "supermatrix_base.cfg")
config_script = os.path.join(phylophlan_dir, "phylophlan_write_config_file.py")
config_cmd = "python {script} \
    -o {out_file} \
    -d a \
    --db_aa diamond \
    --map_dna blastn \
    --map_dna diamond \
    --map_aa diamond \
    --msa mafft \
    --trim trimal \
    --tree1 fasttree \
    --tree2 raxml \
    --overwrite \
    --verbose"
config_cmd = config_cmd.format(script = config_script, out_file = config_file)
default_configs = os.path.join(phylophlan_dir, "phylophlan_write_default_configs.sh")
config_job = 'bash -c "source activate {0}; {1}"'.format(phylophlan_env, config_cmd)
print(config_job)
!$config_job
# I copied the content of the base config file
# added the [map_dna], which is needed to run for some reason

config_text = """
[db_aa]
program_name = diamond
params = makedb
threads = --threads
input = --in
output = --db
version = version
command_line = #program_name# #params# #threads# #input# #output#

[map_dna]
program_name = diamond
params = blastx --quiet --threads 1 --outfmt 6 --more-sensitive --id 50 --max-hsps 35 -k 0
input = --query
database = --db
output = --out
version = version
command_line = #program_name# #params# #input# #database# #output#

[map_aa]
program_name = diamond
params = blastp --quiet --threads 1 --outfmt 6 --more-sensitive --id 50 --max-hsps 35 -k 0
input = --query
database = --db
output = --out
version = version
command_line = #program_name# #params# #input# #database# #output#

[msa]
program_name = mafft
params = --quiet --anysymbol --auto
version = --version
command_line = #program_name# #params# #input# > #output#
environment = TMPDIR=/tmp

[trim]
program_name = trimal
params = -gappyout
input = -in
output = -out
version = --version
command_line = #program_name# #params# #input# #output#

[tree1]
program_name = FastTreeMP
params = -quiet -mlacc 2 -slownni -spr 4 -fastest -mlnni 4 -no2nd -lg
output = -out
command_line = #program_name# #params# #output# #input#
environment = OMP_NUM_THREADS=3

[tree2]
program_name = raxmlHPC-PTHREADS-SSE3
params = -p 1989 -m PROTCATLG
database = -t
input = -s
output_path = -w
output = -n
version = -v
command_line = #program_name# #params# #threads# #database# #output_path# #input# #output#
threads = -T
"""
# Write config file
config_file = os.path.join(phylophlan_dir, "supermatrix_aa.cfg")
with open(config_file, 'w') as outF:
    outF.write(config_text)
# Execute phylophlan
# First time executing the folder .../phylophlan/phylophlan_databases/phylophlan
# has to be renamed to  .../phylophlan/phylophlan_databases/PhyloPhlAn
phylophlan_cmd = "python {script} \
    --input {in_folder} \
    --output {out_folder} \
    --database {database} \
    --db_type a \
    --diversity {div} \
    --config_file {config} \
    --min_num_markers {n_markers} \
    --nproc {cpu} \
    --accurate "
phylophlan_cmd = phylophlan_cmd.format(script = phylophlan_exe,
                                        in_folder = input_dir, 
                                        out_folder = phylophlan_out, 
                                        database = "PhyloPhlAn", 
                                        div = "high ", 
                                        config = config_file, 
                                        n_markers = 50,
                                        cpu = 20)

phylophlan_job = 'bash -c "source activate {0}; {1}"'.format(phylophlan_env, phylophlan_cmd)
print(phylophlan_job)
!$phylophlan_job<jupyter_output>bash -c "source activate py3_phylophlan; python /ebio/abt3_projects/vadinCA11/bin/phylophlan/phylophlan2.py     --input /ebio/abt3_projects/vadinCA11/data/V11/phylophlan     --output /ebio/abt3_projects/vadinCA11/data/V11/phylophlan/trees     --database PhyloPhlAn     --db_type a     --diversity high      --config_file /ebio/abt3_projects/vadinCA11/bin/phylophlan/supermatrix_aa.cfg     --min_num_markers 50     --nproc 20     --accurate "
[w] cannot create database "PhyloPhlAn", section "db_dna" not present in configurations
Loading files from "/ebio/abt3_projects/vadinCA11/data/V11/phylophlan"
Loading files from "/ebio/abt3_projects/vadinCA11/data/V11/phylophlan"
Checking 77 inputs
Checking "/ebio/abt3_projects/vadinCA11/data/V11/phylophlan/GCA_000300255.faa"
Checking "/ebio/abt3_projects/vadinCA11/data/V11/phylophlan/GCA_000308215.faa"
Checking "/ebio/abt3_projects/vadinCA11/data/V11/phylophlan/GCA_000404225.faa"
Checking "/ebio/abt3_projects/vadinCA11/data/V11/phylophlan/GCA_00080080[...]<jupyter_text>Running RAxML to get bootstraped tree<jupyter_code># phylophlan alignment
phylophlan_aln = os.path.join(phylophlan_out, "phylophlan_concatenated.aln")
# Execute RAxML
# Note that the fixed parameters are the same used by phylophlan
RAxML_cmd =  """raxmlHPC-PTHREADS-SSE3 \
    -f a \
    -x {seed} \
    -N {iterations} \
    -m PROTCATLG \
    -p 1989 \
    -T 20 \
    -w {workdir} \
    -s {aln} \
    -n {file_sufix}
"""

RAxML_cmd = RAxML_cmd.format(seed = 2112, 
                             iterations = 1000, 
                             workdir = phylophlan_out, 
                             aln = phylophlan_aln, 
                             file_sufix = "methanomassilii.tre")

RAxML_job = 'bash -c "source activate {0}; {1}"'.format(phylophlan_env, RAxML_cmd)
print(RAxML_job)
!$RAxML_job<jupyter_output>bash -c "source activate py3_phylophlan; raxmlHPC-PTHREADS-SSE3     -f a     -x 2112     -N 1000     -m PROTCATLG     -p 1989     -T 20     -w /ebio/abt3_projects/vadinCA11/data/V11/phylophlan/trees     -s /ebio/abt3_projects/vadinCA11/data/V11/phylophlan/trees/phylophlan_concatenated.aln     -n methanomassilii.tre
"
Warning, you specified a working directory via "-w"
Keep in mind that RAxML only accepts absolute path names, not relative ones!

RAxML can't, parse the alignment file as phylip file 
it will now try to parse it as FASTA file


This is the RAxML Master Pthread

This is RAxML Worker Pthread Number: 1

This is RAxML Worker Pthread Number: 2

This is RAxML Worker Pthread Number: 3

This is RAxML Worker Pthread Number: 10

This is RAxML Worker Pthread Number: 6

This is RAxML Worker Pthread Number: 4

This is RAxML Worker Pthread Number: 5

This is RAxML Worker Pthread Number: 7

This is RAxML Worker Pthread Number: 11

This is RAxML Worker Pthread Number: 12

This is RAxML Wo[...]<jupyter_text>The final tree with bootstrap values is in `/ebio/abt3_projects/vadinCA11/data/V11/phylophlan/trees/RAxML_bipartitions.methanomassilii.tre`  Session Info<jupyter_code>!conda list -n $phylophlan_env<jupyter_output># packages in environment at /ebio/abt3_projects/software/miniconda3_gt4.4/envs/py3_phylophlan:
#
# 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
biopython                 1.72             py36h04863e7_0    anaconda
blas                      1.1                    openblas    conda-forge
blast-plus                2.2.31                        0    biocore
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
bzip2                    [...]
{'fraction_non_alphanumeric': 0.10773809523809524, 'fraction_numerical': 0.04311224489795919, 'mean_word_length': 3.123772791023843, 'pattern_counts': {'":': 0, '<': 21, '<?xml version=': 0, '>': 22, 'https://': 9, 'lorem ipsum': 0, 'www.': 0, 'xml': 4}, 'pii_count': 4, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 21, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '6', 'dir': 'jupyter-structured-clean-dedup', 'id': '658565', 'license': 'MIT', 'n_tokens_mistral': 4730, 'n_tokens_neox': 4228, 'n_words': 894, 'path': 'notebooks/09_phylophlan.ipynb', 'repo_name': 'leylabmpi/Methanomassilii', 'size': '11773'}