Mercurial > repos > davidvanzessen > shm_csr

annotate gene_identification.py @ 87:64bbc6734ec7 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
"planemo upload commit 0dcc3e4e47a066373ff0bd56f16536298f8ac2a0"
author rhpvorderman
date Wed, 27 Oct 2021 11:30:25 +0000
parents 729738462297
children cf8ad181628f
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
81
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
1 import re
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
2 import argparse
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
3 import time
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
4 starttime= int(time.time() * 1000)
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
5
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
6 parser = argparse.ArgumentParser()
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
7 parser.add_argument("--input", help="The 1_Summary file from an IMGT zip file")
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
8 parser.add_argument("--output", help="The annotated output file to be merged back with the summary file")
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
9
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
10 args = parser.parse_args()
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
11
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
12 infile = args.input
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
13 #infile = "test_VH-Ca_Cg_25nt/1_Summary_test_VH-Ca_Cg_25nt_241013.txt"
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
14 output = args.output
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
15 #outfile = "identified.txt"
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
16
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
17 dic = dict()
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
18 total = 0
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
19
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
20
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
21 first = True
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
22 IDIndex = 0
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
23 seqIndex = 0
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
24
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
25 with open(infile, 'r') as f: #read all sequences into a dictionary as key = ID, value = sequence
83
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
26 for line in f:
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
27 total += 1
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
28 linesplt = line.split("\t")
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
29 if first:
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
30 print("linesplt", linesplt)
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
31 IDIndex = linesplt.index("Sequence ID")
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
32 seqIndex = linesplt.index("Sequence")
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
33 first = False
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
34 continue
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
35
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
36 ID = linesplt[IDIndex]
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
37 if len(linesplt) < 28: #weird rows without a sequence
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
38 dic[ID] = ""
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
39 else:
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
40 dic[ID] = linesplt[seqIndex]
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
41
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
42 print("Number of input sequences:", len(dic))
81
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
43
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
44 #old cm sequence: gggagtgcatccgccccaacccttttccccctcgtctcctgtgagaattccc
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
45 #old cg sequence: ctccaccaagggcccatcggtcttccccctggcaccctcctccaagagcacctctgggggcacagcggccctgggctgcctggtcaaggactacttccccgaaccggtgacggtgtcgtggaactcaggcgccctgaccag
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
46
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
47 #lambda/kappa reference sequence
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
48 searchstrings = {"ca": "catccccgaccagccccaaggtcttcccgctgagcctctgcagcacccagccagatgggaacgtggtcatcgcctgcctgg",
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
49 "cg": "ctccaccaagggcccatcggtcttccccctggcaccctcctccaagagcacctctgggggcacagcggcc",
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
50 "ce": "gcctccacacagagcccatccgtcttccccttgacccgctgctgcaaaaacattccctcc",
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
51 "cm": "gggagtgcatccgccccaacc"} #new (shorter) cm sequence
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
52
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
53 compiledregex = {"ca": [],
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
54 "cg": [],
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
55 "ce": [],
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
56 "cm": []}
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
57
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
58 #lambda/kappa reference sequence variable nucleotides
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
59 ca1 = {38: 't', 39: 'g', 48: 'a', 49: 'g', 51: 'c', 68: 'a', 73: 'c'}
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
60 ca2 = {38: 'g', 39: 'a', 48: 'c', 49: 'c', 51: 'a', 68: 'g', 73: 'a'}
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
61 cg1 = {0: 'c', 33: 'a', 38: 'c', 44: 'a', 54: 't', 56: 'g', 58: 'g', 66: 'g', 132: 'c'}
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
62 cg2 = {0: 'c', 33: 'g', 38: 'g', 44: 'g', 54: 'c', 56: 'a', 58: 'a', 66: 'g', 132: 't'}
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
63 cg3 = {0: 't', 33: 'g', 38: 'g', 44: 'g', 54: 't', 56: 'g', 58: 'g', 66: 'g', 132: 'c'}
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
64 cg4 = {0: 't', 33: 'g', 38: 'g', 44: 'g', 54: 'c', 56: 'a', 58: 'a', 66: 'c', 132: 'c'}
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
65
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
66 #remove last snp for shorter cg sequence --- note, also change varsInCG
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
67 del cg1[132]
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
68 del cg2[132]
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
69 del cg3[132]
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
70 del cg4[132]
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
71
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
72 #reference sequences are cut into smaller parts of 'chunklength' length, and with 'chunklength' / 2 overlap
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
73 chunklength = 8
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
74
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
75 #create the chunks of the reference sequence with regular expressions for the variable nucleotides
83
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
76 for i in range(0, len(searchstrings["ca"]) - chunklength, chunklength // 2):
81
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
77 pos = i
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
78 chunk = searchstrings["ca"][i:i+chunklength]
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
79 result = ""
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
80 varsInResult = 0
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
81 for c in chunk:
83
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
82 if pos in list(ca1.keys()):
81
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
83 varsInResult += 1
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
84 result += "[" + ca1[pos] + ca2[pos] + "]"
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
85 else:
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
86 result += c
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
87 pos += 1
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
88 compiledregex["ca"].append((re.compile(result), varsInResult))
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
89
83
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
90 for i in range(0, len(searchstrings["cg"]) - chunklength, chunklength // 2):
81
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
91 pos = i
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
92 chunk = searchstrings["cg"][i:i+chunklength]
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
93 result = ""
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
94 varsInResult = 0
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
95 for c in chunk:
83
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
96 if pos in list(cg1.keys()):
81
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
97 varsInResult += 1
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
98 result += "[" + "".join(set([cg1[pos], cg2[pos], cg3[pos], cg4[pos]])) + "]"
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
99 else:
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
100 result += c
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
101 pos += 1
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
102 compiledregex["cg"].append((re.compile(result), varsInResult))
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
103
83
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
104 for i in range(0, len(searchstrings["cm"]) - chunklength, chunklength // 2):
81
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
105 compiledregex["cm"].append((re.compile(searchstrings["cm"][i:i+chunklength]), False))
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
106
83
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
107 for i in range(0, len(searchstrings["ce"]) - chunklength + 1, chunklength // 2):
81
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
108 compiledregex["ce"].append((re.compile(searchstrings["ce"][i:i+chunklength]), False))
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
109
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
110 def removeAndReturnMaxIndex(x): #simplifies a list comprehension
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
111 m = max(x)
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
112 index = x.index(m)
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
113 x[index] = 0
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
114 return index
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
115
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
116
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
117 start_location = dict()
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
118 hits = dict()
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
119 alltotal = 0
83
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
120 for key in compiledregex: #for ca/cg/cm/ce
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
121 regularexpressions = compiledregex[key] # get the compiled regular expressions
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
122 for ID in list(dic.keys())[0:]: #for every ID
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
123 if ID not in list(hits.keys()): #ensure that the dictionairy that keeps track of the hits for every gene exists
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
124 hits[ID] = {"ca_hits": 0, "cg_hits": 0, "cm_hits": 0, "ce_hits": 0, "ca1": 0, "ca2": 0, "cg1": 0, "cg2": 0, "cg3": 0, "cg4": 0}
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
125 currentIDHits = hits[ID]
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
126 seq = dic[ID]
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
127 lastindex = 0
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
128 start_zero = len(searchstrings[key]) #allows the reference sequence to start before search sequence (start_locations of < 0)
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
129 start = [0] * (len(seq) + start_zero)
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
130 for i, regexp in enumerate(regularexpressions): #for every regular expression
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
131 relativeStartLocation = lastindex - (chunklength // 2) * i
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
132 if relativeStartLocation >= len(seq):
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
133 break
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
134 regex, hasVar = regexp
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
135 matches = regex.finditer(seq[lastindex:])
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
136 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
137 lastindex += match.start()
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
138 start[relativeStartLocation + start_zero] += 1
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
139 if hasVar: #if the regex has a variable nt in it
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
140 chunkstart = chunklength // 2 * i #where in the reference does this chunk start
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
141 chunkend = chunklength // 2 * i + chunklength #where in the reference does this chunk end
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
142 if key == "ca": #just calculate the variable nt score for 'ca', cheaper
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
143 currentIDHits["ca1"] += len([1 for x in ca1 if chunkstart <= x < chunkend and ca1[x] == seq[lastindex + x - chunkstart]])
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
144 currentIDHits["ca2"] += len([1 for x in ca2 if chunkstart <= x < chunkend and ca2[x] == seq[lastindex + x - chunkstart]])
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
145 elif key == "cg": #just calculate the variable nt score for 'cg', cheaper
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
146 currentIDHits["cg1"] += len([1 for x in cg1 if chunkstart <= x < chunkend and cg1[x] == seq[lastindex + x - chunkstart]])
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
147 currentIDHits["cg2"] += len([1 for x in cg2 if chunkstart <= x < chunkend and cg2[x] == seq[lastindex + x - chunkstart]])
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
148 currentIDHits["cg3"] += len([1 for x in cg3 if chunkstart <= x < chunkend and cg3[x] == seq[lastindex + x - chunkstart]])
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
149 currentIDHits["cg4"] += len([1 for x in cg4 if chunkstart <= x < chunkend and cg4[x] == seq[lastindex + x - chunkstart]])
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
150 else: #key == "cm" #no variable regions in 'cm' or 'ce'
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
151 pass
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
152 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
153 else: #only runs if there were no hits
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
154 continue
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
155 #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
156 currentIDHits[key + "_hits"] += 1
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
157 start_location[ID + "_" + key] = str([(removeAndReturnMaxIndex(start) + 1 - start_zero) for x in range(5) if len(start) > 0 and max(start) > 1])
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
158 #start_location[ID + "_" + key] = str(start.index(max(start)))
81
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
159
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
160
83
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
161 varsInCA = float(len(list(ca1.keys())) * 2)
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
162 varsInCG = float(len(list(cg1.keys())) * 2) - 2 # -2 because the sliding window doesn't hit the first and last nt twice
81
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
163 varsInCM = 0
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
164 varsInCE = 0
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
165
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
166 def round_int(val):
83
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
167 return int(round(val))
81
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
168
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
169 first = True
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
170 seq_write_count=0
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
171 with open(infile, 'r') as f: #read all sequences into a dictionary as key = ID, value = sequence
83
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
172 with open(output, 'w') as o:
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
173 for line in f:
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
174 total += 1
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
175 if first:
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
176 o.write("Sequence ID\tbest_match\tnt_hit_percentage\tchunk_hit_percentage\tstart_locations\n")
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
177 first = False
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
178 continue
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
179 linesplt = line.split("\t")
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
180 if linesplt[2] == "No results":
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
181 pass
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
182 ID = linesplt[1]
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
183 currentIDHits = hits[ID]
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
184 possibleca = float(len(compiledregex["ca"]))
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
185 possiblecg = float(len(compiledregex["cg"]))
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
186 possiblecm = float(len(compiledregex["cm"]))
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
187 possiblece = float(len(compiledregex["ce"]))
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
188 cahits = currentIDHits["ca_hits"]
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
189 cghits = currentIDHits["cg_hits"]
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
190 cmhits = currentIDHits["cm_hits"]
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
191 cehits = currentIDHits["ce_hits"]
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
192 if cahits >= cghits and cahits >= cmhits and cahits >= cehits: #its a ca gene
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
193 ca1hits = currentIDHits["ca1"]
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
194 ca2hits = currentIDHits["ca2"]
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
195 if ca1hits >= ca2hits:
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
196 o.write(ID + "\tIGA1\t" + str(round_int(ca1hits / varsInCA * 100)) + "\t" + str(round_int(cahits / possibleca * 100)) + "\t" + start_location[ID + "_ca"] + "\n")
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
197 else:
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
198 o.write(ID + "\tIGA2\t" + str(round_int(ca2hits / varsInCA * 100)) + "\t" + str(round_int(cahits / possibleca * 100)) + "\t" + start_location[ID + "_ca"] + "\n")
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
199 elif cghits >= cahits and cghits >= cmhits and cghits >= cehits: #its a cg gene
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
200 cg1hits = currentIDHits["cg1"]
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
201 cg2hits = currentIDHits["cg2"]
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
202 cg3hits = currentIDHits["cg3"]
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
203 cg4hits = currentIDHits["cg4"]
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
204 if cg1hits >= cg2hits and cg1hits >= cg3hits and cg1hits >= cg4hits: #cg1 gene
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
205 o.write(ID + "\tIGG1\t" + str(round_int(cg1hits / varsInCG * 100)) + "\t" + str(round_int(cghits / possiblecg * 100)) + "\t" + start_location[ID + "_cg"] + "\n")
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
206 elif cg2hits >= cg1hits and cg2hits >= cg3hits and cg2hits >= cg4hits: #cg2 gene
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
207 o.write(ID + "\tIGG2\t" + str(round_int(cg2hits / varsInCG * 100)) + "\t" + str(round_int(cghits / possiblecg * 100)) + "\t" + start_location[ID + "_cg"] + "\n")
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
208 elif cg3hits >= cg1hits and cg3hits >= cg2hits and cg3hits >= cg4hits: #cg3 gene
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
209 o.write(ID + "\tIGG3\t" + str(round_int(cg3hits / varsInCG * 100)) + "\t" + str(round_int(cghits / possiblecg * 100)) + "\t" + start_location[ID + "_cg"] + "\n")
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
210 else: #cg4 gene
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
211 o.write(ID + "\tIGG4\t" + str(round_int(cg4hits / varsInCG * 100)) + "\t" + str(round_int(cghits / possiblecg * 100)) + "\t" + start_location[ID + "_cg"] + "\n")
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
212 else: #its a cm or ce gene
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
213 if cmhits >= cehits:
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
214 o.write(ID + "\tIGM\t100\t" + str(round_int(cmhits / possiblecm * 100)) + "\t" + start_location[ID + "_cm"] + "\n")
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
215 else:
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
216 o.write(ID + "\tIGE\t100\t" + str(round_int(cehits / possiblece * 100)) + "\t" + start_location[ID + "_ce"] + "\n")
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
217 seq_write_count += 1
81
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
218
83
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
219 print("Time: %i" % (int(time.time() * 1000) - starttime))
81
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
220
83
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
221 print("Number of sequences written to file:", seq_write_count)
81
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
222
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
223
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
224
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
225
b6f9a640e098 Uploaded
davidvanzessen
parents:
diff changeset
226