Mercurial > repos > sanbi-uwc > assess_poliovirus_alignment

annotate assess_alignment.py @ 4:1897677e107c draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
author sanbi-uwc
date Fri, 22 Jul 2022 09:55:46 +0000
parents 31ca16290d4f
children 0e556a3f85d6
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
0
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
1 #!/usr/bin/env python
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
2
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
3 import argparse
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
4 import json
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
5 import sys
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
6 from dataclasses import dataclass
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
7 from typing import TextIO
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
8
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
9
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
10 @dataclass
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
11 class Sample:
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
12 name: str
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
13 reference: str
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
14 mismatches: int
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
15 quality: str
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
16
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
17
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
18 def load_json(json_file: TextIO) -> dict:
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
19 data = json.load(json_file)
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
20 if "msa" not in data:
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
21 raise ValueError("MSA missing from JSON, cannot proceed")
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
22 if "gappedTraces" not in data:
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
23 raise ValueError("gappedTraces missing from JSON, cannot proceed")
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
24 return data
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
25
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
26
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
27 def analyse_mismatches(
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
28 json_file: TextIO,
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
29 offset: int,
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
30 length: int,
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
31 vp1only: bool = True,
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
32 sec_is_conflict: bool = False,
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
33 ) -> list:
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
34 data = load_json(json_file)
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
35 msas = [al for al in data["msa"] if not al["reference"]]
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
36 reference = None
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
37 for al in data["msa"]:
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
38 if al["reference"]:
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
39 if reference is None:
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
40 reference = al
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
41 else:
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
42 sys.exit(
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
43 "more than one reference found in JSON MSA list, cannot proceed"
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
44 )
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
45 min_start = min([int(al["leadingGaps"]) for al in msas])
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
46 max_end = max([int(al["leadingGaps"]) + len(al["align"]) for al in msas])
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
47 base_state = ["n"] * len(reference["align"])
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
48 mismatch_bases = {}
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
49 for i, base in enumerate(reference["align"]):
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
50 for k, al in enumerate(msas):
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
51 leading_gaps = int(al["leadingGaps"])
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
52 align_len = len(al["align"])
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
53 if leading_gaps < i and (leading_gaps + align_len) > i:
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
54 vp1pos = i - offset
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
55 if vp1only and vp1pos < 0 or vp1pos > length:
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
56 # skip positions outside of vp1 gene region
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
57 continue
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
58 al_base = al["align"][i - leading_gaps]
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
59 has_secondary_basecall = False
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
60 if sec_is_conflict:
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
61 gappedTrace = data["gappedTraces"][k]
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
62 pos = i - int(gappedTrace["leadingGaps"])
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
63 # print(len(gappedTrace['basecallPos']), pos, k, len(gappedTrace['basecalls']), gappedTrace['basecallPos'][pos])
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
64 basecall_str = gappedTrace["basecalls"][
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
65 str(gappedTrace["basecallPos"][pos])
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
66 ]
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
67 if "|" in basecall_str:
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
68 has_secondary_basecall = True
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
69 # set this position to conflicted
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
70 base_state[i] = "C"
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
71 if al_base != base:
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
72 # let's deal with all the cases where the base state doesn't match the reference
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
73 if base_state[i] == "G":
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
74 # the base state was G (a trace matches reference) and now we see a mismatch
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
75 base_state[i] = "C"
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
76 elif base_state[i] == "C":
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
77 # already marked as conflicting - a mismatch doesn't change that
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
78 pass
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
79 elif base_state[i] == "n" or base_state[i] == "M":
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
80 # we never saw this before or its already marked as a mismatch
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
81 base_state[i] = "M"
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
82 mismatch_bases[i] = al_base
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
83 else:
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
84 sys.exit("unexpected base state: " + base_state[i])
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
85 else:
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
86 if base_state[i] == "G" or base_state[i] == "n":
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
87 # we saw this before and got a match or
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
88 # we never saw this before
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
89 base_state[i] = "G"
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
90 elif base_state[i] == "M":
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
91 # we saw this before but it was a mismatch - mark this as a conflict
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
92 base_state[i] = "C"
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
93 if i in mismatch_bases:
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
94 del mismatch_bases[i]
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
95 elif base_state[i] == "C":
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
96 # we have seen a conflict here before
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
97 pass
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
98 else:
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
99 sys.exit("unexpected base_state: " + base_state[i])
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
100 conflicts = base_state.count("C")
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
101 matches = base_state.count("G")
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
102 mismatches = base_state.count("M")
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
103 mismatch_list = []
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
104 for i, state in enumerate(base_state):
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
105 # i is in zero-based genome coordinates
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
106 if state == "M":
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
107 # for mismatch store [pos_in_genome, pos_in_vp1, reference_base, sequenced_base]
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
108 mismatch_list.append(
4
1897677e107c planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents: 2
diff changeset
109 [i + 1, i - offset + 1, reference["align"][i], mismatch_bases[i]]
0
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
110 )
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
111 return [conflicts, matches, mismatches, mismatch_list]
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
112
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
113
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
114 def analyse_trace_quality(json_file: TextIO) -> float:
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
115 data = load_json(json_file)
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
116
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
117 traces = data["gappedTraces"]
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
118 overall_avg = 0
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
119 for trace in traces:
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
120 start = min(trace["basecallPos"])
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
121 end = max(trace["basecallPos"])
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
122 call_quality = {}
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
123 avg_ratio = 0
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
124 for base in ("A", "C", "G", "T"):
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
125 calls = trace["peak" + base][start : end + 1]
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
126 min_call = min(calls)
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
127 max_call = max(calls)
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
128 avg_call = sum(calls) / len(calls)
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
129 ratio = max_call / avg_call
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
130 call_quality["avg" + base] = avg_call
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
131 call_quality["min" + base] = min_call
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
132 call_quality["max" + base] = max_call
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
133 call_quality["ratio" + base] = ratio
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
134 avg_ratio += ratio
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
135 avg_ratio = avg_ratio / 4
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
136 overall_avg += avg_ratio
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
137 overall_avg = overall_avg / len(traces)
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
138 return overall_avg
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
139
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
140
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
141 def comma_split(args: str) -> list[str]:
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
142 return args.split(",")
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
143
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
144
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
145 if __name__ == "__main__":
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
146 parser = argparse.ArgumentParser()
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
147 parser.add_argument("--output_filename", help="Path to output file")
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
148 parser.add_argument("--sample_name", help="Name of sample being analysed")
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
149 parser.add_argument(
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
150 "--dataset_names", type=comma_split, help="Comma separated names for datasets"
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
151 )
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
152 parser.add_argument("--datasets", nargs="+")
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
153 args = parser.parse_args()
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
154
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
155 offsets = {
2
31ca16290d4f planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents: 0
diff changeset
156 # these are in 0-based coordinates, so off-by-one from NCBI 1-based coordinates
31ca16290d4f planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents: 0
diff changeset
157 "poliovirus1sabin": 2479, # V01150
31ca16290d4f planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents: 0
diff changeset
158 "poliovirus2sabin": 2481, # AY184220
31ca16290d4f planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents: 0
diff changeset
159 "poliovirus3sabin": 2478, # X00925
0
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
160 }
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
161
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
162 lengths = {
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
163 "poliovirus1sabin": 906,
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
164 "poliovirus2sabin": 903,
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
165 "poliovirus3sabin": 900,
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
166 }
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
167
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
168 min_mismatches = None
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
169 for file_index, json_filename in enumerate(args.datasets):
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
170 dataset_name = args.dataset_names[file_index].replace(
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
171 ".json", ""
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
172 ) # take the name but remove any json suffix
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
173 offset = offsets[dataset_name]
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
174 length = lengths[dataset_name]
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
175 (conflicts, matches, mismatches, mismatch_list) = analyse_mismatches(
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
176 open(json_filename), offset, length
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
177 )
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
178 # analyse_mismatches(json_filename, True)
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
179 quality = analyse_trace_quality(open(json_filename))
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
180 if min_mismatches is None or mismatches < min_mismatches:
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
181 min_mismatches = mismatches
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
182 best_match_mismatch_list = mismatch_list
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
183 best_match_quality = quality
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
184 best_match_reference = dataset_name
2
31ca16290d4f planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents: 0
diff changeset
185 percent_mismatches = round(min_mismatches / lengths[best_match_reference] * 100, 2)
0
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
186
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
187 info = {
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
188 "sample_name": args.sample_name,
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
189 "best_reference": best_match_reference,
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
190 "mismatches": min_mismatches,
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
191 "mismatch_list": best_match_mismatch_list,
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
192 "quality": best_match_quality,
2
31ca16290d4f planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents: 0
diff changeset
193 "perc_mismatches": percent_mismatches
0
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
194 }
8bc106442b1a planemo upload for repository https://github.com/pvanheus/polio_report commit a99e10fec2fac5aae70974c977eb3b362a1a8429-dirty
sanbi-uwc
parents:
diff changeset
195 json.dump(info, open(args.output_filename, "w"))