|
0
|
1 #!/usr/bin/env python3
|
|
|
2 ''' TAndem REpeat ANalyzer '''
|
|
|
3 import os
|
|
|
4 import sys
|
|
|
5 import shutil
|
|
|
6 import subprocess
|
|
|
7 import argparse
|
|
|
8 from argparse import RawTextHelpFormatter
|
|
|
9 import logging
|
|
|
10 import shlex
|
|
|
11 import multiprocessing
|
|
|
12 # config must be loaded before seqtools,...
|
|
|
13 import config
|
|
|
14 import re
|
|
|
15 from lib import seqtools, graphtools, utils, assembly_tools
|
|
|
16 from lib import r2py
|
|
|
17
|
|
|
18 REQUIRED_VERSION = (3, 4)
|
|
|
19 if sys.version_info < REQUIRED_VERSION:
|
|
|
20 raise Exception("\n\npython 3.4 or higher is required!\n")
|
|
|
21
|
|
|
22 # append path to louvain clustering and other binaries
|
|
|
23 os.environ['PATH'] = "{}:{}:{}".format(config.BINARIES, config.LOUVAIN,
|
|
|
24 os.environ['PATH'])
|
|
|
25
|
|
|
26 LOGGER = logging.getLogger(__name__)
|
|
|
27
|
|
|
28
|
|
|
29 def get_version(path, tarean_mode):
|
|
|
30 try:
|
|
|
31 branch = subprocess.check_output("git rev-parse --abbrev-ref HEAD",
|
|
|
32 shell=True,
|
|
|
33 cwd=path).decode('ascii').strip()
|
|
|
34 shorthash = subprocess.check_output(
|
|
|
35 "git log --pretty=format:'%h' -n 1 ",
|
|
|
36 shell=True,
|
|
|
37 cwd=path).decode('ascii').strip()
|
|
|
38 revcount = len(subprocess.check_output(
|
|
|
39 "git log --oneline", shell=True,
|
|
|
40 cwd=path).decode('ascii').split())
|
|
|
41 try:
|
|
|
42 tag = subprocess.check_output("git describe --tags --abbrev=0",
|
|
|
43 cwd=path,
|
|
|
44 shell=True).decode('ascii').strip()
|
|
|
45 except subprocess.CalledProcessError:
|
|
|
46 tag = " "
|
|
|
47
|
|
|
48 version_string = (
|
|
|
49 "-------------------------------------"
|
|
|
50 "-------------------------------------\n"
|
|
|
51 "PIPELINE VERSION : "
|
|
|
52 "{branch}-{tag}-{revcount}({shorthash})\n\n"
|
|
|
53 "PROTEIN DATABASE VERSION : {PD}\n"
|
|
|
54 " md5 checksum : {PDmd5}\n\n"
|
|
|
55 "DNA DATABASE VERSION : {DD}\n"
|
|
|
56 " md5 checksum : {DDmd5}\n"
|
|
|
57 "-------------------------------------"
|
|
|
58 "-------------------------------------\n").format(
|
|
|
59 branch=branch,
|
|
|
60 shorthash=shorthash,
|
|
|
61 revcount=revcount,
|
|
|
62 tag=tag,
|
|
|
63 PD=os.path.basename(config.PROTEIN_DATABASE),
|
|
|
64 PDmd5=utils.md5checksum(config.PROTEIN_DATABASE + ".psq", fail_if_missing = not tarean_mode),
|
|
|
65 DD=os.path.basename(config.DNA_DATABASE),
|
|
|
66 DDmd5=utils.md5checksum(config.DNA_DATABASE + ".nsq"))
|
|
|
67
|
|
1
|
68 except:
|
|
0
|
69 version_string = "version of pipeline not available!"
|
|
1
|
70
|
|
0
|
71 LOGGER.info(version_string)
|
|
|
72 return version_string
|
|
|
73
|
|
|
74
|
|
|
75 def valid_database(database_file):
|
|
|
76 with open(database_file, 'r', encoding='ascii') as f:
|
|
|
77 for i in f:
|
|
|
78 if i[0] == ">":
|
|
|
79 if not re.match(">.+#.+/*", i):
|
|
|
80 # TODO - make edits to correct fomating of custom database???
|
|
|
81 return False
|
|
|
82 return True
|
|
|
83
|
|
|
84
|
|
|
85 def add_databases(databases, custom_databases_dir, dbtype='nucl'):
|
|
|
86 '''custom databases are copied to directory tree and blast
|
|
|
87 database is created using makeblastdb
|
|
|
88 '''
|
|
|
89
|
|
|
90 databases_ok = []
|
|
|
91 print(databases)
|
|
|
92 for db_path, db_name in databases:
|
|
|
93 db_destination = "{}/{}".format(custom_databases_dir, db_name)
|
|
|
94 shutil.copyfile(db_path, db_destination)
|
|
|
95 if not valid_database(db_destination):
|
|
|
96 raise ValueError((
|
|
|
97 "\n------------------------------------------------------------\n"
|
|
|
98 "Custom database is not valid!\n"
|
|
|
99 "Custom database of repeats are DNA sequences in fasta format.\n"
|
|
|
100 "The required format for IDs in a custom library is : \n"
|
|
|
101 " '>reapeatname#class/subclass'\n"
|
|
|
102 "Reformat the database and try again!\n"
|
|
|
103 "-------------------------------------------------------------\n\n"
|
|
|
104 ))
|
|
|
105
|
|
|
106 cmd = "makeblastdb -in {0} -out {0} -dbtype {1}".format(db_destination,
|
|
|
107 dbtype)
|
|
|
108 print(cmd)
|
|
|
109 args = shlex.split(cmd)
|
|
|
110 print(args)
|
|
|
111 if subprocess.check_call(args, stderr=sys.stdout):
|
|
|
112 Warning("makeblastdb on {} failed".format(db_name))
|
|
|
113 else:
|
|
|
114 databases_ok.append([db_destination, "custom_db_" + db_name])
|
|
|
115 if len(databases_ok) == 0:
|
|
|
116 return None
|
|
|
117 else:
|
|
|
118 return databases_ok
|
|
|
119
|
|
|
120
|
|
|
121 def meminfo():
|
|
|
122 ''' detect physical memory and memory usage'''
|
|
|
123 info = {}
|
|
|
124 required_fields = [
|
|
|
125 'MemTotal:', 'MemFree:', 'Cached:', 'SwapCached:', 'Buffers:'
|
|
|
126 ]
|
|
|
127 with open('/proc/meminfo', 'r') as f:
|
|
|
128 for i in f:
|
|
|
129 a = i.split()
|
|
|
130 if a[0] in required_fields:
|
|
|
131 info[a[0]] = int(a[1])
|
|
|
132 return info
|
|
|
133
|
|
|
134
|
|
|
135 def dict2lists(d):
|
|
|
136 ''' convert dict to nested list
|
|
|
137 use the funsction to pass dictionary to R function
|
|
|
138 '''
|
|
|
139 values = list(d.values())
|
|
|
140 keys = list(d.keys())
|
|
|
141 return [values, keys]
|
|
|
142
|
|
|
143
|
|
|
144 def show_object(obj):
|
|
|
145 '''
|
|
|
146 helper function for printing all public atributes,
|
|
|
147 does not print callebme atributes e.i. methods..
|
|
|
148 '''
|
|
|
149
|
|
|
150 s = "Configuration--------------->\n"
|
|
|
151 for i in dir(obj):
|
|
|
152 # do not show private
|
|
|
153 if i[:2] != "__":
|
|
|
154 value = getattr(obj, i)
|
|
|
155 if not callable(value):
|
|
|
156 s += "{} : {}\n".format(i, value)
|
|
|
157 s += "<---------------configuration\n"
|
|
|
158 return s
|
|
|
159
|
|
|
160
|
|
|
161 class DataInfo():
|
|
|
162 '''
|
|
|
163 stores information state of clustering and data
|
|
|
164 '''
|
|
|
165
|
|
|
166 def __init__(self, args, paths):
|
|
|
167 LOGGER.info("getting information about input sequences")
|
|
|
168 self.args = args
|
|
|
169 self.working_directory = args.output_dir
|
|
|
170 self.input_sequences = args.sequences.name
|
|
|
171 self.number_of_input_sequences = seqtools.SequenceSet.fasta_length(
|
|
|
172 self.input_sequences)
|
|
|
173 self.paired = args.paired
|
|
|
174 self.prefix_length = args.prefix_length
|
|
|
175 self.physical_memory = meminfo()['MemTotal:']
|
|
|
176 self.edges_max = config.EMAX
|
|
|
177 # set max memory
|
|
|
178 if args.max_memory:
|
|
|
179 self.max_memory = args.max_memory
|
|
|
180 else:
|
|
|
181 self.max_memory = meminfo()["MemTotal:"]
|
|
|
182 # modify initial setup if number of sequences is low
|
|
|
183 if args.automatic_filtering:
|
|
|
184 config.NUMBER_OF_SEQUENCES_FOR_PRERUN = config.NUMBER_OF_SEQUENCES_FOR_PRERUN_WITH_FILTERING
|
|
|
185
|
|
|
186 if self.number_of_input_sequences < config.NUMBER_OF_SEQUENCES_FOR_PRERUN:
|
|
|
187 config.NUMBER_OF_SEQUENCES_FOR_PRERUN = self.number_of_input_sequences
|
|
|
188
|
|
|
189 # is number of input sequences sufficient
|
|
|
190 if self.number_of_input_sequences < config.MINIMUM_NUMBER_OF_INPUT_SEQUENCES:
|
|
|
191 raise WrongInputDataError(
|
|
|
192 "provide more sequences for clustering, minumum {} is .required".format(
|
|
|
193 config.MINIMUM_NUMBER_OF_INPUT_SEQUENCES))
|
|
|
194 # these atribudes will be set later after clustering is done
|
|
|
195 self.max_annotated_clusters = None
|
|
|
196 self.max_annotated_superclusters = None
|
|
|
197 # the atributes will be set after prerun is performed
|
|
|
198 self.prerun_ecount = None
|
|
|
199 self.prerun_ecount_corrected = None
|
|
|
200 self.sample_size = None
|
|
|
201 self.max_number_reads_for_clustering = None
|
|
|
202 self.mincln = None
|
|
|
203 self.number_of_omitted_reads = 0
|
|
|
204 LOGGER.info("sampling sequences for prerun analysis")
|
|
|
205 sample = seqtools.SequenceSet(
|
|
|
206 source=self.input_sequences,
|
|
|
207 sample_size=config.NUMBER_OF_SEQUENCES_FOR_PRERUN,
|
|
|
208 paired=self.paired,
|
|
|
209 filename=paths.sample_db,
|
|
|
210 fasta=paths.sample_fasta,
|
|
|
211 rename=True)
|
|
|
212 sample.makeblastdb(legacy=args.options.legacy_database, lastdb=args.options.lastdb)
|
|
|
213 # preliminary clustering
|
|
|
214 self.prerun_vcount = len(sample)
|
|
|
215 # line count
|
|
|
216 self._prerun(sample, paths)
|
|
|
217 # adjust size of chunks:
|
|
|
218 if self.number_of_reads_for_clustering < config.CHUNK_SIZE * 30:
|
|
|
219 config.CHUNK_SIZE = round(self.number_of_reads_for_clustering / 40)
|
|
|
220
|
|
|
221 def _prerun(self, sample, paths):
|
|
|
222 '''Preliminary characterization sequences using
|
|
|
223 clustering on small dataset - stored as sample '''
|
|
|
224 sample.make_chunks(chunk_size=1000)
|
|
|
225 sample.create_hitsort(options=self.args.options)
|
|
|
226 sample_hitsort = graphtools.Graph(source=sample.hitsort,
|
|
|
227 paired=self.paired,
|
|
|
228 seqids=sample.keys())
|
|
|
229 sample_hitsort.save_indexed_graph()
|
|
|
230 sample_hitsort.louvain_clustering(merge_threshold=0.2)
|
|
|
231 sample_hitsort.export_cls(path=paths.prerun_cls_file)
|
|
|
232 sample.annotate(
|
|
|
233 config.DNA_DATABASE,
|
|
|
234 annotation_name="dna_database",
|
|
|
235 directory=paths.prerun,
|
|
|
236 params=self.args.options.annotation_search_params.blastn)
|
|
|
237
|
|
|
238 selected_tarean_contigs = []
|
|
|
239 ecount_corrected = sample_hitsort.ecount
|
|
|
240 vcount_corrected = sample_hitsort.vcount
|
|
|
241 if self.args.automatic_filtering:
|
|
|
242 prerun_cluster_info = sample_hitsort.export_clusters_files_multiple(
|
|
|
243 min_size=10,
|
|
|
244 directory=paths.prerun_clusters,
|
|
|
245 sequences=sample,
|
|
|
246 tRNA_database_path=config.TRNA_DATABASE,
|
|
|
247 satellite_model_path=config.SATELLITE_MODEL)
|
|
|
248 # check of prerun contain clusters with large number of edges
|
|
|
249 # these sequences can be used for filtering
|
|
|
250 for cl in prerun_cluster_info:
|
|
|
251 print(cl.ecount, cl.vcount, sample_hitsort.ecount,
|
|
|
252 cl.tandem_rank)
|
|
|
253
|
|
|
254 if (cl.tandem_rank in config.TANDEM_RANKS[0:2] and
|
|
|
255 cl.ecount / sample_hitsort.ecount >
|
|
|
256 config.FILTER_MIN_PROP_THRESHOLD and
|
|
|
257 cl.vcount > config.FILTER_MIN_SIZE_THRESHOLD):
|
|
|
258 selected_tarean_contigs.append(cl.tarean_contig_file)
|
|
|
259 ecount_corrected -= cl.ecount
|
|
|
260 vcount_corrected -= cl.vcount
|
|
|
261
|
|
|
262 if selected_tarean_contigs:
|
|
|
263 with open(paths.filter_sequences_file, 'w') as out:
|
|
|
264 for fname in selected_tarean_contigs:
|
|
|
265 with open(fname, 'r') as f:
|
|
|
266 out.write(f.read())
|
|
|
267 self.sequence_fiter = paths.filter_sequences_file
|
|
|
268 else:
|
|
|
269 self.sequence_fiter = None
|
|
|
270
|
|
|
271 self.prerun_ecount = sample_hitsort.ecount
|
|
|
272 self.prerun_ecount_corrected = ecount_corrected
|
|
|
273 self.prerun_vcount_corrected = vcount_corrected
|
|
|
274 self.max_number_reads_for_clustering = round((
|
|
|
275 ((self.edges_max * self.max_memory) /
|
|
|
276 self.prerun_ecount_corrected * self.prerun_vcount**2)**(0.5)) / 2)
|
|
|
277
|
|
|
278 if self.max_number_reads_for_clustering >= self.number_of_input_sequences:
|
|
|
279 self.sample_size = 0
|
|
|
280 else:
|
|
|
281 self.sample_size = self.max_number_reads_for_clustering
|
|
|
282
|
|
|
283 n1 = self.sample_size if self.sample_size != 0 else self.number_of_input_sequences
|
|
|
284 n2 = self.args.sample if self.args.sample != 0 else self.number_of_input_sequences
|
|
|
285 self.number_of_reads_for_clustering = min(n1, n2)
|
|
|
286 # minlcn is set either based on mincl or value specified in config,
|
|
|
287 # whatever is higher
|
|
|
288 self.mincln = int(self.number_of_reads_for_clustering *
|
|
|
289 self.args.mincl / 100)
|
|
|
290 if self.mincln < config.MINIMUM_NUMBER_OF_READS_IN_CLUSTER:
|
|
|
291 self.mincln = config.MINIMUM_NUMBER_OF_READS_IN_CLUSTER
|
|
|
292
|
|
|
293 def __str__(self):
|
|
|
294 s = "Data info------------------->\n"
|
|
|
295 for i in dir(self):
|
|
|
296 # do not show private
|
|
|
297 if i[:2] != "__":
|
|
|
298 value = getattr(self, i)
|
|
|
299 if not callable(value):
|
|
|
300 s += "{} : {}\n".format(i, value)
|
|
|
301 s += "<----------------------Data info\n"
|
|
|
302 return s
|
|
|
303
|
|
|
304
|
|
|
305 class DataFiles(object):
|
|
|
306 '''
|
|
|
307 stores location of data files and create directories ...
|
|
|
308 atributes are:
|
|
|
309 - individual directories
|
|
|
310 - individual files
|
|
|
311 - list of files or directories
|
|
|
312
|
|
|
313 directories are created if does not exist
|
|
|
314 '''
|
|
|
315
|
|
|
316 def __init__(self, working_dir, subdirs, files):
|
|
|
317 LOGGER.info("creating directory structure")
|
|
|
318 self.working_dir = working_dir
|
|
|
319 # add and create directories paths
|
|
|
320 for i in subdirs:
|
|
|
321 d = os.path.join(self.working_dir, subdirs[i])
|
|
|
322 os.makedirs(d, exist_ok=True)
|
|
|
323 setattr(self, i, d)
|
|
|
324 setattr(self, i + "__relative", subdirs[i])
|
|
|
325 # add file paths
|
|
|
326 for i in files:
|
|
|
327 d = os.path.join(self.working_dir, files[i])
|
|
|
328 setattr(self, i, d)
|
|
|
329 setattr(self, i + "__relative", files[i])
|
|
|
330
|
|
|
331 def __str__(self):
|
|
|
332 s = ""
|
|
|
333 for i in dir(self):
|
|
|
334 # do not show private
|
|
|
335 if i[:2] != "__":
|
|
|
336 value = getattr(self, i)
|
|
|
337 if not callable(value):
|
|
|
338 s += "{} : {}\n".format(i, value)
|
|
|
339 return s
|
|
|
340
|
|
|
341 def as_list(self):
|
|
|
342 '''
|
|
|
343 convert attr and vaues to list - suitable for passing values to R functions
|
|
|
344 '''
|
|
|
345 values = list()
|
|
|
346 keys = list()
|
|
|
347 for i in dir(self):
|
|
|
348 # do not show private
|
|
|
349 if i[:2] != "__":
|
|
|
350 value = getattr(self, i)
|
|
|
351 if not callable(value):
|
|
|
352 values.append(value)
|
|
|
353 keys.append(i)
|
|
|
354 return [values, keys]
|
|
|
355
|
|
|
356 def cleanup(self, paths):
|
|
|
357 ''' will remove unnecessary files from working directory '''
|
|
|
358 for i in paths:
|
|
|
359 fn = getattr(self, i)
|
|
|
360 if os.path.exists(fn):
|
|
|
361 if os.path.isdir(fn):
|
|
|
362 shutil.rmtree(fn, ignore_errors=False)
|
|
|
363 else:
|
|
|
364 os.remove(fn)
|
|
|
365
|
|
|
366
|
|
|
367 class WrongInputDataError(Exception):
|
|
|
368 '''Custom exception for wrong input
|
|
|
369 '''
|
|
|
370
|
|
|
371 def __init__(self, arg):
|
|
|
372 super(WrongInputDataError, self).__init__(arg)
|
|
|
373 self.msg = arg
|
|
|
374
|
|
|
375
|
|
|
376 class Range():
|
|
|
377 '''
|
|
|
378 This class is used to check float range in argparse
|
|
|
379 '''
|
|
|
380
|
|
|
381 def __init__(self, start, end):
|
|
|
382 self.start = start
|
|
|
383 self.end = end
|
|
|
384
|
|
|
385 def __eq__(self, other):
|
|
|
386 return self.start <= other <= self.end
|
|
|
387
|
|
|
388 def __str__(self):
|
|
|
389 return "float range {}..{}".format(self.start, self.end)
|
|
|
390
|
|
|
391 def __repr__(self):
|
|
|
392 return "float range {}..{}".format(self.start, self.end)
|
|
|
393
|
|
|
394
|
|
|
395 class DirectoryType(object):
|
|
|
396 '''
|
|
|
397 this class is similar to argparse.FileType
|
|
|
398 for mode 'w' creates and check the access to the directory
|
|
|
399 for mode 'r' check the presence of the dictory and accesibility
|
|
|
400 '''
|
|
|
401
|
|
|
402 def __init__(self, mode='r'):
|
|
|
403 self._mode = mode
|
|
|
404
|
|
|
405 def __call__(self, string):
|
|
|
406 if self._mode == 'w':
|
|
|
407 try:
|
|
|
408 os.makedirs(string, exist_ok=True)
|
|
|
409 except FileExistsError:
|
|
|
410 raise argparse.ArgumentTypeError(
|
|
|
411 "Cannot create directory, '{}' is a file".format(string))
|
|
|
412 if os.access(string, os.W_OK):
|
|
|
413 return string
|
|
|
414 else:
|
|
|
415 raise argparse.ArgumentTypeError(
|
|
|
416 "Directory '{}' is not writable".format(string))
|
|
|
417 if self._mode == 'r':
|
|
|
418 if not os.path.isdir(string):
|
|
|
419 raise argparse.ArgumentTypeError(
|
|
|
420 "'{}' is not a directory".format(string))
|
|
|
421 if os.access(string, os.R_OK):
|
|
|
422 return string
|
|
|
423 else:
|
|
|
424 raise argparse.ArgumentTypeError(
|
|
|
425 "Directory '{}' is not readable".format(string))
|
|
|
426
|
|
|
427
|
|
|
428 def get_cmdline_args():
|
|
|
429 '''seqclust command line parser'''
|
|
|
430
|
|
|
431 description = """RepeatExplorer:
|
|
|
432 Repetitive sequence discovery and clasification from NGS data
|
|
|
433
|
|
|
434 """
|
|
|
435
|
|
|
436 # arguments parsing
|
|
|
437 parser = argparse.ArgumentParser(description=description,
|
|
|
438 formatter_class=RawTextHelpFormatter)
|
|
|
439 parser.add_argument('-p', '--paired', action='store_true', default=False)
|
|
|
440 parser.add_argument('-A',
|
|
|
441 '--automatic_filtering',
|
|
|
442 action='store_true',
|
|
|
443 default=False)
|
|
|
444 parser.add_argument(
|
|
|
445 '-t',
|
|
|
446 '--tarean_mode',
|
|
|
447 action='store_true',
|
|
|
448 default=False,
|
|
|
449 help="analyze only tandem reapeats without additional classification")
|
|
|
450 parser.add_argument('sequences', type=argparse.FileType('r'))
|
|
|
451 parser.add_argument('-l',
|
|
|
452 '--logfile',
|
|
|
453 type=argparse.FileType('w'),
|
|
|
454 default=None,
|
|
|
455 help='log file, logging goes to stdout if not defines')
|
|
|
456 parser.add_argument('-m',
|
|
|
457 '--mincl',
|
|
|
458 type=float,
|
|
|
459 choices=[Range(0.0, 100.0)],
|
|
|
460 default=0.01)
|
|
|
461 parser.add_argument(
|
|
|
462 '-M',
|
|
|
463 '--merge_threshold',
|
|
|
464 type=float,
|
|
|
465 choices=[0, Range(0.1, 1)],
|
|
|
466 default=0,
|
|
|
467 help=
|
|
|
468 "threshold for mate-pair based cluster merging, default 0 - no merging")
|
|
|
469 parser.add_argument(
|
|
|
470 '-o',
|
|
|
471 '--min_lcov',
|
|
|
472 type=float,
|
|
|
473 choices=[Range(30.0, 80.0)],
|
|
|
474 default=55,
|
|
|
475 help=
|
|
|
476 "minimal overlap coverage - relative to longer sequence length, default 55")
|
|
|
477 parser.add_argument('-c',
|
|
|
478 '--cpu',
|
|
|
479 type=int,
|
|
|
480 default=int(os.environ.get('TAREAN_CPU', 0)),
|
|
|
481 help="number of cpu to use, if 0 use max available")
|
|
|
482 parser.add_argument(
|
|
|
483 '-s',
|
|
|
484 '--sample',
|
|
|
485 type=int,
|
|
|
486 default=0,
|
|
|
487 help="use only sample of input data[by default max reads is used")
|
|
|
488 parser.add_argument(
|
|
|
489 '-P',
|
|
|
490 '--prefix_length',
|
|
|
491 type=int,
|
|
|
492 default=0,
|
|
|
493 help=("If you wish to keep part of the sequences name,\n"
|
|
|
494 " enter the number of characters which should be \n"
|
|
|
495 "kept (1-10) instead of zero. Use this setting if\n"
|
|
|
496 " you are doing comparative analysis"))
|
|
|
497 parser.add_argument('-v',
|
|
|
498 '--output_dir',
|
|
|
499 type=DirectoryType('w'),
|
|
|
500 default="clustering_results")
|
|
|
501 parser.add_argument(
|
|
|
502 '-r',
|
|
|
503 '--max_memory',
|
|
|
504 type=int,
|
|
|
505 default=int(os.environ.get('TAREAN_MAX_MEM', 0)),
|
|
|
506 help=("Maximal amount of available RAM in kB if not set\n"
|
|
|
507 "clustering tries to use whole available RAM"))
|
|
|
508 parser.add_argument(
|
|
|
509 '-d',
|
|
|
510 '--database',
|
|
|
511 default=None,
|
|
|
512 help="fasta file with database for annotation and name of database",
|
|
|
513 nargs=2,
|
|
|
514 action='append')
|
|
|
515
|
|
|
516 parser.add_argument(
|
|
|
517 "-C",
|
|
|
518 "--cleanup",
|
|
|
519 default=False,
|
|
|
520 action="store_true",
|
|
|
521 help="remove unncessary large files from working directory")
|
|
|
522
|
|
|
523 parser.add_argument(
|
|
|
524 "-k",
|
|
|
525 "--keep_names",
|
|
|
526 default=False,
|
|
|
527 action="store_true",
|
|
|
528 help="keep sequence names, by default sequences are renamed")
|
|
|
529
|
|
|
530 parser.add_argument(
|
|
|
531 '-a', '--assembly_min',
|
|
|
532 default=5, type=int,
|
|
|
533 choices=[2,3,4,5],
|
|
|
534 help=('Assembly is performed on individual clusters, by default \n'
|
|
|
535 'clusters with size less then 5 are not assembled. If you \n'
|
|
|
536 'want need assembly of smaller cluster set *assmbly_min* \n'
|
|
|
537 'accordingly\n')
|
|
|
538 )
|
|
|
539
|
|
|
540 parser.add_argument('-tax',
|
|
|
541 '--taxon',
|
|
|
542 default=config.PROTEIN_DATABASE_DEFAULT,
|
|
|
543 choices=list(config.PROTEIN_DATABASE_OPTIONS.keys()),
|
|
|
544 help="Select taxon and protein database version"
|
|
|
545 )
|
|
|
546
|
|
|
547 parser.add_argument(
|
|
|
548 '-opt',
|
|
|
549 '--options',
|
|
|
550 default="ILLUMINA",
|
|
|
551 choices=['ILLUMINA','ILLUMINA_DUST_OFF', 'ILLUMINA_SHORT', 'OXFORD_NANOPORE'])
|
|
|
552
|
|
|
553 parser.add_argument(
|
|
|
554 '-D',
|
|
|
555 '--domain_search',
|
|
|
556 default="BLASTX_W3",
|
|
|
557 choices=['BLASTX_W2', 'BLASTX_W3', 'DIAMOND'],
|
|
|
558 help=
|
|
|
559 ('Detection of protein domains can be performed by either blastx or\n'
|
|
|
560 ' diamond" program. options are:\n'
|
|
|
561 ' BLASTX_W2 - blastx with word size 2 (slowest, the most sesitive)\n'
|
|
|
562 ' BLASTX_W3 - blastx with word size 3 (default)\n'
|
|
|
563 ' DIAMOND - diamond program (significantly faster, less sensitive)\n'
|
|
|
564 'To use this option diamond program must be installed in your PATH'))
|
|
|
565
|
|
|
566 args = parser.parse_args()
|
|
|
567
|
|
|
568 # covert option string to namedtuple of options
|
|
|
569 args.options = getattr(config, args.options)
|
|
|
570 # set protein database
|
|
|
571 args.options = args.options._replace(
|
|
|
572 annotation_search_params=
|
|
|
573 args.options.annotation_search_params._replace(blastx=getattr(
|
|
|
574 config, args.domain_search)))
|
|
|
575 return args
|
|
|
576
|
|
|
577
|
|
|
578 def main():
|
|
|
579 '''
|
|
|
580 Perform graph based clustering
|
|
|
581 '''
|
|
|
582 # argument parsing:
|
|
|
583 args = get_cmdline_args()
|
|
|
584 config.ARGS = args
|
|
|
585 logfile = args.logfile.name if args.logfile else None
|
|
|
586 logging.basicConfig(
|
|
|
587 filename=logfile,
|
|
|
588 format='\n%(asctime)s - %(name)s - %(levelname)s -\n%(message)s\n',
|
|
|
589 level=logging.INFO)
|
|
|
590 config.PROTEIN_DATABASE, config.CLASSIFICATION_HIERARCHY = config.PROTEIN_DATABASE_OPTIONS[
|
|
|
591 args.taxon]
|
|
|
592 # number of CPU to use
|
|
|
593 pipeline_version_info = get_version(config.MAIN_DIR, tarean_mode = args.tarean_mode)
|
|
|
594 config.PROC = args.cpu if args.cpu != 0 else multiprocessing.cpu_count()
|
|
|
595 # TODO add kmer range specification to config - based on the technology
|
|
|
596 r2py.create_connection()
|
|
|
597 try:
|
|
|
598 reporting = r2py.R(config.RSOURCE_reporting, verbose=True)
|
|
|
599 create_annotation = r2py.R(config.RSOURCE_create_annotation,
|
|
|
600 verbose=True)
|
|
|
601 LOGGER.info(args)
|
|
|
602 paths = DataFiles(working_dir=args.output_dir,
|
|
|
603 subdirs=config.DIRECTORY_TREE,
|
|
|
604 files=config.FILES)
|
|
|
605 # files to be included in output
|
|
|
606 for src, dest in config.INCLUDE:
|
|
|
607 shutil.copy(src, os.path.join(paths.working_dir, dest))
|
|
|
608 # geting information about data
|
|
|
609 run_info = DataInfo(args, paths)
|
|
|
610 LOGGER.info(run_info)
|
|
|
611 LOGGER.info(show_object(config))
|
|
|
612 # load all sequences or sample
|
|
|
613 sequences = seqtools.SequenceSet(
|
|
|
614 source=run_info.input_sequences,
|
|
|
615 sample_size=run_info.number_of_reads_for_clustering,
|
|
|
616 paired=run_info.paired,
|
|
|
617 filename=paths.sequences_db,
|
|
|
618 fasta=paths.sequences_fasta,
|
|
|
619 prefix_length=run_info.prefix_length,
|
|
|
620 rename=not run_info.args.keep_names)
|
|
|
621 if run_info.sequence_fiter:
|
|
|
622 n = sequences.remove_sequences_using_filter(
|
|
|
623 run_info.sequence_fiter,
|
|
|
624 keep_proportion=config.FILTER_PROPORTION_OF_KEPT,
|
|
|
625 omitted_sequences_file=paths.filter_omitted,
|
|
|
626 kept_sequences_file=paths.filter_kept
|
|
|
627 )
|
|
|
628 run_info.number_of_omitted_reads = n
|
|
|
629 # add custom databases if provided
|
|
|
630 if args.database:
|
|
|
631 config.CUSTOM_DNA_DATABASE = add_databases(
|
|
|
632 args.database,
|
|
|
633 custom_databases_dir=paths.custom_databases)
|
|
|
634 sequences.makeblastdb(legacy=args.options.legacy_database, lastdb=args.options.lastdb)
|
|
|
635 LOGGER.info("chunksize: {}".format(config.CHUNK_SIZE))
|
|
|
636 sequences.make_chunks(chunk_size=config.CHUNK_SIZE)
|
|
|
637 sequences.create_hitsort(output=paths.hitsort, options=args.options)
|
|
|
638 hitsort = graphtools.Graph(filename=paths.hitsort_db,
|
|
|
639 source=paths.hitsort,
|
|
|
640 paired=run_info.paired,
|
|
|
641 seqids=sequences.keys())
|
|
|
642
|
|
|
643 LOGGER.info('hitsort with {} reads and {} edges loaded.'.format(
|
|
|
644 hitsort.vcount, hitsort.ecount))
|
|
|
645
|
|
|
646 hitsort.save_indexed_graph()
|
|
|
647 LOGGER.info('hitsort index created.')
|
|
|
648
|
|
|
649 hitsort.louvain_clustering(merge_threshold=args.merge_threshold,
|
|
|
650 cleanup=args.cleanup)
|
|
|
651 hitsort.export_cls(path=paths.cls_file)
|
|
|
652 hitsort.adjust_cluster_size(config.FILTER_PROPORTION_OF_KEPT,
|
|
|
653 sequences.ids_kept)
|
|
|
654 sequences.annotate(config.DNA_DATABASE,
|
|
|
655 annotation_name="dna_database",
|
|
|
656 directory=paths.blastn,
|
|
|
657 params=args.options.annotation_search_params.blastn)
|
|
|
658
|
|
|
659 if config.CUSTOM_DNA_DATABASE:
|
|
|
660 LOGGER.info('annotating with custom database')
|
|
|
661 for db, db_name in config.CUSTOM_DNA_DATABASE:
|
|
|
662 sequences.annotate(
|
|
|
663 db,
|
|
|
664 annotation_name=db_name,
|
|
|
665 directory=paths.blastn,
|
|
|
666 params=args.options.annotation_search_params.blastn)
|
|
|
667
|
|
|
668 if not args.tarean_mode:
|
|
|
669 # additional analyses - full RE run
|
|
|
670 # this must be finished befor creating clusters_info
|
|
|
671 sequences.annotate(
|
|
|
672 config.PROTEIN_DATABASE,
|
|
|
673 annotation_name="protein_database",
|
|
|
674 directory=paths.blastx,
|
|
|
675 params=args.options.annotation_search_params.blastx)
|
|
|
676
|
|
|
677 ## annotating using customa databasesreplace
|
|
|
678 LOGGER.info('creating cluster graphs')
|
|
|
679 clusters_info = hitsort.export_clusters_files_multiple(
|
|
|
680 min_size=run_info.mincln,
|
|
|
681 directory=paths.clusters,
|
|
|
682 sequences=sequences,
|
|
|
683 tRNA_database_path=config.TRNA_DATABASE,
|
|
|
684 satellite_model_path=config.SATELLITE_MODEL)
|
|
|
685 if not args.tarean_mode:
|
|
|
686 LOGGER.info("assembling..")
|
|
|
687 assembly_tools.assembly(sequences,
|
|
|
688 hitsort,
|
|
|
689 clusters_info,
|
|
|
690 assembly_dir=paths.assembly,
|
|
|
691 contigs_file=paths.contigs,
|
|
|
692 min_size_of_cluster_for_assembly=args.assembly_min)
|
|
|
693
|
|
|
694 LOGGER.info("detecting LTR in assembly..")
|
|
|
695 for i in clusters_info:
|
|
|
696 i.detect_ltr(config.TRNA_DATABASE)
|
|
|
697
|
|
|
698 run_info.max_annotated_clusters = max([i.index for i in clusters_info])
|
|
|
699 run_info.max_annotated_superclusters = max([i.supercluster
|
|
|
700 for i in clusters_info])
|
|
|
701 # make reports
|
|
|
702 cluster_listing = [i.listing() for i in clusters_info]
|
|
|
703 # make path relative to paths.cluster_info
|
|
|
704 utils.save_as_table(cluster_listing, paths.clusters_info)
|
|
|
705 # creates table cluster_info in hitsort database
|
|
|
706 graphtools.Cluster.add_cluster_table_to_database(cluster_listing,
|
|
|
707 paths.hitsort_db)
|
|
|
708 # export files for consensus sequences, one for each ranks
|
|
|
709 consensus_files = []
|
|
|
710 for i in config.TANDEM_RANKS:
|
|
|
711 consensus_files.append(utils.export_tandem_consensus(
|
|
|
712 clusters_info,
|
|
|
713 path=paths.TR_consensus_fasta.format(i),
|
|
|
714 rank=i))
|
|
|
715
|
|
|
716 if not args.tarean_mode:
|
|
|
717 LOGGER.info("Creating report for superclusters")
|
|
|
718 create_annotation.create_all_superclusters_report(
|
|
|
719 max_supercluster=run_info.max_annotated_superclusters,
|
|
|
720 paths=paths.as_list(),
|
|
|
721 libdir=paths.libdir,
|
|
|
722 superclusters_dir=paths.superclusters,
|
|
|
723 seqdb=paths.sequences_db,
|
|
|
724 hitsortdb=paths.hitsort_db,
|
|
|
725 classification_hierarchy_file=config.CLASSIFICATION_HIERARCHY,
|
|
|
726 HTML_LINKS=dict2lists(config.HTML_LINKS))
|
|
|
727
|
|
|
728 LOGGER.info("Creating report for individual clusters")
|
|
|
729 for cluster in clusters_info:
|
|
|
730 create_annotation.create_cluster_report(
|
|
|
731 cluster.index,
|
|
|
732 seqdb=paths.sequences_db,
|
|
|
733 hitsortdb=paths.hitsort_db,
|
|
|
734 classification_hierarchy_file=
|
|
|
735 config.CLASSIFICATION_HIERARCHY,
|
|
|
736 HTML_LINKS=dict2lists(config.HTML_LINKS))
|
|
|
737
|
|
|
738 LOGGER.info("Creating main html report")
|
|
|
739 reporting.create_main_reports(
|
|
|
740 paths=paths.as_list(),
|
|
|
741 N_clustering=run_info.number_of_reads_for_clustering,
|
|
|
742 N_input=run_info.number_of_input_sequences,
|
|
|
743 N_omit=run_info.number_of_omitted_reads,
|
|
|
744 merge_threshold=args.merge_threshold,
|
|
|
745 paired=run_info.paired,
|
|
|
746 consensus_files=consensus_files,
|
|
|
747 custom_db=bool(config.CUSTOM_DNA_DATABASE),
|
|
|
748 tarean_mode=args.tarean_mode,
|
|
|
749 HTML_LINKS=dict2lists(config.HTML_LINKS),
|
|
|
750 pipeline_version_info=pipeline_version_info,
|
|
|
751 max_memory=run_info.max_memory,
|
|
|
752 max_number_reads_for_clustering=run_info.max_number_reads_for_clustering,
|
|
|
753 mincln=run_info.mincln
|
|
|
754 )
|
|
|
755
|
|
|
756 LOGGER.info("Html report reports created")
|
|
|
757
|
|
|
758 except:
|
|
|
759 r2py.shutdown(config.RSERVE_PORT)
|
|
|
760 raise
|
|
|
761 finally:
|
|
|
762 if args.cleanup:
|
|
|
763 paths.cleanup(config.FILES_TO_DISCARD_AT_CLEANUP)
|
|
|
764 else:
|
|
|
765 LOGGER.info("copy databases to working directory")
|
|
|
766 shutil.copy(paths.sequences_db, paths.working_dir)
|
|
|
767 shutil.copy(paths.hitsort_db, paths.working_dir)
|
|
|
768 # copy log file inside working directory
|
|
|
769 if logfile:
|
|
|
770 shutil.copyfile(logfile, paths.logfile)
|
|
|
771
|
|
|
772
|
|
|
773 if __name__ == "__main__":
|
|
|
774 main()
|
|
|
775 # some error handling here:
|
