workflow is updated

Emrah Akkoyun · Emrah Akkoyun · commit 2e226a40febf · 2021-06-02T13:17:32.000+03:00
diff --git a/config/config.yml b/config/config.yml
@@ -1,9 +1,8 @@
 cnfigfile: "config.yml"
 
 #outputs
-workdir: "/cta/users/eakkoyun/WORKFOLDER/TEST/120421_test/phylogeny-snakemake"
-query_ids: ["P35520", "P63000"]
-#query_ids: ["P16157", "P22105"]
+workdir: "/cta/users/eakkoyun/WORKFOLDER/TEST/020621_test/phylogeny-snakemake"
+query_ids: ["P13942", "P16157"]
 
 
 #blast
diff --git a/images/rulegraph.svg b/images/rulegraph.svg
diff --git a/workflow/Snakefile b/workflow/Snakefile
@@ -1,8 +1,8 @@
+
 configfile: "../config/config.yml"
 
 rule all:
     input:
-        expand("{workdir}/results/{query_id}/8_pruned_raxmlng_ancestral_scores/pruned_{query_id}_wol_param_{pattern}.csv",workdir=config["workdir"], query_id=config["query_ids"],pattern=config["pattern"]),
         expand("{workdir}/results/{query_id}/5_raxmlng_ancestral_scores/{query_id}_wl_param_{pattern}.csv",workdir=config["workdir"], query_id=config["query_ids"],pattern=config["pattern"]),
 
 include : "rules/query_fasta.smk"
@@ -12,11 +12,7 @@ include : "rules/msa.smk"
 include : "rules/trim_msa.smk"
 include : "rules/ml_tree.smk"
 include : "rules/remove_gaps.smk"
-include : "rules/remove_outliers.smk"
-include : "rules/prune.smk"
 include : "rules/unroot_tree.smk"
 include : "rules/raxmlng_ancestral.smk"
 include : "rules/compute_raxml_anc_score.smk"
-include : "rules/unroot_pruned_tree.smk"
-include : "rules/pruned_raxmlng_ancestral.smk"
-include : "rules/pruned_raxmlng_ancestral_score.smk"
+
diff --git a/workflow/envs/raxml-ng.yml b/workflow/envs/raxml-ng.yml
@@ -1,5 +1,7 @@
 name: raxml-ng 
 channels:
   - bioconda
+  - conda-forge 
 dependencies:
   - raxml-ng = 1.0.2 
+  - openmpi = 4.0.5 
diff --git a/workflow/rules/prune.smk b/workflow/rules/prune.smk
@@ -1,12 +1,12 @@
 rule prune:
     input:
         blastp_out_file = "{workdir}/results/{query_id}/1_psiblast/{query_id}_blasthits.out",
-        ml_tree = "{workdir}/results/{query_id}/3_mltree/{query_id}.no_outlier.nwk",
-        msa_file = "{workdir}/results/{query_id}/2_msa/{query_id}_no_outlier_no_gap.fasta",
+        ml_tree = "{workdir}/results/{query_id}/3_mltree/{query_id}.raxml.bestTree"
+        msa_file = "{workdir}/results/{query_id}/2_msa/{query_id}_nogap_msa.fasta",
         query_fasta = "{workdir}/results/{query_id}/1_psiblast/{query_id}.fasta",
     output:
-        pruned_tree =  "{workdir}/results/{query_id}/6_pruned_msa_tree/{query_id}.pruned.no_outlier",
-        pruned_msa = "{workdir}/results/{query_id}/6_pruned_msa_tree/{query_id}_pruned_nogap_nooutlier_msa.fasta",
+        pruned_tree =  "{workdir}/results/{query_id}/6_pruned_msa_tree/{query_id}.pruned",
+        pruned_msa = "{workdir}/results/{query_id}/6_pruned_msa_tree/{query_id}_pruned_nogap_msa.fasta",
     log:
         "{workdir}/workflow/logs/rules/{query_id}_pruned.err"
     benchmark:
