annotate baseline/script_imgt.py @ 0:c33d93683a09 draft

Uploaded

1 #import xlrd #avoid dep
2 import argparse
3 import re
4 
5 parser = argparse.ArgumentParser()
6 parser.add_argument("--input", help="Excel input file containing one or more sheets where column G has the gene annotation, H has the sequence id and J has the sequence")
7 parser.add_argument("--ref", help="Reference file")
8 parser.add_argument("--output", help="Output file")
9 parser.add_argument("--id", help="ID to be used at the '>>>' line in the output")
10 
11 args = parser.parse_args()
12 
13 refdic = dict()
14 with open(args.ref, 'r') as ref:
15 	currentSeq = ""
16 	currentId = ""
17 	for line in ref:
18 		if line[0] is ">":
19 			if currentSeq is not "" and currentId is not "":
20 				refdic[currentId[1:]] = currentSeq
21 			currentId = line.rstrip()
22 			currentSeq = ""
23 		else:
24 			currentSeq += line.rstrip()
25 	refdic[currentId[1:]] = currentSeq
26 	
27 
28 vPattern = [r"(IGHV[0-9]-[0-9ab]+-?[0-9]?D?\*\d{1,2})"]#,
29 #						r"(TRBV[0-9]{1,2}-?[0-9]?-?[123]?)",
30 #						r"(IGKV[0-3]D?-[0-9]{1,2})",
31 #						r"(IGLV[0-9]-[0-9]{1,2})",
32 #						r"(TRAV[0-9]{1,2}(-[1-46])?(/DV[45678])?)",
33 #						r"(TRGV[234589])",
34 #						r"(TRDV[1-3])"]
35 
36 #vPattern = re.compile(r"|".join(vPattern))
37 vPattern = re.compile("|".join(vPattern))
38 
39 def filterGene(s, pattern):
40     if type(s) is not str:
41         return None
42     res = pattern.search(s)
43     if res:
44         return res.group(0)
45     return None
46 
47 
48 
49 currentSeq = ""
50 currentId = ""
51 first=True
52 with open(args.input, 'r') as i:
53 	with open(args.output, 'a') as o:
54 		o.write(">>>" + args.id + "\n")
55 		outputdic = dict()
56 		for line in i:
57 			if first:
58 				first = False
59 				continue
60 			linesplt = line.split("\t")
61 			ref = filterGene(linesplt[1], vPattern)
62 			if not ref or not linesplt[2].rstrip():
63 				continue
64 			if ref in outputdic:
65 				outputdic[ref] += [(linesplt[0].replace(">", ""), linesplt[2].replace(">", "").rstrip())]
66 			else:
67 				outputdic[ref] = [(linesplt[0].replace(">", ""), linesplt[2].replace(">", "").rstrip())]
68 		#print outputdic
69 		
70 		for k in outputdic.keys():
71 			if k in refdic:
72 				o.write(">>" + k + "\n")
73 				o.write(refdic[k] + "\n")
74 				for seq in outputdic[k]:
75 					#print seq
76 					o.write(">" + seq[0] + "\n")
77 					o.write(seq[1] + "\n")
78 			else:
79 				print k + " not in reference, skipping " + k