Mercurial > repos > davidvanzessen > shm_csr

annotate gene_identification.py @ 98:d714f5ea83d7 draft default tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload commit 1a01065a084a817382872154f779b94090a35ebf
author rhpvorderman
date Wed, 10 Jan 2024 12:32:47 +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()