@@ -17,3 +17,6 @@ rule prune:
     shell:
         "python3 scripts/prune_msa_tree.py {input.blastp_out_file} {config[blast_hit_number]} {config[max_e_value]} {config[min_identity]} {config[max_identity]} {input.ml_tree} {output.pruned_tree} {input.msa_file} {output.pruned_msa} {input.query_fasta}"
 
+
+
+
diff --git a/workflow/rules/pruned_raxmlng_ancestral.smk b/workflow/rules/pruned_raxmlng_ancestral.smk
@@ -1,8 +1,7 @@
-
 rule pruned_raxmlng_ancestral:
     input:
-        pruned_nogap_msa = "{workdir}/results/{query_id}/6_pruned_msa_tree/{query_id}_pruned_nogap_nooutlier_msa.fasta",
-        pruned_unrooted_tree = "{workdir}/results/{query_id}/6_pruned_msa_tree/{query_id}.pruned.no_outlier_unrooted",        
+        pruned_nogap_msa = "{workdir}/results/{query_id}/6_pruned_msa_tree/{query_id}_pruned_nogap_msa.fasta",
+        pruned_unrooted_tree = "{workdir}/results/{query_id}/6_pruned_msa_tree/{query_id}.pruned_unrooted",        
     output:
         pruned_ancestral_probabilities = "{workdir}/results/{query_id}/7_pruned_raxmlng_ancestral/pruned_{query_id}.raxml.ancestralProbs",
         pruned_ancestral_states = "{workdir}/results/{query_id}/7_pruned_raxmlng_ancestral/pruned_{query_id}.raxml.ancestralStates",
@@ -19,3 +18,5 @@ rule pruned_raxmlng_ancestral:
         "{workdir}/workflow/logs/benchmarks/{query_id}_pruned_raxmlng_ancestral.out"
     shell:
         "raxml-ng --ancestral --msa {input.pruned_nogap_msa} --tree {input.pruned_unrooted_tree} --model {config[raxmlng_ancestral_model]} --prefix {params.pruned_raxml_ancestral_out_name}  --threads {config[raxml_ancestral_threads]}"
+
+
diff --git a/workflow/rules/pruned_raxmlng_ancestral_score.smk b/workflow/rules/pruned_raxmlng_ancestral_score.smk
@@ -7,7 +7,7 @@ rule pruned_raxml_anc_score:
         "{workdir}/results/{query_id}/8_pruned_raxmlng_ancestral_scores/pruned_{query_id}_wl_param_{pattern}.csv",
     params:
         out = "{workdir}/results/{query_id}/8_pruned_raxmlng_ancestral_scores/pruned_{query_id}",
-        fasta = "{workdir}/results/{query_id}/6_pruned_msa_tree/{query_id}_pruned_nogap_nooutlier_msa.fasta",
+        fasta = "{workdir}/results/{query_id}/6_pruned_msa_tree/{query_id}_pruned_nogap_msa.fasta",
         query_fasta = "{workdir}/results/{query_id}/1_psiblast/{query_id}.fasta",
     log:
         "{workdir}/workflow/logs/rules/{query_id}_pruned_raxmlanc_{pattern}_compute_score.err"
@@ -19,3 +19,4 @@ rule pruned_raxml_anc_score:
     shell:
         "query=`python scripts/get_query.py {params.query_fasta}` && Rscript scripts/compute_score_RaxmlNg_Final.R {input.tree_file} {input.probabilities} {params.fasta} {params.out} $query {config[weights]} 2>{log}"
 
+
diff --git a/workflow/rules/raxmlng_ancestral.smk b/workflow/rules/raxmlng_ancestral.smk
@@ -1,7 +1,7 @@
 rule raxmlng_ancestral:
     input:
-        no_outlier_msa = "{workdir}/results/{query_id}/2_msa/{query_id}_no_outlier_no_gap.fasta",
-        unrooted_tree = "{workdir}/results/{query_id}/3_mltree/{query_id}.no_outlier.nwk_unrooted",        
+        no_gap_msa = "{workdir}/results/{query_id}/2_msa/{query_id}_nogap_msa.fasta",
+        unrooted_tree = "{workdir}/results/{query_id}/3_mltree/{query_id}.raxml.bestTree_unrooted",        
     output:
         ancestral_probabilities = "{workdir}/results/{query_id}/4_raxmlng_ancestral/{query_id}.raxml.ancestralProbs",
         ancestral_states = "{workdir}/results/{query_id}/4_raxmlng_ancestral/{query_id}.raxml.ancestralStates",
