Mercurial > repos > davidvanzessen > shm_csr

annotate gene_identification.py @ 94:84e9e5c8c101 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
"planemo upload commit d4be85014b638f1d50b318d4b735be7f6e973140"
author rhpvorderman
date Fri, 24 Mar 2023 16:58:28 +0000
parents cf8ad181628f
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
92
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
1 #!/usr/bin/env python3
81
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
2
92
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
3 import argparse
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
4 import re
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
5 from typing import Dict, Iterator, List, Tuple
81
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
6
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
7
92
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
8 def generate_sequence_and_id_from_summary(summary_file: str
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
9 ) -> Iterator[Tuple[str, str]]:
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
10 with open(summary_file, "rt") as summary:
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
11 header = next(summary)
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
12 column_names = header.strip("\n").split("\t")
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
13 id_column = column_names.index("Sequence ID")
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
14 sequence_column = column_names.index("Sequence")
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
15 for line in summary:
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
16 values = line.strip("\n").split("\t")
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
17 id = values[id_column]
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
18 try:
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
19 sequence = values[sequence_column]
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
20 except IndexError: # weird rows without a sequence
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
21 sequence = ""
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
22 yield id, sequence
81
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
23
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
24
92
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
25 # old cm sequence: gggagtgcatccgccccaacccttttccccctcgtctcctgtgagaattccc
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
26 # old cg sequence: ctccaccaagggcccatcggtcttccccctggcaccctcctccaagagcacctctg
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
27 # ggggcacagcggccctgggctgcctggtcaaggactacttccccgaaccggtgacggtgtcgtggaactcagg
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
28 # cgccctgaccag
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
29 SEARCHSTRINGS = {"ca": "catccccgaccagccccaaggtcttcccgctgagcctctgcagcacccagccag"
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
30 "atgggaacgtggtcatcgcctgcctgg",
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
31 "cg": "ctccaccaagggcccatcggtcttccccctggcaccctcctccaagagcacctc"
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
32 "tgggggcacagcggcc",
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
33 "ce": "gcctccacacagagcccatccgtcttccccttgacccgctgctgcaaaaacatt"
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
34 "ccctcc",
81
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
35 "cm": "gggagtgcatccgccccaacc"} #new (shorter) cm sequence
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
36
92
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
37 #lambda/kappa referesearchstringsnce sequence variable nucleotides
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
38 CA1_MUTATIONS = {38: 't', 39: 'g', 48: 'a', 49: 'g', 51: 'c', 68: 'a', 73: 'c'}
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
39 CA2_MUTATIONS = {38: 'g', 39: 'a', 48: 'c', 49: 'c', 51: 'a', 68: 'g', 73: 'a'}
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
40 CG1_MUTATIONS = {0: 'c', 33: 'a', 38: 'c', 44: 'a', 54: 't', 56: 'g', 58: 'g', 66: 'g', 132: 'c'}
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
41 CG2_MUTATIONS = {0: 'c', 33: 'g', 38: 'g', 44: 'g', 54: 'c', 56: 'a', 58: 'a', 66: 'g', 132: 't'}
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
42 CG3_MUTATIONS = {0: 't', 33: 'g', 38: 'g', 44: 'g', 54: 't', 56: 'g', 58: 'g', 66: 'g', 132: 'c'}
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
43 CG4_MUTATIONS = {0: 't', 33: 'g', 38: 'g', 44: 'g', 54: 'c', 56: 'a', 58: 'a', 66: 'c', 132: 'c'}
81
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
44
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
45 #remove last snp for shorter cg sequence --- note, also change varsInCG
92
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
46 del CG1_MUTATIONS[132]
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
47 del CG2_MUTATIONS[132]
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
48 del CG3_MUTATIONS[132]
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
49 del CG4_MUTATIONS[132]
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
50
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
51 # reference sequences are cut into smaller parts of 'chunklength' length,
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
52 # and with 'chunklength' / 2 overlap
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
53 CHUNK_LENGTH = 8
81
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
54
92
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
55
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
56 def create_compiled_regexes() -> Dict[str, List[Tuple[re.Pattern, int]]]:
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
57
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
58 compiledregex: Dict[str, List[Tuple[re.Pattern, int]]] = {
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
59 "ca": [],
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
60 "cg": [],
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
61 "ce": [],
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
62 "cm": []
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
63 }
81
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
64
92
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
65 for i in range(0, len(SEARCHSTRINGS["ca"]) - CHUNK_LENGTH, CHUNK_LENGTH // 2):
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
66 pos = i
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
67 chunk = SEARCHSTRINGS["ca"][i:i + CHUNK_LENGTH]
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
68 result = ""
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
69 varsInResult = 0
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
70 for c in chunk:
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
71 if pos in list(CA1_MUTATIONS.keys()):
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
72 varsInResult += 1
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
73 result += "[" + CA1_MUTATIONS[pos] + CA2_MUTATIONS[pos] + "]"
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
74 else:
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
75 result += c
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
76 pos += 1
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
77 compiledregex["ca"].append((re.compile(result), varsInResult))
81
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
78
92
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
79 for i in range(0, len(SEARCHSTRINGS["cg"]) - CHUNK_LENGTH, CHUNK_LENGTH // 2):
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
80 pos = i
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
81 chunk = SEARCHSTRINGS["cg"][i:i + CHUNK_LENGTH]
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
82 result = ""
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
83 varsInResult = 0
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
84 for c in chunk:
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
85 if pos in list(CG1_MUTATIONS.keys()):
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
86 varsInResult += 1
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
87 result += "[" + "".join(set([CG1_MUTATIONS[pos], CG2_MUTATIONS[pos], CG3_MUTATIONS[pos], CG4_MUTATIONS[pos]])) + "]"
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
88 else:
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
89 result += c
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
90 pos += 1
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
91 compiledregex["cg"].append((re.compile(result), varsInResult))
81
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
92
92
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
93 for i in range(0, len(SEARCHSTRINGS["cm"]) - CHUNK_LENGTH, CHUNK_LENGTH // 2):
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
94 compiledregex["cm"].append((re.compile(SEARCHSTRINGS["cm"][i:i + CHUNK_LENGTH]), 0))
81
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
95
92
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
96 for i in range(0, len(SEARCHSTRINGS["ce"]) - CHUNK_LENGTH + 1, CHUNK_LENGTH // 2):
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
97 compiledregex["ce"].append((re.compile(SEARCHSTRINGS["ce"][i:i + CHUNK_LENGTH]), 0))
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
98
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
99 return compiledregex
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
100
81
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
101
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
102 def removeAndReturnMaxIndex(x): #simplifies a list comprehension
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
103 m = max(x)
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
104 index = x.index(m)
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
105 x[index] = 0
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
106 return index
92
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
107
81
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
108
92
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
109 def match_sequence(seq, compiledregex):
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
110 currentIDHits = {"ca_hits": 0, "cg_hits": 0, "cm_hits": 0, "ce_hits": 0,
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
111 "ca1": 0, "ca2": 0, "cg1": 0, "cg2": 0, "cg3": 0, "cg4": 0}
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
112 alltotal = 0
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
113 start_location = dict()
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
114 for key in compiledregex: # for ca/cg/cm/ce
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
115 regularexpressions = compiledregex[key]
83
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
116 lastindex = 0
92
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
117 start_zero = len(SEARCHSTRINGS[key]) #allows the reference sequence to start before search sequence (start_locations of < 0)
83
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
118 start = [0] * (len(seq) + start_zero)
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
119 for i, regexp in enumerate(regularexpressions): #for every regular expression
92
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
120 relativeStartLocation = lastindex - (CHUNK_LENGTH // 2) * i
83
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
121 if relativeStartLocation >= len(seq):
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
122 break
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
123 regex, hasVar = regexp
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
124 matches = regex.finditer(seq[lastindex:])
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
125 for match in matches: #for every match with the current regex, only uses the first hit because of the break at the end of this loop
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
126 lastindex += match.start()
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
127 start[relativeStartLocation + start_zero] += 1
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
128 if hasVar: #if the regex has a variable nt in it
92
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
129 chunkstart = CHUNK_LENGTH // 2 * i #where in the reference does this chunk start
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
130 chunkend = CHUNK_LENGTH // 2 * i + CHUNK_LENGTH #where in the reference does this chunk end
83
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
131 if key == "ca": #just calculate the variable nt score for 'ca', cheaper
92
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
132 currentIDHits["ca1"] += len([1 for x in CA1_MUTATIONS if chunkstart <= x < chunkend and CA1_MUTATIONS[x] == seq[lastindex + x - chunkstart]])
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
133 currentIDHits["ca2"] += len([1 for x in CA2_MUTATIONS if chunkstart <= x < chunkend and CA2_MUTATIONS[x] == seq[lastindex + x - chunkstart]])
83
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
134 elif key == "cg": #just calculate the variable nt score for 'cg', cheaper
92
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
135 currentIDHits["cg1"] += len([1 for x in CG1_MUTATIONS if chunkstart <= x < chunkend and CG1_MUTATIONS[x] == seq[lastindex + x - chunkstart]])
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
136 currentIDHits["cg2"] += len([1 for x in CG2_MUTATIONS if chunkstart <= x < chunkend and CG2_MUTATIONS[x] == seq[lastindex + x - chunkstart]])
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
137 currentIDHits["cg3"] += len([1 for x in CG3_MUTATIONS if chunkstart <= x < chunkend and CG3_MUTATIONS[x] == seq[lastindex + x - chunkstart]])
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
138 currentIDHits["cg4"] += len([1 for x in CG4_MUTATIONS if chunkstart <= x < chunkend and CG4_MUTATIONS[x] == seq[lastindex + x - chunkstart]])
83
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
139 else: #key == "cm" #no variable regions in 'cm' or 'ce'
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
140 pass
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
141 break #this only breaks when there was a match with the regex, breaking means the 'else:' clause is skipped
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
142 else: #only runs if there were no hits
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
143 continue
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
144 #print "found ", regex.pattern , "at", lastindex, "adding one to", (lastindex - chunklength / 2 * i), "to the start array of", ID, "gene", key, "it's now:", start[lastindex - chunklength / 2 * i]
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
145 currentIDHits[key + "_hits"] += 1
92
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
146 start_location[key] = str([(removeAndReturnMaxIndex(start) + 1 - start_zero) for x in range(5) if len(start) > 0 and max(start) > 1])
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
147
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
148 cahits = currentIDHits["ca_hits"]
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
149 cghits = currentIDHits["cg_hits"]
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
150 cmhits = currentIDHits["cm_hits"]
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
151 cehits = currentIDHits["ce_hits"]
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
152 if cahits >= cghits and cahits >= cmhits and cahits >= cehits: # its a ca gene
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
153 ca1hits = currentIDHits["ca1"]
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
154 ca2hits = currentIDHits["ca2"]
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
155 if ca1hits >= ca2hits:
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
156 # TODO: All variants with 0 matched are matched to IGA1 with 0 hits
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
157 # TODO: these are later turned into unmatched by the merge_and_filter.R
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
158 # TODO: script
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
159 return "IGA1", ca1hits, cahits, start_location["ca"]
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
160 else:
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
161 return "IGA2", ca2hits, cahits, start_location["ca"]
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
162 elif cghits >= cahits and cghits >= cmhits and cghits >= cehits: # its a cg gene
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
163 cg1hits = currentIDHits["cg1"]
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
164 cg2hits = currentIDHits["cg2"]
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
165 cg3hits = currentIDHits["cg3"]
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
166 cg4hits = currentIDHits["cg4"]
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
167 if cg1hits >= cg2hits and cg1hits >= cg3hits and cg1hits >= cg4hits: # cg1 gene
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
168 return "IGG1", cg1hits, cghits, start_location["cg"]
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
169 elif cg2hits >= cg1hits and cg2hits >= cg3hits and cg2hits >= cg4hits: # cg2 gene
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
170 return "IGG2", cg2hits, cghits, start_location["cg"]
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
171 elif cg3hits >= cg1hits and cg3hits >= cg2hits and cg3hits >= cg4hits: # cg3 gene
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
172 return "IGG3", cg3hits, cghits, start_location["cg"]
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
173 else: # cg4 gene
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
174 return "IGG4", cg4hits, cghits, start_location["cg"]
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
175 else: # its a cm or ce gene
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
176 if cmhits >= cehits:
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
177 return "IGM", 0, cmhits, start_location["cm"]
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
178 else:
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
179 return "IGE", 0, cehits, start_location["ce"]
81
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
180
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
181
92
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
182 def main():
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
183 parser = argparse.ArgumentParser()
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
184 parser.add_argument("--input",
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
185 help="The 1_Summary file from an IMGT zip file")
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
186 parser.add_argument("--output",
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
187 help="The annotated output file to be merged back "
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
188 "with the summary file")
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
189 args = parser.parse_args()
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
190 varsInCA = float(len(list(CA1_MUTATIONS.keys())) * 2)
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
191 varsInCG = float(len(list(
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
192 CG1_MUTATIONS.keys())) * 2) - 2 # -2 because the sliding window doesn't hit the first and last nt twice
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
193 subclass_vars = {
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
194 "IGA1": varsInCA, "IGA2": varsInCA,
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
195 "IGG1": varsInCG, "IGG2": varsInCG, "IGG3": varsInCG, "IGG4": varsInCG,
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
196 "IGE": 0,
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
197 "IGM": 0,
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
198 }
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
199 compiledregex = create_compiled_regexes()
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
200 possibleca = float(len(compiledregex["ca"]))
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
201 possiblecg = float(len(compiledregex["cg"]))
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
202 possiblecm = float(len(compiledregex["cm"]))
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
203 possiblece = float(len(compiledregex["ce"]))
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
204 class_chunks = {
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
205 "IGA1": possibleca, "IGA2": possibleca,
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
206 "IGE": possiblece,
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
207 "IGG1": possiblecg, "IGG2": possiblecg, "IGG3": possiblecg,
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
208 "IGG4": possiblecg,
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
209 "IGM": possiblecm
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
210 }
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
211 with open(args.output, "wt") as output:
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
212 output.write("Sequence ID\tbest_match\tnt_hit_percentage\t"
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
213 "chunk_hit_percentage\tstart_locations\n")
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
214 for id, sequence in generate_sequence_and_id_from_summary(args.input):
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
215 best_match, subclass_hits, class_hits, start_locations = \
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
216 match_sequence(sequence, compiledregex)
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
217 variable_nucs = subclass_vars[best_match]
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
218 if variable_nucs:
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
219 subclass_percentage = round(subclass_hits * 100 /
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
220 variable_nucs)
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
221 else:
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
222 subclass_percentage = 100
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
223 class_percentage = round(class_hits * 100 / class_chunks[best_match])
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
224 output.write(f"{id}\t{best_match}\t{subclass_percentage}\t"
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
225 f"{class_percentage}\t{start_locations}\n")
81
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
226
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
227
92
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
228 if __name__ == "__main__":
cf8ad181628f planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
rhpvorderman
parents: 83
diff changeset
229 main()