Mercurial > repos > davidvanzessen > shm_csr

comparison baseline/script_imgt.py @ 0:c33d93683a09 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author davidvanzessen
date Thu, 13 Oct 2016 10:52:24 -0400
parents
children 4c5ba6b5d10d
comparison
equal deleted inserted replaced
-1:000000000000 0:c33d93683a09
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