Mercurial > repos > rhpvorderman > shm_csr

annotate shm_csr.py @ 6:ea9d5fc4c001 draft default tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
"planemo upload commit 9ada186a78831ca2618ec817a23a77de6adf1a5d"
author rhpvorderman
date Wed, 22 Dec 2021 11:29:16 +0000
parents 64d74ba01a7c
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
0
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
1 import argparse
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
2 import logging
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
3 import sys
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
4 import os
6
ea9d5fc4c001 "planemo upload commit 9ada186a78831ca2618ec817a23a77de6adf1a5d"
rhpvorderman
parents: 0
diff changeset
5 import typing
ea9d5fc4c001 "planemo upload commit 9ada186a78831ca2618ec817a23a77de6adf1a5d"
rhpvorderman
parents: 0
diff changeset
6 from typing import Optional
0
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
7
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
8 from collections import defaultdict
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
9
6
ea9d5fc4c001 "planemo upload commit 9ada186a78831ca2618ec817a23a77de6adf1a5d"
rhpvorderman
parents: 0
diff changeset
10 REGION_FILTERS = ("leader", "FR1", "CDR1", "FR2", "CDR2")
ea9d5fc4c001 "planemo upload commit 9ada186a78831ca2618ec817a23a77de6adf1a5d"
rhpvorderman
parents: 0
diff changeset
11
ea9d5fc4c001 "planemo upload commit 9ada186a78831ca2618ec817a23a77de6adf1a5d"
rhpvorderman
parents: 0
diff changeset
12
ea9d5fc4c001 "planemo upload commit 9ada186a78831ca2618ec817a23a77de6adf1a5d"
rhpvorderman
parents: 0
diff changeset
13 class Mutation(typing.NamedTuple):
ea9d5fc4c001 "planemo upload commit 9ada186a78831ca2618ec817a23a77de6adf1a5d"
rhpvorderman
parents: 0
diff changeset
14 """Represent a mutation type as a tuple"""
ea9d5fc4c001 "planemo upload commit 9ada186a78831ca2618ec817a23a77de6adf1a5d"
rhpvorderman
parents: 0
diff changeset
15 frm: str # 'from' is a reserved python keyword.
ea9d5fc4c001 "planemo upload commit 9ada186a78831ca2618ec817a23a77de6adf1a5d"
rhpvorderman
parents: 0
diff changeset
16 where: int
ea9d5fc4c001 "planemo upload commit 9ada186a78831ca2618ec817a23a77de6adf1a5d"
rhpvorderman
parents: 0
diff changeset
17 to: str
ea9d5fc4c001 "planemo upload commit 9ada186a78831ca2618ec817a23a77de6adf1a5d"
rhpvorderman
parents: 0
diff changeset
18 frmAA: Optional[str] = None
ea9d5fc4c001 "planemo upload commit 9ada186a78831ca2618ec817a23a77de6adf1a5d"
rhpvorderman
parents: 0
diff changeset
19 whereAA: Optional[int] = None
ea9d5fc4c001 "planemo upload commit 9ada186a78831ca2618ec817a23a77de6adf1a5d"
rhpvorderman
parents: 0
diff changeset
20 toAA: Optional[str] = None
ea9d5fc4c001 "planemo upload commit 9ada186a78831ca2618ec817a23a77de6adf1a5d"
rhpvorderman
parents: 0
diff changeset
21 thing: Optional[str] = None # '(---)' or '(+-+)' etc. No idea
ea9d5fc4c001 "planemo upload commit 9ada186a78831ca2618ec817a23a77de6adf1a5d"
rhpvorderman
parents: 0
diff changeset
22
ea9d5fc4c001 "planemo upload commit 9ada186a78831ca2618ec817a23a77de6adf1a5d"
rhpvorderman
parents: 0
diff changeset
23 @classmethod
ea9d5fc4c001 "planemo upload commit 9ada186a78831ca2618ec817a23a77de6adf1a5d"
rhpvorderman
parents: 0
diff changeset
24 def from_string(cls, string: str):
ea9d5fc4c001 "planemo upload commit 9ada186a78831ca2618ec817a23a77de6adf1a5d"
rhpvorderman
parents: 0
diff changeset
25 # Complete mutation example: a88>g,I30>V(+ - +)
ea9d5fc4c001 "planemo upload commit 9ada186a78831ca2618ec817a23a77de6adf1a5d"
rhpvorderman
parents: 0
diff changeset
26 # Only nucleotide example: g303>t
ea9d5fc4c001 "planemo upload commit 9ada186a78831ca2618ec817a23a77de6adf1a5d"
rhpvorderman
parents: 0
diff changeset
27 if ',' in string:
ea9d5fc4c001 "planemo upload commit 9ada186a78831ca2618ec817a23a77de6adf1a5d"
rhpvorderman
parents: 0
diff changeset
28 nucleotide_change, aa_change = string.split(',', maxsplit=1) # type: str, Optional[str]
ea9d5fc4c001 "planemo upload commit 9ada186a78831ca2618ec817a23a77de6adf1a5d"
rhpvorderman
parents: 0
diff changeset
29 else:
ea9d5fc4c001 "planemo upload commit 9ada186a78831ca2618ec817a23a77de6adf1a5d"
rhpvorderman
parents: 0
diff changeset
30 nucleotide_change = string
ea9d5fc4c001 "planemo upload commit 9ada186a78831ca2618ec817a23a77de6adf1a5d"
rhpvorderman
parents: 0
diff changeset
31 aa_change = None
ea9d5fc4c001 "planemo upload commit 9ada186a78831ca2618ec817a23a77de6adf1a5d"
rhpvorderman
parents: 0
diff changeset
32 frm_part, to = nucleotide_change.split('>', maxsplit=1)
ea9d5fc4c001 "planemo upload commit 9ada186a78831ca2618ec817a23a77de6adf1a5d"
rhpvorderman
parents: 0
diff changeset
33 frm = frm_part[0]
ea9d5fc4c001 "planemo upload commit 9ada186a78831ca2618ec817a23a77de6adf1a5d"
rhpvorderman
parents: 0
diff changeset
34 where = int(frm_part[1:])
ea9d5fc4c001 "planemo upload commit 9ada186a78831ca2618ec817a23a77de6adf1a5d"
rhpvorderman
parents: 0
diff changeset
35
ea9d5fc4c001 "planemo upload commit 9ada186a78831ca2618ec817a23a77de6adf1a5d"
rhpvorderman
parents: 0
diff changeset
36 if aa_change is None:
ea9d5fc4c001 "planemo upload commit 9ada186a78831ca2618ec817a23a77de6adf1a5d"
rhpvorderman
parents: 0
diff changeset
37 return cls(frm, where, to)
ea9d5fc4c001 "planemo upload commit 9ada186a78831ca2618ec817a23a77de6adf1a5d"
rhpvorderman
parents: 0
diff changeset
38
ea9d5fc4c001 "planemo upload commit 9ada186a78831ca2618ec817a23a77de6adf1a5d"
rhpvorderman
parents: 0
diff changeset
39 frmAA_part, toAA_part = aa_change.split('>', maxsplit=1) # type: str, str
ea9d5fc4c001 "planemo upload commit 9ada186a78831ca2618ec817a23a77de6adf1a5d"
rhpvorderman
parents: 0
diff changeset
40 frmAA = frmAA_part[0]
ea9d5fc4c001 "planemo upload commit 9ada186a78831ca2618ec817a23a77de6adf1a5d"
rhpvorderman
parents: 0
diff changeset
41 whereAA = int(frmAA_part[1:])
ea9d5fc4c001 "planemo upload commit 9ada186a78831ca2618ec817a23a77de6adf1a5d"
rhpvorderman
parents: 0
diff changeset
42 brace_start = toAA_part.index('(')
ea9d5fc4c001 "planemo upload commit 9ada186a78831ca2618ec817a23a77de6adf1a5d"
rhpvorderman
parents: 0
diff changeset
43 toAA = toAA_part[:brace_start]
ea9d5fc4c001 "planemo upload commit 9ada186a78831ca2618ec817a23a77de6adf1a5d"
rhpvorderman
parents: 0
diff changeset
44 thing = toAA_part[brace_start:]
ea9d5fc4c001 "planemo upload commit 9ada186a78831ca2618ec817a23a77de6adf1a5d"
rhpvorderman
parents: 0
diff changeset
45 return cls(frm, where, to, frmAA, whereAA, toAA, thing)
ea9d5fc4c001 "planemo upload commit 9ada186a78831ca2618ec817a23a77de6adf1a5d"
rhpvorderman
parents: 0
diff changeset
46
ea9d5fc4c001 "planemo upload commit 9ada186a78831ca2618ec817a23a77de6adf1a5d"
rhpvorderman
parents: 0
diff changeset
47
ea9d5fc4c001 "planemo upload commit 9ada186a78831ca2618ec817a23a77de6adf1a5d"
rhpvorderman
parents: 0
diff changeset
48 class Hotspot(typing.NamedTuple):
ea9d5fc4c001 "planemo upload commit 9ada186a78831ca2618ec817a23a77de6adf1a5d"
rhpvorderman
parents: 0
diff changeset
49 start: int
ea9d5fc4c001 "planemo upload commit 9ada186a78831ca2618ec817a23a77de6adf1a5d"
rhpvorderman
parents: 0
diff changeset
50 end: int
ea9d5fc4c001 "planemo upload commit 9ada186a78831ca2618ec817a23a77de6adf1a5d"
rhpvorderman
parents: 0
diff changeset
51 region: str
ea9d5fc4c001 "planemo upload commit 9ada186a78831ca2618ec817a23a77de6adf1a5d"
rhpvorderman
parents: 0
diff changeset
52
ea9d5fc4c001 "planemo upload commit 9ada186a78831ca2618ec817a23a77de6adf1a5d"
rhpvorderman
parents: 0
diff changeset
53 @classmethod
ea9d5fc4c001 "planemo upload commit 9ada186a78831ca2618ec817a23a77de6adf1a5d"
rhpvorderman
parents: 0
diff changeset
54 def from_string(cls, string):
ea9d5fc4c001 "planemo upload commit 9ada186a78831ca2618ec817a23a77de6adf1a5d"
rhpvorderman
parents: 0
diff changeset
55 # Example: aa,40-41(FR1)
ea9d5fc4c001 "planemo upload commit 9ada186a78831ca2618ec817a23a77de6adf1a5d"
rhpvorderman
parents: 0
diff changeset
56 sequence, rest = string.split(',') # type: str, str
ea9d5fc4c001 "planemo upload commit 9ada186a78831ca2618ec817a23a77de6adf1a5d"
rhpvorderman
parents: 0
diff changeset
57 brace_pos = rest.index('(')
ea9d5fc4c001 "planemo upload commit 9ada186a78831ca2618ec817a23a77de6adf1a5d"
rhpvorderman
parents: 0
diff changeset
58 numbers = rest[:brace_pos]
ea9d5fc4c001 "planemo upload commit 9ada186a78831ca2618ec817a23a77de6adf1a5d"
rhpvorderman
parents: 0
diff changeset
59 start, end = numbers.split('-')
ea9d5fc4c001 "planemo upload commit 9ada186a78831ca2618ec817a23a77de6adf1a5d"
rhpvorderman
parents: 0
diff changeset
60 region = rest[brace_pos + 1:-1] # Remove the braces
ea9d5fc4c001 "planemo upload commit 9ada186a78831ca2618ec817a23a77de6adf1a5d"
rhpvorderman
parents: 0
diff changeset
61 return cls(int(start), int(end), region)
ea9d5fc4c001 "planemo upload commit 9ada186a78831ca2618ec817a23a77de6adf1a5d"
rhpvorderman
parents: 0
diff changeset
62
ea9d5fc4c001 "planemo upload commit 9ada186a78831ca2618ec817a23a77de6adf1a5d"
rhpvorderman
parents: 0
diff changeset
63
0
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
64 def main():
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
65 parser = argparse.ArgumentParser()
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
66 parser.add_argument("--input", help="The '7_V-REGION-mutation-and-AA-change-table' and '10_V-REGION-mutation-hotspots' merged together, with an added 'best_match' annotation")
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
67 parser.add_argument("--genes", help="The genes available in the 'best_match' column")
6
ea9d5fc4c001 "planemo upload commit 9ada186a78831ca2618ec817a23a77de6adf1a5d"
rhpvorderman
parents: 0
diff changeset
68 parser.add_argument("--empty_region_filter", help="Where does the sequence start?", choices=REGION_FILTERS)
0
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
69 parser.add_argument("--output", help="Output file")
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
70
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
71 args = parser.parse_args()
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
72
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
73 infile = args.input
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
74 genes = str(args.genes).split(",")
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
75 empty_region_filter = args.empty_region_filter
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
76 outfile = args.output
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
77
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
78 genedic = dict()
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
79
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
80 mutationdic = dict()
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
81 NAMatchResult = (None, None, None, None, None, None, '')
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
82 linecount = 0
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
83
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
84 IDIndex = 0
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
85 best_matchIndex = 0
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
86 fr1Index = 0
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
87 cdr1Index = 0
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
88 fr2Index = 0
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
89 cdr2Index = 0
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
90 fr3Index = 0
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
91 first = True
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
92 IDlist = []
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
93 mutationList = []
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
94 mutationListByID = {}
6
ea9d5fc4c001 "planemo upload commit 9ada186a78831ca2618ec817a23a77de6adf1a5d"
rhpvorderman
parents: 0
diff changeset
95 cdr1AALengthDic = {}
ea9d5fc4c001 "planemo upload commit 9ada186a78831ca2618ec817a23a77de6adf1a5d"
rhpvorderman
parents: 0
diff changeset
96 cdr2AALengthDic = {}
0
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
97
6
ea9d5fc4c001 "planemo upload commit 9ada186a78831ca2618ec817a23a77de6adf1a5d"
rhpvorderman
parents: 0
diff changeset
98 LengthDic = {}
0
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
99
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
100 cdr1LengthIndex = 0
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
101 cdr2LengthIndex = 0
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
102
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
103 tandem_sum_by_class = defaultdict(int)
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
104 expected_tandem_sum_by_class = defaultdict(float)
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
105
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
106 with open(infile, 'r') as i:
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
107 for line in i:
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
108 if first:
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
109 linesplt = line.split("\t")
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
110 IDIndex = linesplt.index("Sequence.ID")
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
111 best_matchIndex = linesplt.index("best_match")
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
112 fr1Index = linesplt.index("FR1.IMGT")
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
113 cdr1Index = linesplt.index("CDR1.IMGT")
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
114 fr2Index = linesplt.index("FR2.IMGT")
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
115 cdr2Index = linesplt.index("CDR2.IMGT")
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
116 fr3Index = linesplt.index("FR3.IMGT")
6
ea9d5fc4c001 "planemo upload commit 9ada186a78831ca2618ec817a23a77de6adf1a5d"
rhpvorderman
parents: 0
diff changeset
117 fr1LengthIndex = linesplt.index("FR1.IMGT.Nb.of.nucleotides")
ea9d5fc4c001 "planemo upload commit 9ada186a78831ca2618ec817a23a77de6adf1a5d"
rhpvorderman
parents: 0
diff changeset
118 fr2LengthIndex = linesplt.index("FR2.IMGT.Nb.of.nucleotides")
ea9d5fc4c001 "planemo upload commit 9ada186a78831ca2618ec817a23a77de6adf1a5d"
rhpvorderman
parents: 0
diff changeset
119 fr3LengthIndex = linesplt.index("FR3.IMGT.Nb.of.nucleotides")
ea9d5fc4c001 "planemo upload commit 9ada186a78831ca2618ec817a23a77de6adf1a5d"
rhpvorderman
parents: 0
diff changeset
120 cdr1LengthIndex = linesplt.index("CDR1.IMGT.Nb.of.nucleotides")
ea9d5fc4c001 "planemo upload commit 9ada186a78831ca2618ec817a23a77de6adf1a5d"
rhpvorderman
parents: 0
diff changeset
121 cdr2LengthIndex = linesplt.index("CDR2.IMGT.Nb.of.nucleotides")
ea9d5fc4c001 "planemo upload commit 9ada186a78831ca2618ec817a23a77de6adf1a5d"
rhpvorderman
parents: 0
diff changeset
122 cdr1AALengthIndex = linesplt.index("CDR1.IMGT.length")
ea9d5fc4c001 "planemo upload commit 9ada186a78831ca2618ec817a23a77de6adf1a5d"
rhpvorderman
parents: 0
diff changeset
123 cdr2AALengthIndex = linesplt.index("CDR2.IMGT.length")
0
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
124 first = False
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
125 continue
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
126 linecount += 1
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
127 linesplt = line.split("\t")
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
128 ID = linesplt[IDIndex]
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
129 genedic[ID] = linesplt[best_matchIndex]
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
130
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
131 mutationdic[ID + "_FR1"] = []
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
132 if len(linesplt[fr1Index]) > 5 and empty_region_filter == "leader":
6
ea9d5fc4c001 "planemo upload commit 9ada186a78831ca2618ec817a23a77de6adf1a5d"
rhpvorderman
parents: 0
diff changeset
133 mutationdic[ID + "_FR1"] = [Mutation.from_string(x) for x in linesplt[fr1Index].split("|") if x]
0
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
134
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
135 mutationdic[ID + "_CDR1"] = []
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
136 if len(linesplt[cdr1Index]) > 5 and empty_region_filter in ["leader", "FR1"]:
6
ea9d5fc4c001 "planemo upload commit 9ada186a78831ca2618ec817a23a77de6adf1a5d"
rhpvorderman
parents: 0
diff changeset
137 mutationdic[ID + "_CDR1"] = [Mutation.from_string(x) for x in linesplt[cdr1Index].split("|") if x]
0
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
138
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
139 mutationdic[ID + "_FR2"] = []
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
140 if len(linesplt[fr2Index]) > 5 and empty_region_filter in ["leader", "FR1", "CDR1"]:
6
ea9d5fc4c001 "planemo upload commit 9ada186a78831ca2618ec817a23a77de6adf1a5d"
rhpvorderman
parents: 0
diff changeset
141 mutationdic[ID + "_FR2"] = [Mutation.from_string(x) for x in linesplt[fr2Index].split("|") if x]
0
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
142
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
143 mutationdic[ID + "_CDR2"] = []
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
144 if len(linesplt[cdr2Index]) > 5:
6
ea9d5fc4c001 "planemo upload commit 9ada186a78831ca2618ec817a23a77de6adf1a5d"
rhpvorderman
parents: 0
diff changeset
145 mutationdic[ID + "_CDR2"] = [Mutation.from_string(x) for x in linesplt[cdr2Index].split("|") if x]
0
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
146
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
147 mutationdic[ID + "_FR2-CDR2"] = mutationdic[ID + "_FR2"] + mutationdic[ID + "_CDR2"]
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
148
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
149 mutationdic[ID + "_FR3"] = []
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
150 if len(linesplt[fr3Index]) > 5:
6
ea9d5fc4c001 "planemo upload commit 9ada186a78831ca2618ec817a23a77de6adf1a5d"
rhpvorderman
parents: 0
diff changeset
151 mutationdic[ID + "_FR3"] = [Mutation.from_string(x) for x in linesplt[fr3Index].split("|") if x]
0
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
152
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
153 mutationList += mutationdic[ID + "_FR1"] + mutationdic[ID + "_CDR1"] + mutationdic[ID + "_FR2"] + mutationdic[ID + "_CDR2"] + mutationdic[ID + "_FR3"]
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
154 mutationListByID[ID] = mutationdic[ID + "_FR1"] + mutationdic[ID + "_CDR1"] + mutationdic[ID + "_FR2"] + mutationdic[ID + "_CDR2"] + mutationdic[ID + "_FR3"]
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
155
6
ea9d5fc4c001 "planemo upload commit 9ada186a78831ca2618ec817a23a77de6adf1a5d"
rhpvorderman
parents: 0
diff changeset
156 fr1Length = int(linesplt[fr1LengthIndex])
ea9d5fc4c001 "planemo upload commit 9ada186a78831ca2618ec817a23a77de6adf1a5d"
rhpvorderman
parents: 0
diff changeset
157 fr2Length = int(linesplt[fr2LengthIndex])
ea9d5fc4c001 "planemo upload commit 9ada186a78831ca2618ec817a23a77de6adf1a5d"
rhpvorderman
parents: 0
diff changeset
158 fr3Length = int(linesplt[fr3LengthIndex])
ea9d5fc4c001 "planemo upload commit 9ada186a78831ca2618ec817a23a77de6adf1a5d"
rhpvorderman
parents: 0
diff changeset
159 cdr1Length = int(linesplt[cdr1LengthIndex])
ea9d5fc4c001 "planemo upload commit 9ada186a78831ca2618ec817a23a77de6adf1a5d"
rhpvorderman
parents: 0
diff changeset
160 cdr2Length = int(linesplt[cdr2LengthIndex])
ea9d5fc4c001 "planemo upload commit 9ada186a78831ca2618ec817a23a77de6adf1a5d"
rhpvorderman
parents: 0
diff changeset
161 LengthDic[ID] = (fr1Length, cdr1Length, fr2Length, cdr2Length, fr3Length)
0
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
162
6
ea9d5fc4c001 "planemo upload commit 9ada186a78831ca2618ec817a23a77de6adf1a5d"
rhpvorderman
parents: 0
diff changeset
163 cdr1AALengthDic[ID] = int(linesplt[cdr1AALengthIndex])
ea9d5fc4c001 "planemo upload commit 9ada186a78831ca2618ec817a23a77de6adf1a5d"
rhpvorderman
parents: 0
diff changeset
164 cdr2AALengthDic[ID] = int(linesplt[cdr2AALengthIndex])
0
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
165
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
166 IDlist += [ID]
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
167 print("len(mutationdic) =", len(mutationdic))
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
168
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
169 with open(os.path.join(os.path.dirname(os.path.abspath(infile)), "mutationdict.txt"), 'w') as out_handle:
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
170 for ID, lst in mutationdic.items():
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
171 for mut in lst:
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
172 out_handle.write("{0}\t{1}\n".format(ID, "\t".join([str(x) for x in mut])))
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
173
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
174 #tandem mutation stuff
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
175 tandem_frequency = defaultdict(int)
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
176 mutation_frequency = defaultdict(int)
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
177
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
178 mutations_by_id_dic = {}
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
179 first = True
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
180 mutation_by_id_file = os.path.join(os.path.dirname(outfile), "mutation_by_id.txt")
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
181 with open(mutation_by_id_file, 'r') as mutation_by_id:
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
182 for l in mutation_by_id:
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
183 if first:
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
184 first = False
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
185 continue
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
186 splt = l.split("\t")
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
187 mutations_by_id_dic[splt[0]] = int(splt[1])
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
188
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
189 tandem_file = os.path.join(os.path.dirname(outfile), "tandems_by_id.txt")
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
190 with open(tandem_file, 'w') as o:
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
191 highest_tandem_length = 0
6
ea9d5fc4c001 "planemo upload commit 9ada186a78831ca2618ec817a23a77de6adf1a5d"
rhpvorderman
parents: 0
diff changeset
192 # LengthDic stores length as a tuple
ea9d5fc4c001 "planemo upload commit 9ada186a78831ca2618ec817a23a77de6adf1a5d"
rhpvorderman
parents: 0
diff changeset
193 # (fr1Length, cdr1Length, fr2Length, cdr2Length, fr3Length)
ea9d5fc4c001 "planemo upload commit 9ada186a78831ca2618ec817a23a77de6adf1a5d"
rhpvorderman
parents: 0
diff changeset
194 # To get the total length, we can sum(region_lengths)
ea9d5fc4c001 "planemo upload commit 9ada186a78831ca2618ec817a23a77de6adf1a5d"
rhpvorderman
parents: 0
diff changeset
195 # To get the total length for leader:
ea9d5fc4c001 "planemo upload commit 9ada186a78831ca2618ec817a23a77de6adf1a5d"
rhpvorderman
parents: 0
diff changeset
196 # sum(region_lengths[0:]) (Equivalent to everything)
ea9d5fc4c001 "planemo upload commit 9ada186a78831ca2618ec817a23a77de6adf1a5d"
rhpvorderman
parents: 0
diff changeset
197 # sum(region_lengths[1:]) Gets everything except FR1 etc.
ea9d5fc4c001 "planemo upload commit 9ada186a78831ca2618ec817a23a77de6adf1a5d"
rhpvorderman
parents: 0
diff changeset
198 # We determine the position to start summing below.
ea9d5fc4c001 "planemo upload commit 9ada186a78831ca2618ec817a23a77de6adf1a5d"
rhpvorderman
parents: 0
diff changeset
199 # This returns 0 for leader, 1 for FR1 etc.
ea9d5fc4c001 "planemo upload commit 9ada186a78831ca2618ec817a23a77de6adf1a5d"
rhpvorderman
parents: 0
diff changeset
200 length_start_pos = REGION_FILTERS.index(empty_region_filter)
0
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
201
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
202 o.write("Sequence.ID\tnumber_of_mutations\tnumber_of_tandems\tregion_length\texpected_tandems\tlongest_tandem\ttandems\n")
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
203 for ID in IDlist:
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
204 mutations = mutationListByID[ID]
6
ea9d5fc4c001 "planemo upload commit 9ada186a78831ca2618ec817a23a77de6adf1a5d"
rhpvorderman
parents: 0
diff changeset
205 region_length = sum(LengthDic[ID][length_start_pos:])
0
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
206 if len(mutations) == 0:
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
207 continue
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
208 last_mut = max(mutations, key=lambda x: int(x[1]))
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
209
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
210 last_mut_pos = int(last_mut[1])
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
211
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
212 mut_positions = [False] * (last_mut_pos + 1)
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
213
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
214 for mutation in mutations:
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
215 frm, where, to, frmAA, whereAA, toAA, thing = mutation
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
216 where = int(where)
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
217 mut_positions[where] = True
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
218
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
219 tandem_muts = []
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
220 tandem_start = -1
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
221 tandem_length = 0
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
222 for i in range(len(mut_positions)):
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
223 if mut_positions[i]:
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
224 if tandem_start == -1:
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
225 tandem_start = i
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
226 tandem_length += 1
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
227 #print "".join(["1" if x else "0" for x in mut_positions[:i+1]])
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
228 else:
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
229 if tandem_length > 1:
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
230 tandem_muts.append((tandem_start, tandem_length))
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
231 #print "{0}{1} {2}:{3}".format(" " * (i - tandem_length), "^" * tandem_length, tandem_start, tandem_length)
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
232 tandem_start = -1
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
233 tandem_length = 0
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
234 if tandem_length > 1: # if the sequence ends with a tandem mutation
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
235 tandem_muts.append((tandem_start, tandem_length))
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
236
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
237 if len(tandem_muts) > 0:
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
238 if highest_tandem_length < len(tandem_muts):
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
239 highest_tandem_length = len(tandem_muts)
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
240
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
241 longest_tandem = max(tandem_muts, key=lambda x: x[1]) if len(tandem_muts) else (0, 0)
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
242 num_mutations = mutations_by_id_dic[ID] # len(mutations)
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
243 f_num_mutations = float(num_mutations)
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
244 num_tandem_muts = len(tandem_muts)
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
245 expected_tandem_muts = f_num_mutations * (f_num_mutations - 1.0) / float(region_length)
6
ea9d5fc4c001 "planemo upload commit 9ada186a78831ca2618ec817a23a77de6adf1a5d"
rhpvorderman
parents: 0
diff changeset
246 # String format and round disagree slightly (see 3.605).
ea9d5fc4c001 "planemo upload commit 9ada186a78831ca2618ec817a23a77de6adf1a5d"
rhpvorderman
parents: 0
diff changeset
247 # So round before formatting.
ea9d5fc4c001 "planemo upload commit 9ada186a78831ca2618ec817a23a77de6adf1a5d"
rhpvorderman
parents: 0
diff changeset
248 o.write(f"{ID}\t{num_mutations}\t{num_tandem_muts}\t{region_length}\t"
ea9d5fc4c001 "planemo upload commit 9ada186a78831ca2618ec817a23a77de6adf1a5d"
rhpvorderman
parents: 0
diff changeset
249 f"{round(expected_tandem_muts, 2):.2f}\t"
ea9d5fc4c001 "planemo upload commit 9ada186a78831ca2618ec817a23a77de6adf1a5d"
rhpvorderman
parents: 0
diff changeset
250 f"{longest_tandem[1]}\t{tandem_muts}\n")
0
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
251 gene = genedic[ID]
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
252 if gene.find("unmatched") == -1:
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
253 tandem_sum_by_class[gene] += num_tandem_muts
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
254 expected_tandem_sum_by_class[gene] += expected_tandem_muts
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
255
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
256 tandem_sum_by_class["all"] += num_tandem_muts
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
257 expected_tandem_sum_by_class["all"] += expected_tandem_muts
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
258
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
259 gene = gene[:3]
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
260 if gene in ["IGA", "IGG"]:
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
261 tandem_sum_by_class[gene] += num_tandem_muts
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
262 expected_tandem_sum_by_class[gene] += expected_tandem_muts
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
263 else:
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
264 tandem_sum_by_class["unmatched"] += num_tandem_muts
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
265 expected_tandem_sum_by_class["unmatched"] += expected_tandem_muts
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
266
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
267
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
268 for tandem_mut in tandem_muts:
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
269 tandem_frequency[str(tandem_mut[1])] += 1
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
270 #print "\t".join([ID, str(len(tandem_muts)), str(longest_tandem[1]) , str(tandem_muts)])
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
271
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
272 tandem_freq_file = os.path.join(os.path.dirname(outfile), "tandem_frequency.txt")
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
273 with open(tandem_freq_file, 'w') as o:
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
274 for frq in sorted([int(x) for x in list(tandem_frequency.keys())]):
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
275 o.write("{0}\t{1}\n".format(frq, tandem_frequency[str(frq)]))
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
276
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
277 tandem_row = []
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
278 genes_extra = list(genes)
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
279 genes_extra.append("all")
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
280 for x, y, in zip([tandem_sum_by_class[x] for x in genes_extra], [expected_tandem_sum_by_class[x] for x in genes_extra]):
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
281 if y != 0:
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
282 tandem_row += [x, round(y, 2), round(x / y, 2)]
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
283 else:
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
284 tandem_row += [x, round(y, 2), 0]
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
285
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
286 tandem_freq_file = os.path.join(os.path.dirname(outfile), "shm_overview_tandem_row.txt")
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
287 with open(tandem_freq_file, 'w') as o:
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
288 o.write("Tandems/Expected (ratio),{0}\n".format(",".join([str(x) for x in tandem_row])))
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
289
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
290 #print mutationList, linecount
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
291
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
292 AALength = (int(max(mutationList, key=lambda i: int(i[4]) if i[4] and i[5] != ";" else 0)[4]) + 1) # [4] is the position of the AA mutation, None if silent
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
293 if AALength < 60:
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
294 AALength = 64
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
295
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
296 AA_mutation = [0] * AALength
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
297 AA_mutation_dic = {"IGA": AA_mutation[:], "IGG": AA_mutation[:], "IGM": AA_mutation[:], "IGE": AA_mutation[:], "unm": AA_mutation[:], "all": AA_mutation[:]}
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
298 AA_mutation_empty = AA_mutation[:]
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
299
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
300 print("AALength:", AALength)
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
301 aa_mutations_by_id_file = outfile[:outfile.rindex("/")] + "/aa_id_mutations.txt"
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
302 with open(aa_mutations_by_id_file, 'w') as o:
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
303 o.write("ID\tbest_match\t" + "\t".join([str(x) for x in range(1,AALength)]) + "\n")
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
304 for ID in list(mutationListByID.keys()):
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
305 AA_mutation_for_ID = AA_mutation_empty[:]
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
306 for mutation in mutationListByID[ID]:
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
307 if mutation[4] and mutation[5] != ";":
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
308 AA_mutation_position = int(mutation[4])
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
309 try:
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
310 AA_mutation[AA_mutation_position] += 1
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
311 AA_mutation_for_ID[AA_mutation_position] += 1
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
312 except Exception as e:
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
313 print(e)
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
314 print(mutation)
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
315 sys.exit()
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
316 clss = genedic[ID][:3]
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
317 AA_mutation_dic[clss][AA_mutation_position] += 1
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
318 o.write(ID + "\t" + genedic[ID] + "\t" + "\t".join([str(x) for x in AA_mutation_for_ID[1:]]) + "\n")
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
319
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
320
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
321
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
322 #absent AA stuff
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
323 absentAACDR1Dic = defaultdict(list)
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
324 absentAACDR1Dic[5] = list(range(29,36))
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
325 absentAACDR1Dic[6] = list(range(29,35))
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
326 absentAACDR1Dic[7] = list(range(30,35))
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
327 absentAACDR1Dic[8] = list(range(30,34))
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
328 absentAACDR1Dic[9] = list(range(31,34))
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
329 absentAACDR1Dic[10] = list(range(31,33))
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
330 absentAACDR1Dic[11] = [32]
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
331
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
332 absentAACDR2Dic = defaultdict(list)
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
333 absentAACDR2Dic[0] = list(range(55,65))
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
334 absentAACDR2Dic[1] = list(range(56,65))
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
335 absentAACDR2Dic[2] = list(range(56,64))
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
336 absentAACDR2Dic[3] = list(range(57,64))
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
337 absentAACDR2Dic[4] = list(range(57,63))
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
338 absentAACDR2Dic[5] = list(range(58,63))
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
339 absentAACDR2Dic[6] = list(range(58,62))
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
340 absentAACDR2Dic[7] = list(range(59,62))
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
341 absentAACDR2Dic[8] = list(range(59,61))
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
342 absentAACDR2Dic[9] = [60]
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
343
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
344 absentAA = [len(IDlist)] * (AALength-1)
6
ea9d5fc4c001 "planemo upload commit 9ada186a78831ca2618ec817a23a77de6adf1a5d"
rhpvorderman
parents: 0
diff changeset
345 for k, cdr1Length in cdr1AALengthDic.items():
0
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
346 for c in absentAACDR1Dic[cdr1Length]:
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
347 absentAA[c] -= 1
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
348
6
ea9d5fc4c001 "planemo upload commit 9ada186a78831ca2618ec817a23a77de6adf1a5d"
rhpvorderman
parents: 0
diff changeset
349 for k, cdr2Length in cdr2AALengthDic.items():
0
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
350 for c in absentAACDR2Dic[cdr2Length]:
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
351 absentAA[c] -= 1
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
352
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
353
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
354 aa_mutations_by_id_file = outfile[:outfile.rindex("/")] + "/absent_aa_id.txt"
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
355 with open(aa_mutations_by_id_file, 'w') as o:
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
356 o.write("ID\tcdr1length\tcdr2length\tbest_match\t" + "\t".join([str(x) for x in range(1,AALength)]) + "\n")
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
357 for ID in IDlist:
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
358 absentAAbyID = [1] * (AALength-1)
6
ea9d5fc4c001 "planemo upload commit 9ada186a78831ca2618ec817a23a77de6adf1a5d"
rhpvorderman
parents: 0
diff changeset
359 cdr1Length = cdr1AALengthDic[ID]
0
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
360 for c in absentAACDR1Dic[cdr1Length]:
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
361 absentAAbyID[c] -= 1
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
362
6
ea9d5fc4c001 "planemo upload commit 9ada186a78831ca2618ec817a23a77de6adf1a5d"
rhpvorderman
parents: 0
diff changeset
363 cdr2Length = cdr2AALengthDic[ID]
0
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
364 for c in absentAACDR2Dic[cdr2Length]:
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
365 absentAAbyID[c] -= 1
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
366 o.write(ID + "\t" + str(cdr1Length) + "\t" + str(cdr2Length) + "\t" + genedic[ID] + "\t" + "\t".join([str(x) for x in absentAAbyID]) + "\n")
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
367
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
368 if linecount == 0:
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
369 print("No data, exiting")
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
370 with open(outfile, 'w') as o:
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
371 o.write("RGYW (%)," + ("0,0,0\n" * len(genes)))
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
372 o.write("WRCY (%)," + ("0,0,0\n" * len(genes)))
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
373 o.write("WA (%)," + ("0,0,0\n" * len(genes)))
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
374 o.write("TW (%)," + ("0,0,0\n" * len(genes)))
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
375 sys.exit()
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
376
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
377 RGYWCount = {}
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
378 WRCYCount = {}
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
379 WACount = {}
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
380 TWCount = {}
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
381
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
382 #IDIndex = 0
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
383 ataIndex = 0
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
384 tatIndex = 0
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
385 aggctatIndex = 0
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
386 atagcctIndex = 0
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
387 first = True
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
388 with open(infile, 'r') as i:
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
389 for line in i:
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
390 if first:
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
391 linesplt = line.split("\t")
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
392 ataIndex = linesplt.index("X.a.t.a")
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
393 tatIndex = linesplt.index("t.a.t.")
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
394 aggctatIndex = linesplt.index("X.a.g.g.c.t..a.t.")
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
395 atagcctIndex = linesplt.index("X.a.t..a.g.c.c.t.")
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
396 first = False
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
397 continue
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
398 linesplt = line.split("\t")
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
399 gene = linesplt[best_matchIndex]
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
400 ID = linesplt[IDIndex]
6
ea9d5fc4c001 "planemo upload commit 9ada186a78831ca2618ec817a23a77de6adf1a5d"
rhpvorderman
parents: 0
diff changeset
401 RGYW = [Hotspot.from_string(x) for x in linesplt[aggctatIndex].split("|") if x]
ea9d5fc4c001 "planemo upload commit 9ada186a78831ca2618ec817a23a77de6adf1a5d"
rhpvorderman
parents: 0
diff changeset
402 WRCY = [Hotspot.from_string(x) for x in linesplt[atagcctIndex].split("|") if x]
ea9d5fc4c001 "planemo upload commit 9ada186a78831ca2618ec817a23a77de6adf1a5d"
rhpvorderman
parents: 0
diff changeset
403 WA = [Hotspot.from_string(x) for x in linesplt[ataIndex].split("|") if x]
ea9d5fc4c001 "planemo upload commit 9ada186a78831ca2618ec817a23a77de6adf1a5d"
rhpvorderman
parents: 0
diff changeset
404 TW = [Hotspot.from_string(x) for x in linesplt[tatIndex].split("|") if x]
0
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
405 RGYWCount[ID], WRCYCount[ID], WACount[ID], TWCount[ID] = 0, 0, 0, 0
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
406
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
407 with open(os.path.join(os.path.dirname(os.path.abspath(infile)), "RGYW.txt"), 'a') as out_handle:
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
408 for hotspot in RGYW:
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
409 out_handle.write("{0}\t{1}\n".format(ID, "\t".join([str(x) for x in hotspot])))
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
410
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
411 mutationList = mutationdic[ID + "_FR1"] + mutationdic[ID + "_CDR1"] + mutationdic[ID + "_FR2"] + mutationdic[ID + "_CDR2"] + mutationdic[ID + "_FR3"]
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
412 for mutation in mutationList:
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
413 frm, where, to, AAfrm, AAwhere, AAto, junk = mutation
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
414 mutation_in_RGYW = any(((start <= int(where) <= end) for (start, end, region) in RGYW))
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
415 mutation_in_WRCY = any(((start <= int(where) <= end) for (start, end, region) in WRCY))
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
416 mutation_in_WA = any(((start <= int(where) <= end) for (start, end, region) in WA))
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
417 mutation_in_TW = any(((start <= int(where) <= end) for (start, end, region) in TW))
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
418
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
419 in_how_many_motifs = sum([mutation_in_RGYW, mutation_in_WRCY, mutation_in_WA, mutation_in_TW])
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
420
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
421 if in_how_many_motifs > 0:
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
422 RGYWCount[ID] += (1.0 * int(mutation_in_RGYW)) / in_how_many_motifs
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
423 WRCYCount[ID] += (1.0 * int(mutation_in_WRCY)) / in_how_many_motifs
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
424 WACount[ID] += (1.0 * int(mutation_in_WA)) / in_how_many_motifs
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
425 TWCount[ID] += (1.0 * int(mutation_in_TW)) / in_how_many_motifs
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
426
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
427 mutations_in_motifs_file = os.path.join(os.path.dirname(os.path.abspath(infile)), "mutation_in_motifs.txt")
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
428 if not os.path.exists(mutation_by_id_file):
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
429 with open(mutations_in_motifs_file, 'w') as out_handle:
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
430 out_handle.write("{0}\n".format("\t".join([
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
431 "Sequence.ID",
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
432 "mutation_position",
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
433 "region",
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
434 "from_nt",
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
435 "to_nt",
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
436 "mutation_position_AA",
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
437 "from_AA",
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
438 "to_AA",
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
439 "motif",
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
440 "motif_start_nt",
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
441 "motif_end_nt",
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
442 "rest"
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
443 ])))
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
444
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
445 with open(mutations_in_motifs_file, 'a') as out_handle:
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
446 motif_dic = {"RGYW": RGYW, "WRCY": WRCY, "WA": WA, "TW": TW}
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
447 for mutation in mutationList:
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
448 frm, where, to, AAfrm, AAwhere, AAto, junk = mutation
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
449 for motif in list(motif_dic.keys()):
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
450
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
451 for start, end, region in motif_dic[motif]:
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
452 if start <= int(where) <= end:
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
453 out_handle.write("{0}\n".format(
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
454 "\t".join([
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
455 ID,
6
ea9d5fc4c001 "planemo upload commit 9ada186a78831ca2618ec817a23a77de6adf1a5d"
rhpvorderman
parents: 0
diff changeset
456 str(where),
0
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
457 region,
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
458 frm,
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
459 to,
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
460 str(AAwhere),
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
461 str(AAfrm),
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
462 str(AAto),
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
463 motif,
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
464 str(start),
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
465 str(end),
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
466 str(junk)
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
467 ])
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
468 ))
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
469
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
470
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
471
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
472 def mean(lst):
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
473 return (float(sum(lst)) / len(lst)) if len(lst) > 0 else 0.0
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
474
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
475
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
476 def median(lst):
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
477 lst = sorted(lst)
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
478 l = len(lst)
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
479 if l == 0:
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
480 return 0
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
481 if l == 1:
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
482 return lst[0]
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
483
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
484 l = int(l / 2)
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
485
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
486 if len(lst) % 2 == 0:
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
487 return float(lst[l] + lst[(l - 1)]) / 2.0
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
488 else:
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
489 return lst[l]
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
490
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
491 funcs = {"mean": mean, "median": median, "sum": sum}
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
492
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
493 directory = outfile[:outfile.rfind("/") + 1]
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
494 value = 0
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
495 valuedic = dict()
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
496
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
497 for fname in list(funcs.keys()):
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
498 for gene in genes:
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
499 with open(directory + gene + "_" + fname + "_value.txt", 'r') as v:
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
500 valuedic[gene + "_" + fname] = float(v.readlines()[0].rstrip())
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
501 with open(directory + "all_" + fname + "_value.txt", 'r') as v:
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
502 valuedic["total_" + fname] = float(v.readlines()[0].rstrip())
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
503
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
504
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
505 def get_xyz(lst, gene, f, fname):
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
506 x = round(round(f(lst), 1))
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
507 y = valuedic[gene + "_" + fname]
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
508 z = str(round(x / float(y) * 100, 1)) if y != 0 else "0"
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
509 return (str(x), str(y), z)
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
510
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
511 dic = {"RGYW": RGYWCount, "WRCY": WRCYCount, "WA": WACount, "TW": TWCount}
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
512 arr = ["RGYW", "WRCY", "WA", "TW"]
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
513
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
514 for fname in list(funcs.keys()):
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
515 func = funcs[fname]
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
516 foutfile = outfile[:outfile.rindex("/")] + "/hotspot_analysis_" + fname + ".txt"
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
517 with open(foutfile, 'w') as o:
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
518 for typ in arr:
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
519 o.write(typ + " (%)")
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
520 curr = dic[typ]
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
521 for gene in genes:
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
522 if valuedic[gene + "_" + fname] is 0:
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
523 o.write(",0,0,0")
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
524 else:
6
ea9d5fc4c001 "planemo upload commit 9ada186a78831ca2618ec817a23a77de6adf1a5d"
rhpvorderman
parents: 0
diff changeset
525 x, y, z = get_xyz([curr[x] for x in [y for y, z in genedic.items() if z.startswith(gene)]], gene, func, fname)
0
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
526 o.write("," + x + "," + y + "," + z)
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
527 x, y, z = get_xyz([y for x, y in curr.items() if not genedic[x].startswith("unmatched")], "total", func, fname)
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
528 #x, y, z = get_xyz([y for x, y in curr.iteritems()], "total", func, fname)
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
529 o.write("," + x + "," + y + "," + z + "\n")
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
530
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
531
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
532 # for testing
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
533 seq_motif_file = outfile[:outfile.rindex("/")] + "/motif_per_seq.txt"
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
534 with open(seq_motif_file, 'w') as o:
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
535 o.write("ID\tRGYW\tWRCY\tWA\tTW\n")
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
536 for ID in IDlist:
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
537 #o.write(ID + "\t" + str(round(RGYWCount[ID], 2)) + "\t" + str(round(WRCYCount[ID], 2)) + "\t" + str(round(WACount[ID], 2)) + "\t" + str(round(TWCount[ID], 2)) + "\n")
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
538 o.write(ID + "\t" + str(RGYWCount[ID]) + "\t" + str(WRCYCount[ID]) + "\t" + str(WACount[ID]) + "\t" + str(TWCount[ID]) + "\n")
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
539
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
540 if __name__ == "__main__":
64d74ba01a7c "planemo upload commit 78d1fae87dbcf490e49a9f99e7a06de7328e16d4"
rhpvorderman
parents:
diff changeset
541 main()