Modification for experimentations on meta-consensus threshold

src/exonerate_stats.sh

 #!/bin/bash
 
-printf "output_file, percent_identity, percent_similarity, total_equivalence, total_mismatches\n" > $1
+printf "" > $1
 for FILE in $@ ;
 do
   if [ $FILE != $1 ] && [ $FILE != $2 ] ;
   then
-  printf " %q :" $FILE  >> $1
+  printf "%q," $FILE  >> $1
     exonerate --bestn 1 -Q dna -E -m a:g --showalignment false --showsugar false --showvulgar false --showcigar false  --ryo "%pi, %ps, %et, %em\n"  --verbose 0 -q $2 -t $FILE >> $1 2>/dev/null ;
   fi;
 done

src/plot_results.py

+import matplotlib as mpl
+import matplotlib.pyplot as plt
+import numpy as np
+import re
+import sys
+import os
+
+
+def parse_file(filename, x_scale):
+    identity, similarity, equivalence, mismatches = [ [], [] , [], []]
+    value = []
+    with open(filename, 'r') as f:
+        cnt= 0
+        for line in f:
+            i, s, e, m = map(float, line.split(',')[1:])
+            identity.append(i)
+            similarity.append(s)
+            equivalence.append(e)
+            mismatches.append(m)
+            cnt+=1
+    return identity, similarity, equivalence, mismatches
+
+
+th = [30, 40, 50, 60, 70, 80]
+ls, li, le, lm = [], [], [], []
+for tool_file in os.scandir(sys.argv[1]):
+    if ".txt" in tool_file.name:
+        i, s, e, m = map(np.array, parse_file(tool_file, th))
+        li.append(i)
+        ls.append(s)
+        le.append(e)
+        lm.append(m)
+
+
+ll = [ls, li, le, lm]
+cnt = 0
+target_att = ["identity", "similarity", "equivalence", "mismatches"]
+for graph in ll:
+    plt.clf()
+    for tab in graph:
+        if (len(tab) > 1):
+            plt.plot(th, tab, 'bo', label=cnt)
+        else:
+            plt.plot(th,[tab for x in range(len(th))])
+    plt.savefig(sys.argv[2]+"_" + target_att[cnt] + ".png")
+    cnt+=1

src/workflow/Snakefile

@@ -70,10 +70,31 @@ rule exonerate_ref_result:
           files=expand(rules.final_consensus.output, align="muscle", allow_missing=True) ,
           ref=expand("{data}/seq_selectes_region/region_seq_r{region}.fasta", data=DATA, allow_missing=True),
           consensus_per_tool=expand(rules.create_consensus_per_tool.output, tool=TOOL, allow_missing=True)
-     output: "consensus/stats/{region}/aln_all_consensus_r{region}_d{depth}_t1_{threshold}_t2_{value}.txt"
+     output: "consensus/stats/{region}/aln_all_consensus_r{region}_d{depth}_t1{threshold}_t2_{value}.txt"
      resources: format="\"%pi, %ps, %et, %em\n\""
      conda: "envs/exonerate.yaml"
      shell: "./exonerate_stats.sh {output} {input.ref} {input.files} {input.consensus_per_tool}"
 #          "echo {STAT_HEADER} >> {output} \n for FILE in {input.files}; do exonerate --bestn 1 -Q dna -E -m a:g --showalignment false --showsugar false --showvulgar false --showcigar false  --ryo {resources.format} --verbose 0 -q {input.ref} -t $FILE >> {output}  2>/dev/null; done" #
 
 
+
+
+rule plottable_exonerate_meta_t2:
+     input: files=expand(rules.final_consensus.output, align="muscle", value=THRESHOLD, allow_missing=True) , ref=expand("{data}/seq_selectes_region/region_seq_r{region}.fasta", data=DATA, allow_missing=True)
+     output: "consensus/plot/{region}/{depth}/{threshold}/aln_stats_meta.txt"
+     conda: "envs/exonerate.yaml"
+     shell: "./exonerate_stats.sh {output} {input.ref} {input.files}"
+
+rule plottable_exonerate_tools_t2:
+     input:
+          ref=expand("{data}/seq_selectes_region/region_seq_r{region}.fasta", data=DATA, allow_missing=True),
+          consensus_per_tool=expand(rules.create_consensus_per_tool.output, value=THRESHOLD, allow_missing=True)
+     output: "consensus/plot/{region}/{depth}/{threshold}/aln_stats_{tool}.txt"
+     conda: "envs/exonerate.yaml"
+     shell: "./exonerate_stats.sh {output} {input.ref} {input.consensus_per_tool}"
+
+rule plot_t2:
+     input: expand(rules.plottable_exonerate_tools_t2.output, region=REGION, depth=DEPTH, threshold=THRESHOLD, tool=TOOL),
+     	     expand(rules.plottable_exonerate_meta_t2.output, region=REGION, depth=DEPTH, threshold=THRESHOLD, tool=TOOL)
+
+