|
0
|
1 #!/usr/bin/env python3
|
|
|
2
|
|
|
3 import time
|
|
|
4 import configuration
|
|
|
5 import os
|
|
|
6 import textwrap
|
|
|
7 import subprocess
|
|
|
8 from tempfile import NamedTemporaryFile
|
|
|
9 from collections import defaultdict
|
|
|
10
|
|
|
11
|
|
|
12 class Range():
|
|
|
13 '''
|
|
|
14 This class is used to check float range in argparse
|
|
|
15 '''
|
|
|
16
|
|
|
17 def __init__(self, start, end):
|
|
|
18 self.start = start
|
|
|
19 self.end = end
|
|
|
20
|
|
|
21 def __eq__(self, other):
|
|
|
22 return self.start <= other <= self.end
|
|
|
23
|
|
|
24 def __str__(self):
|
|
|
25 return "float range {}..{}".format(self.start, self.end)
|
|
|
26
|
|
|
27 def __repr__(self):
|
|
|
28 return "float range {}..{}".format(self.start, self.end)
|
|
|
29
|
|
|
30
|
|
|
31 def check_file_start(gff_file):
|
|
|
32 count_comment = 0
|
|
|
33 with open(gff_file, "r") as gff_all:
|
|
|
34 line = gff_all.readline()
|
|
|
35 while line.startswith("#"):
|
|
|
36 line = gff_all.readline()
|
|
|
37 count_comment += 1
|
|
|
38 return count_comment
|
|
|
39
|
|
|
40
|
|
|
41 def write_info(filt_dom_tmp, FILT_DOM_GFF, orig_class_dict, filt_class_dict,
|
|
|
42 dom_dict, version_lines):
|
|
|
43 '''
|
|
|
44 Write domains statistics in beginning of filtered GFF
|
|
|
45 '''
|
|
|
46 with open(FILT_DOM_GFF, "w") as filt_gff:
|
|
|
47 for line in version_lines:
|
|
|
48 filt_gff.write(line)
|
|
|
49 filt_gff.write("##CLASSIFICATION\tORIGINAL_COUNTS\tFILTERED_COUNTS\n")
|
|
|
50 if not orig_class_dict:
|
|
|
51 filt_gff.write("##NO DOMAINS CLASSIFICATIONS\n")
|
|
|
52 for classification in sorted(orig_class_dict.keys()):
|
|
|
53 if classification in filt_class_dict.keys():
|
|
|
54 filt_gff.write("##{}\t{}\t{}\n".format(
|
|
|
55 classification, orig_class_dict[
|
|
|
56 classification], filt_class_dict[classification]))
|
|
|
57 else:
|
|
|
58 filt_gff.write("##{}\t{}\t{}\n".format(
|
|
|
59 classification, orig_class_dict[classification], 0))
|
|
|
60 filt_gff.write("##-----------------------------------------------\n"
|
|
|
61 "##SEQ\tDOMAIN\tCOUNTS\n")
|
|
|
62 if not dom_dict:
|
|
|
63 filt_gff.write("##NO DOMAINS\n")
|
|
|
64 for seq in sorted(dom_dict.keys()):
|
|
|
65 for dom, count in sorted(dom_dict[seq].items()):
|
|
|
66 filt_gff.write("##{}\t{}\t{}\n".format(seq, dom, count))
|
|
|
67 filt_gff.write("##-----------------------------------------------\n")
|
|
|
68 with open(filt_dom_tmp.name, "r") as filt_tmp:
|
|
|
69 for line in filt_tmp:
|
|
|
70 filt_gff.write(line)
|
|
|
71
|
|
|
72
|
|
|
73 def get_file_start(gff_file):
|
|
|
74 count_comment = 0
|
|
|
75 lines = []
|
|
|
76 with open(gff_file, "r") as gff_all:
|
|
|
77 line = gff_all.readline()
|
|
|
78 while line.startswith("#"):
|
|
|
79 lines.append(line)
|
|
|
80 line = gff_all.readline()
|
|
|
81 count_comment += 1
|
|
|
82 return count_comment, lines
|
|
|
83
|
|
|
84
|
|
15
|
85 def parse_gff_line(line):
|
|
|
86 '''Return dictionary with gff fields and atributers
|
|
|
87 Note - type of fields is strings
|
|
|
88 '''
|
|
|
89 # order of first 9 column is fixed
|
|
|
90 gff_line = dict(
|
|
|
91 zip(
|
|
|
92 ['seqid', 'source', 'type', 'start', 'end',
|
|
|
93 'score', 'strand', 'phase', 'attributes'],
|
|
|
94 line.split("\t")
|
|
|
95 )
|
|
|
96 )
|
|
|
97 # split attributes and replace:
|
|
|
98 gff_line['attributes'] = dict([i.split("=") for i in gff_line['attributes'].split(";")])
|
|
|
99 return gff_line
|
|
|
100
|
|
0
|
101 def filter_qual_dom(DOM_GFF, FILT_DOM_GFF, TH_IDENTITY, TH_SIMILARITY,
|
|
|
102 TH_LENGTH, TH_INTERRUPT, TH_LEN_RATIO, SELECTED_DOM,
|
|
|
103 ELEMENT):
|
|
|
104 ''' Filter gff output based on domain and quality of alignment '''
|
|
|
105 [count_comment, version_lines] = get_file_start(DOM_GFF)
|
|
|
106 filt_dom_tmp = NamedTemporaryFile(delete=False)
|
|
|
107 with open(DOM_GFF, "r") as gff_all, open(filt_dom_tmp.name,
|
|
|
108 "w") as gff_filtered:
|
|
15
|
109 for _ in range(count_comment):
|
|
0
|
110 next(gff_all)
|
|
|
111 dom_dict = defaultdict(lambda: defaultdict(int))
|
|
|
112 orig_class_dict = defaultdict(int)
|
|
|
113 filt_class_dict = defaultdict(int)
|
|
|
114 seq_ids_all = []
|
|
|
115 xminimals = []
|
|
|
116 xmaximals = []
|
|
|
117 domains = []
|
|
|
118 xminimals_all = []
|
|
|
119 xmaximals_all = []
|
|
|
120 domains_all = []
|
|
|
121 start = True
|
|
|
122 for line in gff_all:
|
|
17
|
123 gff_line = parse_gff_line(line)
|
|
|
124 classification = gff_line['attributes']['Final_Classification']
|
|
0
|
125 orig_class_dict[classification] += 1
|
|
|
126 ## ambiguous domains filtered out automatically
|
|
|
127 if classification != configuration.AMBIGUOUS_TAG:
|
|
15
|
128 al_identity = float(gff_line['attributes']['Identity'])
|
|
|
129 al_similarity = float(gff_line['attributes']['Similarity'])
|
|
|
130 al_length = float(gff_line['attributes']['Relat_Length'])
|
|
|
131 relat_interrupt = float(gff_line['attributes']['Relat_Interruptions'])
|
|
|
132 db_len_proportion = float(gff_line['attributes']['Hit_to_DB_Length'])
|
|
17
|
133 dom_type = gff_line['attributes']['Name']
|
|
15
|
134 seq_id = gff_line['seqid']
|
|
|
135 xminimal = int(gff_line['start'])
|
|
|
136 xmaximal = int(gff_line['end'])
|
|
|
137 c1 = al_identity >= TH_IDENTITY
|
|
|
138 c2 = al_similarity >= TH_SIMILARITY
|
|
|
139 if (c1 and c2 and al_length >= TH_LENGTH and relat_interrupt <= TH_INTERRUPT and
|
|
|
140 db_len_proportion <= TH_LEN_RATIO and
|
|
|
141 (dom_type == SELECTED_DOM or SELECTED_DOM == "All") and
|
|
|
142 (ELEMENT in classification)):
|
|
0
|
143 gff_filtered.writelines(line)
|
|
|
144 filt_class_dict[classification] += 1
|
|
|
145 dom_dict[seq_id][dom_type] += 1
|
|
|
146 if start:
|
|
|
147 seq_ids_all.append(line.split("\t")[0])
|
|
|
148 start = False
|
|
|
149 if seq_id != seq_ids_all[-1]:
|
|
|
150 seq_ids_all.append(seq_id)
|
|
|
151 xminimals_all.append(xminimals)
|
|
|
152 xmaximals_all.append(xmaximals)
|
|
|
153 domains_all.append(domains)
|
|
|
154 xminimals = []
|
|
|
155 xmaximals = []
|
|
|
156 domains = []
|
|
|
157 xminimals.append(xminimal)
|
|
|
158 xmaximals.append(xmaximal)
|
|
|
159 domains.append(dom_type)
|
|
|
160 path = os.path.dirname(os.path.realpath(__file__))
|
|
|
161 write_info(filt_dom_tmp, FILT_DOM_GFF, orig_class_dict, filt_class_dict,
|
|
|
162 dom_dict, version_lines)
|
|
|
163 os.unlink(filt_dom_tmp.name)
|
|
|
164 xminimals_all.append(xminimals)
|
|
|
165 xmaximals_all.append(xmaximals)
|
|
|
166 domains_all.append(domains)
|
|
|
167 return xminimals_all, xmaximals_all, domains_all, seq_ids_all
|
|
|
168
|
|
|
169
|
|
|
170 def get_domains_protseq(FILT_DOM_GFF, DOMAIN_PROT_SEQ):
|
|
|
171 ''' Get the translated protein sequence of original DNA seq for all the filtered domains regions
|
|
|
172 The translated sequences are taken from alignment reported by LASTAL (Query_Seq attribute in GFF)
|
|
|
173 '''
|
|
|
174 count_comment = check_file_start(FILT_DOM_GFF)
|
|
|
175 with open(FILT_DOM_GFF, "r") as filt_gff:
|
|
|
176 for comment_idx in range(count_comment):
|
|
|
177 next(filt_gff)
|
|
|
178 with open(DOMAIN_PROT_SEQ, "w") as dom_prot_file:
|
|
|
179 for line in filt_gff:
|
|
|
180 attributes = line.rstrip().split("\t")[8]
|
|
|
181 positions = attributes.split(";")[3].split("=")[1].split(":")[
|
|
|
182 -1].split("[")[0]
|
|
|
183 dom = attributes.split(";")[0].split("=")[1]
|
|
|
184 dom_class = attributes.split(";")[1].split("=")[1]
|
|
|
185 seq_id = line.rstrip().split("\t")[0]
|
|
|
186 prot_seq_align = line.rstrip().split("\t")[8].split(";")[
|
|
|
187 6].split("=")[1]
|
|
|
188 prot_seq = prot_seq_align.translate({ord(i): None
|
|
|
189 for i in '/\\-'})
|
|
|
190 header_prot_seq = ">{}:{} {} {}".format(seq_id, positions, dom,
|
|
|
191 dom_class)
|
|
|
192 dom_prot_file.write("{}\n{}\n".format(
|
|
|
193 header_prot_seq, textwrap.fill(prot_seq,
|
|
|
194 configuration.FASTA_LINE)))
|
|
|
195
|
|
|
196
|
|
|
197 def main(args):
|
|
|
198
|
|
|
199 t = time.time()
|
|
|
200
|
|
|
201 DOM_GFF = args.dom_gff
|
|
|
202 DOMAIN_PROT_SEQ = args.domains_prot_seq
|
|
|
203 TH_IDENTITY = args.th_identity
|
|
|
204 TH_LENGTH = args.th_length
|
|
|
205 TH_INTERRUPT = args.interruptions
|
|
|
206 TH_SIMILARITY = args.th_similarity
|
|
|
207 TH_LEN_RATIO = args.max_len_proportion
|
|
|
208 FILT_DOM_GFF = args.domains_filtered
|
|
|
209 SELECTED_DOM = args.selected_dom
|
|
|
210 OUTPUT_DIR = args.output_dir
|
|
|
211 # DELETE : ELEMENT = args.element_type.replace("_pipe_", "|")
|
|
|
212 ELEMENT = args.element_type
|
|
|
213
|
|
|
214 if DOMAIN_PROT_SEQ is None:
|
|
|
215 DOMAIN_PROT_SEQ = configuration.DOM_PROT_SEQ
|
|
|
216 if FILT_DOM_GFF is None:
|
|
|
217 FILT_DOM_GFF = configuration.FILT_DOM_GFF
|
|
|
218
|
|
|
219 if OUTPUT_DIR and not os.path.exists(OUTPUT_DIR):
|
|
|
220 os.makedirs(OUTPUT_DIR)
|
|
|
221
|
|
|
222 if not os.path.isabs(FILT_DOM_GFF):
|
|
|
223 if OUTPUT_DIR is None:
|
|
|
224 OUTPUT_DIR = os.path.dirname(os.path.abspath(DOM_GFF))
|
|
|
225 FILT_DOM_GFF = os.path.join(OUTPUT_DIR, os.path.basename(FILT_DOM_GFF))
|
|
|
226 DOMAIN_PROT_SEQ = os.path.join(OUTPUT_DIR,
|
|
|
227 os.path.basename(DOMAIN_PROT_SEQ))
|
|
|
228
|
|
|
229 [xminimals_all, xmaximals_all, domains_all, seq_ids_all] = filter_qual_dom(
|
|
|
230 DOM_GFF, FILT_DOM_GFF, TH_IDENTITY, TH_SIMILARITY, TH_LENGTH,
|
|
|
231 TH_INTERRUPT, TH_LEN_RATIO, SELECTED_DOM, ELEMENT)
|
|
|
232 get_domains_protseq(FILT_DOM_GFF, DOMAIN_PROT_SEQ)
|
|
|
233
|
|
|
234 print("ELAPSED_TIME_DOMAINS = {} s".format(time.time() - t))
|
|
|
235
|
|
|
236
|
|
|
237 if __name__ == "__main__":
|
|
|
238 import argparse
|
|
|
239 from argparse import RawDescriptionHelpFormatter
|
|
|
240
|
|
|
241 class CustomFormatter(argparse.ArgumentDefaultsHelpFormatter,
|
|
|
242 argparse.RawDescriptionHelpFormatter):
|
|
|
243 pass
|
|
|
244
|
|
|
245 parser = argparse.ArgumentParser(
|
|
|
246 description=
|
|
|
247 '''The script performs DANTE's output filtering for quality and/or extracting specific type of protein domain or mobile elements of origin. For the filtered domains it reports their translated protein sequence of original DNA.
|
|
|
248 WHEN NO PARAMETERS GIVEN, IT PERFORMS QUALITY FILTERING USING THE DEFAULT PARAMETRES (optimized for Viridiplantae species)
|
|
|
249
|
|
|
250 INPUTS:
|
|
|
251 - GFF3 file produced by protein_domains.py OR already filtered GFF3
|
|
|
252
|
|
|
253 FILTERING OPTIONS:
|
|
|
254 > QUALITY: - Min relative length of alignemnt to the protein domain from DB (without gaps)
|
|
|
255 - Identity
|
|
|
256 - Similarity (scoring matrix: BLOSUM82)
|
|
|
257 - Interruption in the reading frame (frameshifts + stop codons) per every starting 100 AA
|
|
|
258 - Max alignment proportion to the original length of database domain sequence
|
|
|
259 > DOMAIN TYPE: choose from choices ('Name' attribute in GFF)
|
|
|
260 Records for ambiguous domain type (e.g. INT/RH) are filtered out automatically
|
|
|
261
|
|
|
262 > MOBILE ELEMENT TYPE:
|
|
|
263 arbitrary substring of the element classification ('Final_Classification' attribute in GFF)
|
|
|
264
|
|
|
265 OUTPUTS:
|
|
|
266 - filtered GFF3 file
|
|
|
267 - fasta file of translated protein sequences (from original DNA) for the aligned domains that match the filtering criteria
|
|
|
268
|
|
|
269 DEPENDENCIES:
|
|
|
270 - python 3.4 or higher
|
|
|
271 > ProfRep modules:
|
|
|
272 - configuration.py
|
|
|
273
|
|
|
274 EXAMPLE OF USAGE:
|
|
|
275 Getting quality filtered integrase(INT) domains of all gypsy transposable elements:
|
|
|
276 ./domains_filtering.py -dom_gff PATH_TO_INPUT_GFF -pdb PATH_TO_PROTEIN_DB -cs PATH_TO_CLASSIFICATION_FILE --selected_dom INT --element_type Ty3/gypsy
|
|
|
277
|
|
|
278 ''',
|
|
|
279 epilog="""""",
|
|
|
280 formatter_class=CustomFormatter)
|
|
|
281 requiredNamed = parser.add_argument_group('required named arguments')
|
|
|
282 requiredNamed.add_argument("-dg",
|
|
|
283 "--dom_gff",
|
|
|
284 type=str,
|
|
|
285 required=True,
|
|
|
286 help="basic unfiltered gff file of all domains")
|
|
|
287 parser.add_argument("-ouf",
|
|
|
288 "--domains_filtered",
|
|
|
289 type=str,
|
|
|
290 help="output filtered domains gff file")
|
|
|
291 parser.add_argument("-dps",
|
|
|
292 "--domains_prot_seq",
|
|
|
293 type=str,
|
|
|
294 help="output file containg domains protein sequences")
|
|
|
295 parser.add_argument("-thl",
|
|
|
296 "--th_length",
|
|
|
297 type=float,
|
|
|
298 choices=[Range(0.0, 1.0)],
|
|
|
299 default=0.8,
|
|
|
300 help="proportion of alignment length threshold")
|
|
|
301 parser.add_argument("-thi",
|
|
|
302 "--th_identity",
|
|
|
303 type=float,
|
|
|
304 choices=[Range(0.0, 1.0)],
|
|
|
305 default=0.35,
|
|
|
306 help="proportion of alignment identity threshold")
|
|
|
307 parser.add_argument("-ths",
|
|
|
308 "--th_similarity",
|
|
|
309 type=float,
|
|
|
310 choices=[Range(0.0, 1.0)],
|
|
|
311 default=0.45,
|
|
|
312 help="threshold for alignment proportional similarity")
|
|
|
313 parser.add_argument(
|
|
|
314 "-ir",
|
|
|
315 "--interruptions",
|
|
|
316 type=int,
|
|
|
317 default=3,
|
|
|
318 help=
|
|
|
319 "interruptions (frameshifts + stop codons) tolerance threshold per 100 AA")
|
|
|
320 parser.add_argument(
|
|
|
321 "-mlen",
|
|
|
322 "--max_len_proportion",
|
|
|
323 type=float,
|
|
|
324 default=1.2,
|
|
|
325 help=
|
|
|
326 "maximal proportion of alignment length to the original length of protein domain from database")
|
|
|
327 parser.add_argument(
|
|
|
328 "-sd",
|
|
|
329 "--selected_dom",
|
|
|
330 type=str,
|
|
|
331 default="All",
|
|
|
332 choices=[
|
|
|
333 "All", "GAG", "INT", "PROT", "RH", "RT", "aRH", "CHDCR", "CHDII",
|
|
|
334 "TPase", "YR", "HEL1", "HEL2", "ENDO"
|
|
|
335 ],
|
|
|
336 help="filter output domains based on the domain type")
|
|
|
337 parser.add_argument(
|
|
|
338 "-el",
|
|
|
339 "--element_type",
|
|
|
340 type=str,
|
|
|
341 default="",
|
|
|
342 help="filter output domains by typing substring from classification")
|
|
|
343 parser.add_argument(
|
|
|
344 "-dir",
|
|
|
345 "--output_dir",
|
|
|
346 type=str,
|
|
|
347 default=None,
|
|
|
348 help="specify if you want to change the output directory")
|
|
|
349 args = parser.parse_args()
|
|
|
350 main(args)
|