@@ -18,5 +18,4 @@ rule raxmlng_ancestral:
     resources:
         time_min=7200,cpus=4
     shell:
-        "raxml-ng --ancestral --msa {input.no_outlier_msa} --tree {input.unrooted_tree} --model {config[raxmlng_ancestral_model]} --prefix {params.raxml_ancestral_out_name} --threads {resources.cpus}"
-
+        "raxml-ng --ancestral --msa {input.no_gap_msa} --tree {input.unrooted_tree} --model {config[raxmlng_ancestral_model]} --prefix {params.raxml_ancestral_out_name} --threads {resources.cpus}"
diff --git a/workflow/rules/unroot_pruned_tree.smk b/workflow/rules/unroot_pruned_tree.smk
@@ -1,8 +1,8 @@
 rule unroot_pruned_tree:
     input:
-        input_tree = "{workdir}/results/{query_id}/6_pruned_msa_tree/{query_id}.pruned.no_outlier",
+        input_tree = "{workdir}/results/{query_id}/6_pruned_msa_tree/{query_id}.pruned",
     output:
-        "{workdir}/results/{query_id}/6_pruned_msa_tree/{query_id}.pruned.no_outlier_unrooted",
+        "{workdir}/results/{query_id}/6_pruned_msa_tree/{query_id}.pruned_unrooted",
     log:
         "{workdir}/workflow/logs/rules/{query_id}_pruned_unroot_tree.err"
     benchmark:
diff --git a/workflow/rules/unroot_tree.smk b/workflow/rules/unroot_tree.smk
@@ -1,8 +1,8 @@
 rule unroot_tree:
     input:
-        input_tree = "{workdir}/results/{query_id}/3_mltree/{query_id}.no_outlier.nwk",
+        input_tree = "{workdir}/results/{query_id}/3_mltree/{query_id}.raxml.bestTree"
     output:
-      "{workdir}/results/{query_id}/3_mltree/{query_id}.no_outlier.nwk_unrooted",
+      "{workdir}/results/{query_id}/3_mltree/{query_id}.raxml.bestTree_unrooted",
     log:
         "{workdir}/workflow/logs/rules/{query_id}_unroot_tree.err"
     benchmark:
diff --git a/workflow/scripts/raxml-ng.sh b/workflow/scripts/raxml-ng.sh
@@ -10,6 +10,12 @@ tree_num=$6
 # recommended # threads
 thread_num=`raxml-ng --parse --msa $trimmed_msa --model $raxml_model --nofile|grep "Recommended number of threads"|awk -F: {'print $2'}| xargs`
 echo "Recommended number of thread: $thread_num"
+echo "Number of CPUs on the allocated nodes: $SLURM_CPUS_ON_NODE" 
+
+if (("$thread_num" <= "$SLURM_CPUS_ON_NODE")); then
+	folder=`dirname $output_tree`
+	raxml-ng  --msa $trimmed_msa --data-type AA --model $raxml_model --prefix $folder/$raxml_out_name --seed $raxml_seed --threads $thread_num --tree rand{$tree_num} pars{$tree_num}
+else
+	echo "ERROR: Required CPU number: $thread_num. Assign more CPUs for the next submission"
+fi
 
-folder=`dirname $output_tree`
-raxml-ng  --msa $trimmed_msa --data-type AA --model $raxml_model --prefix $folder/$raxml_out_name --seed $raxml_seed --threads $thread_num --tree rand{$tree_num} pars{$tree_num}
diff --git a/workflow/scripts/remove_paralogs.py b/workflow/scripts/remove_paralogs.py
@@ -1,6 +1,4 @@
 import sys
-import os
-import re
 import ete3 
 from ete3 import Tree, PhyloTree
 from fasta_dict import *
