Mercurial > repos > davidvanzessen > shm_csr
diff gene_identification.py @ 83:729738462297 draft
"planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
author | rhpvorderman |
---|---|
date | Wed, 15 Sep 2021 12:24:06 +0000 |
parents | b6f9a640e098 |
children | cf8ad181628f |
line wrap: on
line diff
--- a/gene_identification.py Thu Feb 25 10:32:32 2021 +0000 +++ b/gene_identification.py Wed Sep 15 12:24:06 2021 +0000 @@ -23,23 +23,23 @@ seqIndex = 0 with open(infile, 'r') as f: #read all sequences into a dictionary as key = ID, value = sequence - for line in f: - total += 1 - linesplt = line.split("\t") - if first: - print "linesplt", linesplt - IDIndex = linesplt.index("Sequence ID") - seqIndex = linesplt.index("Sequence") - first = False - continue - - ID = linesplt[IDIndex] - if len(linesplt) < 28: #weird rows without a sequence - dic[ID] = "" - else: - dic[ID] = linesplt[seqIndex] - -print "Number of input sequences:", len(dic) + for line in f: + total += 1 + linesplt = line.split("\t") + if first: + print("linesplt", linesplt) + IDIndex = linesplt.index("Sequence ID") + seqIndex = linesplt.index("Sequence") + first = False + continue + + ID = linesplt[IDIndex] + if len(linesplt) < 28: #weird rows without a sequence + dic[ID] = "" + else: + dic[ID] = linesplt[seqIndex] + +print("Number of input sequences:", len(dic)) #old cm sequence: gggagtgcatccgccccaacccttttccccctcgtctcctgtgagaattccc #old cg sequence: ctccaccaagggcccatcggtcttccccctggcaccctcctccaagagcacctctgggggcacagcggccctgggctgcctggtcaaggactacttccccgaaccggtgacggtgtcgtggaactcaggcgccctgaccag @@ -73,13 +73,13 @@ chunklength = 8 #create the chunks of the reference sequence with regular expressions for the variable nucleotides -for i in range(0, len(searchstrings["ca"]) - chunklength, chunklength / 2): +for i in range(0, len(searchstrings["ca"]) - chunklength, chunklength // 2): pos = i chunk = searchstrings["ca"][i:i+chunklength] result = "" varsInResult = 0 for c in chunk: - if pos in ca1.keys(): + if pos in list(ca1.keys()): varsInResult += 1 result += "[" + ca1[pos] + ca2[pos] + "]" else: @@ -87,13 +87,13 @@ pos += 1 compiledregex["ca"].append((re.compile(result), varsInResult)) -for i in range(0, len(searchstrings["cg"]) - chunklength, chunklength / 2): +for i in range(0, len(searchstrings["cg"]) - chunklength, chunklength // 2): pos = i chunk = searchstrings["cg"][i:i+chunklength] result = "" varsInResult = 0 for c in chunk: - if pos in cg1.keys(): + if pos in list(cg1.keys()): varsInResult += 1 result += "[" + "".join(set([cg1[pos], cg2[pos], cg3[pos], cg4[pos]])) + "]" else: @@ -101,10 +101,10 @@ pos += 1 compiledregex["cg"].append((re.compile(result), varsInResult)) -for i in range(0, len(searchstrings["cm"]) - chunklength, chunklength / 2): +for i in range(0, len(searchstrings["cm"]) - chunklength, chunklength // 2): compiledregex["cm"].append((re.compile(searchstrings["cm"][i:i+chunklength]), False)) -for i in range(0, len(searchstrings["ce"]) - chunklength + 1, chunklength / 2): +for i in range(0, len(searchstrings["ce"]) - chunklength + 1, chunklength // 2): compiledregex["ce"].append((re.compile(searchstrings["ce"][i:i+chunklength]), False)) def removeAndReturnMaxIndex(x): #simplifies a list comprehension @@ -117,108 +117,108 @@ start_location = dict() hits = dict() alltotal = 0 -for key in compiledregex.keys(): #for ca/cg/cm/ce - regularexpressions = compiledregex[key] #get the compiled regular expressions - for ID in dic.keys()[0:]: #for every ID - if ID not in hits.keys(): #ensure that the dictionairy that keeps track of the hits for every gene exists - 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} - currentIDHits = hits[ID] - seq = dic[ID] - lastindex = 0 - start_zero = len(searchstrings[key]) #allows the reference sequence to start before search sequence (start_locations of < 0) - start = [0] * (len(seq) + start_zero) - for i, regexp in enumerate(regularexpressions): #for every regular expression - relativeStartLocation = lastindex - (chunklength / 2) * i - if relativeStartLocation >= len(seq): - break - regex, hasVar = regexp - matches = regex.finditer(seq[lastindex:]) - 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 - lastindex += match.start() - start[relativeStartLocation + start_zero] += 1 - if hasVar: #if the regex has a variable nt in it - chunkstart = chunklength / 2 * i #where in the reference does this chunk start - chunkend = chunklength / 2 * i + chunklength #where in the reference does this chunk end - if key == "ca": #just calculate the variable nt score for 'ca', cheaper - currentIDHits["ca1"] += len([1 for x in ca1 if chunkstart <= x < chunkend and ca1[x] == seq[lastindex + x - chunkstart]]) - currentIDHits["ca2"] += len([1 for x in ca2 if chunkstart <= x < chunkend and ca2[x] == seq[lastindex + x - chunkstart]]) - elif key == "cg": #just calculate the variable nt score for 'cg', cheaper - currentIDHits["cg1"] += len([1 for x in cg1 if chunkstart <= x < chunkend and cg1[x] == seq[lastindex + x - chunkstart]]) - currentIDHits["cg2"] += len([1 for x in cg2 if chunkstart <= x < chunkend and cg2[x] == seq[lastindex + x - chunkstart]]) - currentIDHits["cg3"] += len([1 for x in cg3 if chunkstart <= x < chunkend and cg3[x] == seq[lastindex + x - chunkstart]]) - currentIDHits["cg4"] += len([1 for x in cg4 if chunkstart <= x < chunkend and cg4[x] == seq[lastindex + x - chunkstart]]) - else: #key == "cm" #no variable regions in 'cm' or 'ce' - pass - break #this only breaks when there was a match with the regex, breaking means the 'else:' clause is skipped - else: #only runs if there were no hits - continue - #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] - currentIDHits[key + "_hits"] += 1 - start_location[ID + "_" + key] = str([(removeAndReturnMaxIndex(start) + 1 - start_zero) for x in range(5) if len(start) > 0 and max(start) > 1]) - #start_location[ID + "_" + key] = str(start.index(max(start))) +for key in compiledregex: #for ca/cg/cm/ce + regularexpressions = compiledregex[key] # get the compiled regular expressions + for ID in list(dic.keys())[0:]: #for every ID + if ID not in list(hits.keys()): #ensure that the dictionairy that keeps track of the hits for every gene exists + 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} + currentIDHits = hits[ID] + seq = dic[ID] + lastindex = 0 + start_zero = len(searchstrings[key]) #allows the reference sequence to start before search sequence (start_locations of < 0) + start = [0] * (len(seq) + start_zero) + for i, regexp in enumerate(regularexpressions): #for every regular expression + relativeStartLocation = lastindex - (chunklength // 2) * i + if relativeStartLocation >= len(seq): + break + regex, hasVar = regexp + matches = regex.finditer(seq[lastindex:]) + 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 + lastindex += match.start() + start[relativeStartLocation + start_zero] += 1 + if hasVar: #if the regex has a variable nt in it + chunkstart = chunklength // 2 * i #where in the reference does this chunk start + chunkend = chunklength // 2 * i + chunklength #where in the reference does this chunk end + if key == "ca": #just calculate the variable nt score for 'ca', cheaper + currentIDHits["ca1"] += len([1 for x in ca1 if chunkstart <= x < chunkend and ca1[x] == seq[lastindex + x - chunkstart]]) + currentIDHits["ca2"] += len([1 for x in ca2 if chunkstart <= x < chunkend and ca2[x] == seq[lastindex + x - chunkstart]]) + elif key == "cg": #just calculate the variable nt score for 'cg', cheaper + currentIDHits["cg1"] += len([1 for x in cg1 if chunkstart <= x < chunkend and cg1[x] == seq[lastindex + x - chunkstart]]) + currentIDHits["cg2"] += len([1 for x in cg2 if chunkstart <= x < chunkend and cg2[x] == seq[lastindex + x - chunkstart]]) + currentIDHits["cg3"] += len([1 for x in cg3 if chunkstart <= x < chunkend and cg3[x] == seq[lastindex + x - chunkstart]]) + currentIDHits["cg4"] += len([1 for x in cg4 if chunkstart <= x < chunkend and cg4[x] == seq[lastindex + x - chunkstart]]) + else: #key == "cm" #no variable regions in 'cm' or 'ce' + pass + break #this only breaks when there was a match with the regex, breaking means the 'else:' clause is skipped + else: #only runs if there were no hits + continue + #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] + currentIDHits[key + "_hits"] += 1 + start_location[ID + "_" + key] = str([(removeAndReturnMaxIndex(start) + 1 - start_zero) for x in range(5) if len(start) > 0 and max(start) > 1]) + #start_location[ID + "_" + key] = str(start.index(max(start))) -varsInCA = float(len(ca1.keys()) * 2) -varsInCG = float(len(cg1.keys()) * 2) - 2 # -2 because the sliding window doesn't hit the first and last nt twice +varsInCA = float(len(list(ca1.keys())) * 2) +varsInCG = float(len(list(cg1.keys())) * 2) - 2 # -2 because the sliding window doesn't hit the first and last nt twice varsInCM = 0 varsInCE = 0 def round_int(val): - return int(round(val)) + return int(round(val)) first = True seq_write_count=0 with open(infile, 'r') as f: #read all sequences into a dictionary as key = ID, value = sequence - with open(output, 'w') as o: - for line in f: - total += 1 - if first: - o.write("Sequence ID\tbest_match\tnt_hit_percentage\tchunk_hit_percentage\tstart_locations\n") - first = False - continue - linesplt = line.split("\t") - if linesplt[2] == "No results": - pass - ID = linesplt[1] - currentIDHits = hits[ID] - possibleca = float(len(compiledregex["ca"])) - possiblecg = float(len(compiledregex["cg"])) - possiblecm = float(len(compiledregex["cm"])) - possiblece = float(len(compiledregex["ce"])) - cahits = currentIDHits["ca_hits"] - cghits = currentIDHits["cg_hits"] - cmhits = currentIDHits["cm_hits"] - cehits = currentIDHits["ce_hits"] - if cahits >= cghits and cahits >= cmhits and cahits >= cehits: #its a ca gene - ca1hits = currentIDHits["ca1"] - ca2hits = currentIDHits["ca2"] - if ca1hits >= ca2hits: - o.write(ID + "\tIGA1\t" + str(round_int(ca1hits / varsInCA * 100)) + "\t" + str(round_int(cahits / possibleca * 100)) + "\t" + start_location[ID + "_ca"] + "\n") - else: - o.write(ID + "\tIGA2\t" + str(round_int(ca2hits / varsInCA * 100)) + "\t" + str(round_int(cahits / possibleca * 100)) + "\t" + start_location[ID + "_ca"] + "\n") - elif cghits >= cahits and cghits >= cmhits and cghits >= cehits: #its a cg gene - cg1hits = currentIDHits["cg1"] - cg2hits = currentIDHits["cg2"] - cg3hits = currentIDHits["cg3"] - cg4hits = currentIDHits["cg4"] - if cg1hits >= cg2hits and cg1hits >= cg3hits and cg1hits >= cg4hits: #cg1 gene - o.write(ID + "\tIGG1\t" + str(round_int(cg1hits / varsInCG * 100)) + "\t" + str(round_int(cghits / possiblecg * 100)) + "\t" + start_location[ID + "_cg"] + "\n") - elif cg2hits >= cg1hits and cg2hits >= cg3hits and cg2hits >= cg4hits: #cg2 gene - o.write(ID + "\tIGG2\t" + str(round_int(cg2hits / varsInCG * 100)) + "\t" + str(round_int(cghits / possiblecg * 100)) + "\t" + start_location[ID + "_cg"] + "\n") - elif cg3hits >= cg1hits and cg3hits >= cg2hits and cg3hits >= cg4hits: #cg3 gene - o.write(ID + "\tIGG3\t" + str(round_int(cg3hits / varsInCG * 100)) + "\t" + str(round_int(cghits / possiblecg * 100)) + "\t" + start_location[ID + "_cg"] + "\n") - else: #cg4 gene - o.write(ID + "\tIGG4\t" + str(round_int(cg4hits / varsInCG * 100)) + "\t" + str(round_int(cghits / possiblecg * 100)) + "\t" + start_location[ID + "_cg"] + "\n") - else: #its a cm or ce gene - if cmhits >= cehits: - o.write(ID + "\tIGM\t100\t" + str(round_int(cmhits / possiblecm * 100)) + "\t" + start_location[ID + "_cm"] + "\n") - else: - o.write(ID + "\tIGE\t100\t" + str(round_int(cehits / possiblece * 100)) + "\t" + start_location[ID + "_ce"] + "\n") - seq_write_count += 1 + with open(output, 'w') as o: + for line in f: + total += 1 + if first: + o.write("Sequence ID\tbest_match\tnt_hit_percentage\tchunk_hit_percentage\tstart_locations\n") + first = False + continue + linesplt = line.split("\t") + if linesplt[2] == "No results": + pass + ID = linesplt[1] + currentIDHits = hits[ID] + possibleca = float(len(compiledregex["ca"])) + possiblecg = float(len(compiledregex["cg"])) + possiblecm = float(len(compiledregex["cm"])) + possiblece = float(len(compiledregex["ce"])) + cahits = currentIDHits["ca_hits"] + cghits = currentIDHits["cg_hits"] + cmhits = currentIDHits["cm_hits"] + cehits = currentIDHits["ce_hits"] + if cahits >= cghits and cahits >= cmhits and cahits >= cehits: #its a ca gene + ca1hits = currentIDHits["ca1"] + ca2hits = currentIDHits["ca2"] + if ca1hits >= ca2hits: + o.write(ID + "\tIGA1\t" + str(round_int(ca1hits / varsInCA * 100)) + "\t" + str(round_int(cahits / possibleca * 100)) + "\t" + start_location[ID + "_ca"] + "\n") + else: + o.write(ID + "\tIGA2\t" + str(round_int(ca2hits / varsInCA * 100)) + "\t" + str(round_int(cahits / possibleca * 100)) + "\t" + start_location[ID + "_ca"] + "\n") + elif cghits >= cahits and cghits >= cmhits and cghits >= cehits: #its a cg gene + cg1hits = currentIDHits["cg1"] + cg2hits = currentIDHits["cg2"] + cg3hits = currentIDHits["cg3"] + cg4hits = currentIDHits["cg4"] + if cg1hits >= cg2hits and cg1hits >= cg3hits and cg1hits >= cg4hits: #cg1 gene + o.write(ID + "\tIGG1\t" + str(round_int(cg1hits / varsInCG * 100)) + "\t" + str(round_int(cghits / possiblecg * 100)) + "\t" + start_location[ID + "_cg"] + "\n") + elif cg2hits >= cg1hits and cg2hits >= cg3hits and cg2hits >= cg4hits: #cg2 gene + o.write(ID + "\tIGG2\t" + str(round_int(cg2hits / varsInCG * 100)) + "\t" + str(round_int(cghits / possiblecg * 100)) + "\t" + start_location[ID + "_cg"] + "\n") + elif cg3hits >= cg1hits and cg3hits >= cg2hits and cg3hits >= cg4hits: #cg3 gene + o.write(ID + "\tIGG3\t" + str(round_int(cg3hits / varsInCG * 100)) + "\t" + str(round_int(cghits / possiblecg * 100)) + "\t" + start_location[ID + "_cg"] + "\n") + else: #cg4 gene + o.write(ID + "\tIGG4\t" + str(round_int(cg4hits / varsInCG * 100)) + "\t" + str(round_int(cghits / possiblecg * 100)) + "\t" + start_location[ID + "_cg"] + "\n") + else: #its a cm or ce gene + if cmhits >= cehits: + o.write(ID + "\tIGM\t100\t" + str(round_int(cmhits / possiblecm * 100)) + "\t" + start_location[ID + "_cm"] + "\n") + else: + o.write(ID + "\tIGE\t100\t" + str(round_int(cehits / possiblece * 100)) + "\t" + start_location[ID + "_ce"] + "\n") + seq_write_count += 1 -print "Time: %i" % (int(time.time() * 1000) - starttime) +print("Time: %i" % (int(time.time() * 1000) - starttime)) -print "Number of sequences written to file:", seq_write_count +print("Number of sequences written to file:", seq_write_count)