Mercurial > repos > fubar > jbrowse2dev
comparison jbrowse2/jbrowse2_jbrowsenotjson.py @ 6:88b9b105c09b draft
Uploaded
author | fubar |
---|---|
date | Fri, 05 Jan 2024 01:58:02 +0000 |
parents | |
children |
comparison
equal
deleted
inserted
replaced
5:42ca8804cd93 | 6:88b9b105c09b |
---|---|
1 #!/usr/bin/env python | |
2 # change to accumulating all configuration for config.json based on the default from the clone | |
3 import argparse | |
4 import datetime | |
5 import hashlib | |
6 import json | |
7 import logging | |
8 import os | |
9 import shutil | |
10 import subprocess | |
11 import tempfile | |
12 import xml.etree.ElementTree as ET | |
13 from collections import defaultdict | |
14 | |
15 logging.basicConfig(level=logging.INFO) | |
16 log = logging.getLogger("jbrowse") | |
17 TODAY = datetime.datetime.now().strftime("%Y-%m-%d") | |
18 GALAXY_INFRASTRUCTURE_URL = None | |
19 mapped_chars = { | |
20 ">": "__gt__", | |
21 "<": "__lt__", | |
22 "'": "__sq__", | |
23 '"': "__dq__", | |
24 "[": "__ob__", | |
25 "]": "__cb__", | |
26 "{": "__oc__", | |
27 "}": "__cc__", | |
28 "@": "__at__", | |
29 "#": "__pd__", | |
30 "": "__cn__", | |
31 } | |
32 | |
33 | |
34 def etree_to_dict(t): | |
35 if t is None: | |
36 return {} | |
37 | |
38 d = {t.tag: {} if t.attrib else None} | |
39 children = list(t) | |
40 if children: | |
41 dd = defaultdict(list) | |
42 for dc in map(etree_to_dict, children): | |
43 for k, v in dc.items(): | |
44 dd[k].append(v) | |
45 d = {t.tag: {k: v[0] if len(v) == 1 else v for k, v in dd.items()}} | |
46 if t.attrib: | |
47 d[t.tag].update(("@" + k, v) for k, v in t.attrib.items()) | |
48 if t.text: | |
49 text = t.text.strip() | |
50 if children or t.attrib: | |
51 if text: | |
52 d[t.tag]["#text"] = text | |
53 else: | |
54 d[t.tag] = text | |
55 return d | |
56 | |
57 | |
58 INSTALLED_TO = os.path.dirname(os.path.realpath(__file__)) | |
59 | |
60 | |
61 def metadata_from_node(node): | |
62 metadata = {} | |
63 try: | |
64 if len(node.findall("dataset")) != 1: | |
65 # exit early | |
66 return metadata | |
67 except Exception: | |
68 return {} | |
69 | |
70 for (key, value) in node.findall("dataset")[0].attrib.items(): | |
71 metadata["dataset_%s" % key] = value | |
72 | |
73 for (key, value) in node.findall("history")[0].attrib.items(): | |
74 metadata["history_%s" % key] = value | |
75 | |
76 for (key, value) in node.findall("metadata")[0].attrib.items(): | |
77 metadata["metadata_%s" % key] = value | |
78 | |
79 for (key, value) in node.findall("tool")[0].attrib.items(): | |
80 metadata["tool_%s" % key] = value | |
81 | |
82 # Additional Mappings applied: | |
83 metadata[ | |
84 "dataset_edam_format" | |
85 ] = '<a target="_blank" href="http://edamontology.org/{0}">{1}</a>'.format( | |
86 metadata["dataset_edam_format"], metadata["dataset_file_ext"] | |
87 ) | |
88 metadata["history_user_email"] = '<a href="mailto:{0}">{0}</a>'.format( | |
89 metadata["history_user_email"] | |
90 ) | |
91 metadata["hist_name"] = metadata["history_display_name"] | |
92 metadata[ | |
93 "history_display_name" | |
94 ] = '<a target="_blank" href="{galaxy}/history/view/{encoded_hist_id}">{hist_name}</a>'.format( | |
95 galaxy=GALAXY_INFRASTRUCTURE_URL, | |
96 encoded_hist_id=metadata["history_id"], | |
97 hist_name=metadata["history_display_name"], | |
98 ) | |
99 metadata[ | |
100 "tool_tool" | |
101 ] = '<a target="_blank" href="{galaxy}/datasets/{encoded_id}/show_params">{tool_id}</a>'.format( | |
102 galaxy=GALAXY_INFRASTRUCTURE_URL, | |
103 encoded_id=metadata["dataset_id"], | |
104 tool_id=metadata["tool_tool_id"], | |
105 # tool_version=metadata['tool_tool_version'], | |
106 ) | |
107 return metadata | |
108 | |
109 | |
110 class JbrowseConnector(object): | |
111 def __init__(self, jbrowse, outdir, genomes, standalone=None): | |
112 self.debug = False | |
113 self.usejson = True | |
114 self.giURL = GALAXY_INFRASTRUCTURE_URL | |
115 self.jbrowse = jbrowse | |
116 self.outdir = outdir | |
117 os.makedirs(self.outdir, exist_ok=True) | |
118 self.genome_paths = genomes | |
119 self.standalone = standalone | |
120 self.trackIdlist = [] | |
121 self.tracksToAdd = [] | |
122 self.config_json = {} | |
123 self.config_json_file = os.path.realpath(os.path.join(outdir, "config.json")) | |
124 if standalone == "complete": | |
125 self.clone_jbrowse(self.jbrowse, self.outdir) | |
126 elif standalone == "minimal": | |
127 self.clone_jbrowse(self.jbrowse, self.outdir, minimal=True) | |
128 | |
129 def subprocess_check_call(self, command, output=None): | |
130 if output: | |
131 if self.debug: | |
132 log.debug("cd %s && %s > %s", self.outdir, " ".join(command), output) | |
133 subprocess.check_call(command, cwd=self.outdir, stdout=output) | |
134 else: | |
135 log.debug("cd %s && %s", self.outdir, " ".join(command)) | |
136 subprocess.check_call(command, cwd=self.outdir) | |
137 | |
138 def subprocess_popen(self, command): | |
139 if self.debug: | |
140 log.debug("cd %s && %s", self.outdir, command) | |
141 p = subprocess.Popen( | |
142 command, | |
143 shell=True, | |
144 stdin=subprocess.PIPE, | |
145 stdout=subprocess.PIPE, | |
146 stderr=subprocess.PIPE, | |
147 ) | |
148 output, err = p.communicate() | |
149 retcode = p.returncode | |
150 if retcode != 0: | |
151 log.error("cd %s && %s", self.outdir, command) | |
152 log.error(output) | |
153 log.error(err) | |
154 raise RuntimeError("Command failed with exit code %s" % (retcode)) | |
155 | |
156 def subprocess_check_output(self, command): | |
157 if self.debug: | |
158 log.debug("cd %s && %s", self.outdir, " ".join(command)) | |
159 return subprocess.check_output(command, cwd=self.outdir) | |
160 | |
161 def _jbrowse_bin(self, command): | |
162 return os.path.realpath(os.path.join(self.jbrowse, "bin", command)) | |
163 | |
164 def symlink_or_copy(self, src, dest): | |
165 if "GALAXY_JBROWSE_SYMLINKS" in os.environ and bool( | |
166 os.environ["GALAXY_JBROWSE_SYMLINKS"] | |
167 ): | |
168 cmd = ["ln", "-s", src, dest] | |
169 else: | |
170 cmd = ["cp", src, dest] | |
171 | |
172 return self.subprocess_check_call(cmd) | |
173 | |
174 def process_genomes(self): | |
175 assemblies = [] | |
176 for i, genome_node in enumerate(self.genome_paths): | |
177 log.info("genome_node=%s" % str(genome_node)) | |
178 # We only expect one input genome per run. This for loop is just | |
179 # easier to write than the alternative / catches any possible | |
180 # issues. | |
181 genome_name = genome_node["meta"]["dataset_dname"] | |
182 dsId = genome_node["meta"]["dataset_id"] | |
183 fapath = genome_node["path"] | |
184 faname = genome_name + ".fa.gz" | |
185 faind = os.path.realpath(os.path.join(self.outdir, faname + ".gzi")) | |
186 if self.standalone == "complete": | |
187 fadest = os.path.realpath(os.path.join(self.outdir, faname)) | |
188 cmd = "bgzip -i -c %s > %s && samtools faidx %s" % ( | |
189 fapath, | |
190 fadest, | |
191 fadest, | |
192 ) | |
193 self.subprocess_popen(cmd) | |
194 adapter = { | |
195 "type": "BgzipFastaAdapter", | |
196 "fastaLocation": { | |
197 "uri": faname, | |
198 }, | |
199 "faiLocation": { | |
200 "uri": faname + ".fai", | |
201 }, | |
202 "gziLocation": { | |
203 "uri": faname + ".gzi", | |
204 }, | |
205 } | |
206 else: | |
207 faurl = "%s/api/datasets/%s/display" % (self.giURL, dsId) | |
208 fastalocation = { | |
209 "uri": faurl, | |
210 } | |
211 failocation = { | |
212 "uri": faname + ".fai", | |
213 } | |
214 adapter = { | |
215 "type": "IndexedFastaAdapter", | |
216 "fastaLocation": fastalocation, | |
217 "faiLocation": failocation, | |
218 } | |
219 | |
220 cmd = ["samtools", "faidx", fapath, "--fai-idx", faind] | |
221 self.subprocess_check_call(cmd) | |
222 trackDict = { | |
223 "name": genome_name, | |
224 "sequence": { | |
225 "type": "ReferenceSequenceTrack", | |
226 "trackId": genome_name, | |
227 "adapter": adapter, | |
228 }, | |
229 "rendering": {"type": "DivSequenceRenderer"}, | |
230 } | |
231 assemblies.append(trackDict) | |
232 self.genome_name = genome_name | |
233 if self.usejson: | |
234 self.config_json["assemblies"] = assemblies | |
235 else: | |
236 cmd = [ | |
237 "jbrowse", | |
238 "add-assembly", | |
239 faname, | |
240 "-t", | |
241 "bgzipFasta", | |
242 "-n", | |
243 genome_name, | |
244 "--load", | |
245 "inPlace", | |
246 "--faiLocation", | |
247 faname + ".fai", | |
248 "--gziLocation", | |
249 faname + ".gzi", | |
250 "--target", | |
251 self.outdir, | |
252 ] | |
253 self.subprocess_check_call(cmd) | |
254 | |
255 def add_default_view(self): | |
256 cmd = [ | |
257 "jbrowse", | |
258 "set-default-session", | |
259 "-s", | |
260 self.config_json_file, | |
261 "-t", | |
262 ",".join(self.trackIdlist), | |
263 "-n", | |
264 "JBrowse2 in Galaxy", | |
265 "--target", | |
266 self.config_json_file, | |
267 "-v", | |
268 " LinearGenomeView", | |
269 ] | |
270 if True or self.debug: | |
271 log.info("### calling set-default-session with cmd=%s" % " ".join(cmd)) | |
272 self.subprocess_check_call(cmd) | |
273 | |
274 def write_config(self): | |
275 with open(self.config_json_file, "w") as fp: | |
276 json.dump(self.config_json, fp) | |
277 | |
278 def add_hic(self, data, trackData): | |
279 """ | |
280 HiC adapter. | |
281 https://github.com/aidenlab/hic-format/blob/master/HiCFormatV9.md | |
282 for testing locally, these work: | |
283 HiC data is from https://s3.amazonaws.com/igv.broadinstitute.org/data/hic/intra_nofrag_30.hic | |
284 using hg19 reference track as a | |
285 'BgzipFastaAdapter' | |
286 fastaLocation: | |
287 uri: 'https://s3.amazonaws.com/jbrowse.org/genomes/GRCh38/fasta/GRCh38.fa.gz', | |
288 faiLocation: | |
289 uri: 'https://s3.amazonaws.com/jbrowse.org/genomes/GRCh38/fasta/GRCh38.fa.gz.fai', | |
290 gziLocation: | |
291 uri: 'https://s3.amazonaws.com/jbrowse.org/genomes/GRCh38/fasta/GRCh38.fa.gz.gzi', | |
292 Cool will not be likely to be a good fit - see discussion at https://github.com/GMOD/jbrowse-components/issues/2438 | |
293 """ | |
294 log.info("#### trackData=%s" % trackData) | |
295 tId = trackData["label"] | |
296 dsId = trackData["metadata"]["dataset_id"] | |
297 url = "%s/api/datasets/%s/display?to_ext=hic " % ( | |
298 self.giURL, | |
299 dsId, | |
300 ) | |
301 hname = trackData["name"] | |
302 if self.standalone == "complete": | |
303 dest = os.path.realpath(os.path.join(self.outdir, hname)) | |
304 url = hname | |
305 cmd = ["cp", data, dest] | |
306 self.subprocess_check_call(cmd) | |
307 floc = { | |
308 "uri": hname, | |
309 } | |
310 else: | |
311 url = "%s/api/datasets/%s/display?to_ext=hic" % (self.giURL, dsId) | |
312 floc = { | |
313 "uri": url, | |
314 } | |
315 trackDict = { | |
316 "type": "HicTrack", | |
317 "trackId": tId, | |
318 "name": hname, | |
319 "assemblyNames": [self.genome_name], | |
320 "adapter": { | |
321 "type": "HicAdapter", | |
322 "hicLocation": floc, | |
323 }, | |
324 } | |
325 if self.usejson: | |
326 self.tracksToAdd.append(trackDict) | |
327 self.trackIdlist.append(tId) | |
328 else: | |
329 cmd = [ | |
330 "jbrowse", | |
331 "add-track", | |
332 url, | |
333 "-t", | |
334 "HicTrack", | |
335 "-a", | |
336 self.genome_name, | |
337 "-n", | |
338 hname, | |
339 "--load", | |
340 "inPlace", | |
341 "--target", | |
342 self.outdir, | |
343 ] | |
344 self.subprocess_check_call(cmd) | |
345 | |
346 def add_maf(self, data, trackData): | |
347 """ | |
348 from https://github.com/cmdcolin/maf2bed | |
349 Note: Both formats start with a MAF as input, and note that your MAF file should contain the species name and chromosome name | |
350 e.g. hg38.chr1 in the sequence identifiers. | |
351 need the reference id - eg hg18, for maf2bed.pl as the first parameter | |
352 """ | |
353 mafPlugin = { | |
354 "plugins": [ | |
355 { | |
356 "name": "MafViewer", | |
357 "url": "https://unpkg.com/jbrowse-plugin-mafviewer/dist/jbrowse-plugin-mafviewer.umd.production.min.js", | |
358 } | |
359 ] | |
360 } | |
361 tId = trackData["label"] | |
362 fname = "%s.bed" % tId | |
363 dest = os.path.realpath("%s/%s" % (self.outdir, fname)) | |
364 # self.symlink_or_copy(data, dest) | |
365 # Process MAF to bed-like. Need build to munge chromosomes | |
366 gname = self.genome_name | |
367 cmd = [ | |
368 "bash", | |
369 os.path.join(INSTALLED_TO, "convertMAF.sh"), | |
370 data, | |
371 gname, | |
372 INSTALLED_TO, | |
373 dest, | |
374 ] | |
375 self.subprocess_check_call(cmd) | |
376 if True or self.debug: | |
377 log.info("### convertMAF.sh called as %s" % " ".join(cmd)) | |
378 # Construct samples list | |
379 # We could get this from galaxy metadata, not sure how easily. | |
380 ps = subprocess.Popen(["grep", "^s [^ ]*", "-o", data], stdout=subprocess.PIPE) | |
381 output = subprocess.check_output(("sort", "-u"), stdin=ps.stdout) | |
382 ps.wait() | |
383 outp = output.decode("ascii") | |
384 soutp = outp.split("\n") | |
385 samp = [x.split("s ")[1] for x in soutp if x.startswith("s ")] | |
386 samples = [x.split(".")[0] for x in samp] | |
387 if self.debug: | |
388 log.info("### got samples = %s " % (samples)) | |
389 trackDict = { | |
390 "type": "MafTrack", | |
391 "trackId": tId, | |
392 "name": trackData["name"], | |
393 "adapter": { | |
394 "type": "MafTabixAdapter", | |
395 "samples": samples, | |
396 "bedGzLocation": { | |
397 "uri": fname + ".sorted.bed.gz", | |
398 }, | |
399 "index": { | |
400 "location": { | |
401 "uri": fname + ".sorted.bed.gz.tbi", | |
402 }, | |
403 }, | |
404 }, | |
405 "assemblyNames": [self.genome_name], | |
406 } | |
407 self.tracksToAdd.append(trackDict) | |
408 self.trackIdlist.append(tId) | |
409 if self.config_json.get("plugins", None): | |
410 self.config_json["plugins"].append(mafPlugin[0]) | |
411 else: | |
412 self.config_json.update(mafPlugin) | |
413 | |
414 def _blastxml_to_gff3(self, xml, min_gap=10): | |
415 gff3_unrebased = tempfile.NamedTemporaryFile(delete=False) | |
416 cmd = [ | |
417 "python", | |
418 os.path.join(INSTALLED_TO, "blastxml_to_gapped_gff3.py"), | |
419 "--trim", | |
420 "--trim_end", | |
421 "--include_seq", | |
422 "--min_gap", | |
423 str(min_gap), | |
424 xml, | |
425 ] | |
426 subprocess.check_call(cmd, cwd=self.outdir, stdout=gff3_unrebased) | |
427 gff3_unrebased.close() | |
428 return gff3_unrebased.name | |
429 | |
430 def add_blastxml(self, data, trackData, blastOpts, **kwargs): | |
431 gff3 = self._blastxml_to_gff3(data, min_gap=blastOpts["min_gap"]) | |
432 | |
433 if "parent" in blastOpts and blastOpts["parent"] != "None": | |
434 gff3_rebased = tempfile.NamedTemporaryFile(delete=False) | |
435 cmd = ["python", os.path.join(INSTALLED_TO, "gff3_rebase.py")] | |
436 if blastOpts.get("protein", "false") == "true": | |
437 cmd.append("--protein2dna") | |
438 cmd.extend([os.path.realpath(blastOpts["parent"]), gff3]) | |
439 subprocess.check_call(cmd, cwd=self.outdir, stdout=gff3_rebased) | |
440 gff3_rebased.close() | |
441 | |
442 # Replace original gff3 file | |
443 shutil.copy(gff3_rebased.name, gff3) | |
444 os.unlink(gff3_rebased.name) | |
445 url = "%s.gff3" % trackData["label"] | |
446 dest = os.path.realpath("%s/%s" % (self.outdir, url)) | |
447 self._sort_gff(gff3, dest) | |
448 url = url + ".gz" | |
449 tId = trackData["label"] | |
450 trackDict = { | |
451 "type": "FeatureTrack", | |
452 "trackId": tId, | |
453 "name": trackData["name"], | |
454 "assemblyNames": [self.genome_name], | |
455 "adapter": { | |
456 "type": "Gff3TabixAdapter", | |
457 "gffGzLocation": { | |
458 "uri": url, | |
459 }, | |
460 "index": { | |
461 "location": { | |
462 "uri": url + ".tbi", | |
463 } | |
464 }, | |
465 }, | |
466 "displays": [ | |
467 { | |
468 "type": "LinearBasicDisplay", | |
469 "displayId": "%s-LinearBasicDisplay" % tId, | |
470 }, | |
471 {"type": "LinearArcDisplay", "displayId": "%s-LinearArcDisplay" % tId}, | |
472 ], | |
473 } | |
474 if self.usejson: | |
475 self.tracksToAdd.append(trackDict) | |
476 self.trackIdlist.append(tId) | |
477 else: | |
478 cmd = [ | |
479 "jbrowse", | |
480 "add-track", | |
481 url, | |
482 "-t", | |
483 "FeatureTrack", | |
484 "-a", | |
485 self.genome_name, | |
486 "--indexFile", | |
487 url + ".tbi", | |
488 "-n", | |
489 trackData["name"], | |
490 "--load", | |
491 "inPlace", | |
492 "--target", | |
493 self.outdir, | |
494 ] | |
495 self.subprocess_check_call(cmd) | |
496 os.unlink(gff3) | |
497 | |
498 def add_bigwig(self, data, trackData): | |
499 url = "%s.bw" % trackData["name"] | |
500 if self.standalone == "complete": | |
501 dest = os.path.realpath(os.path.join(self.outdir, url)) | |
502 cmd = ["cp", data, dest] | |
503 self.subprocess_check_call(cmd) | |
504 bwloc = {"uri": url} | |
505 else: | |
506 dsId = trackData["metadata"]["dataset_id"] | |
507 url = "%s/api/datasets/%s/display?to_ext=fasta" % (self.giURL, dsId) | |
508 bwloc = {"uri": url} | |
509 tId = trackData["label"] | |
510 trackDict = { | |
511 "type": "QuantitativeTrack", | |
512 "trackId": tId, | |
513 "name": url, | |
514 "assemblyNames": [ | |
515 self.genome_name, | |
516 ], | |
517 "adapter": { | |
518 "type": "BigWigAdapter", | |
519 "bigWigLocation": bwloc, | |
520 }, | |
521 "displays": [ | |
522 { | |
523 "type": "LinearWiggleDisplay", | |
524 "displayId": "%s-LinearWiggleDisplay" % tId, | |
525 } | |
526 ], | |
527 } | |
528 if self.usejson: | |
529 self.tracksToAdd.append(trackDict) | |
530 self.trackIdlist.append(tId) | |
531 else: | |
532 cmd = [ | |
533 "jbrowse", | |
534 "add-track", | |
535 url, | |
536 "-t", | |
537 "QuantitativeTrack", | |
538 "-a", | |
539 self.genome_name, | |
540 "-n", | |
541 trackData["name"], | |
542 "--load", | |
543 "inPlace", | |
544 "--target", | |
545 self.outdir, | |
546 ] | |
547 self.subprocess_check_call(cmd) | |
548 | |
549 def add_bam(self, data, trackData, bamOpts, bam_index=None, **kwargs): | |
550 tId = trackData["label"] | |
551 fname = "%s.bam" % trackData["label"] | |
552 dest = os.path.realpath("%s/%s" % (self.outdir, fname)) | |
553 if self.standalone == "complete": | |
554 url = fname | |
555 self.subprocess_check_call(["cp", data, dest]) | |
556 log.info("### copied %s to %s" % (data, dest)) | |
557 bloc = {"uri": url} | |
558 else: | |
559 dsId = trackData["metadata"]["dataset_id"] | |
560 url = "%s/api/datasets/%s/display?to_ext=bam" % (self.giURL, dsId) | |
561 bloc = {"uri": url} | |
562 if bam_index is not None and os.path.exists(os.path.realpath(bam_index)): | |
563 # bai most probably made by galaxy and stored in galaxy dirs, need to copy it to dest | |
564 self.subprocess_check_call( | |
565 ["cp", os.path.realpath(bam_index), dest + ".bai"] | |
566 ) | |
567 else: | |
568 # Can happen in exotic condition | |
569 # e.g. if bam imported as symlink with datatype=unsorted.bam, then datatype changed to bam | |
570 # => no index generated by galaxy, but there might be one next to the symlink target | |
571 # this trick allows to skip the bam sorting made by galaxy if already done outside | |
572 if os.path.exists(os.path.realpath(data) + ".bai"): | |
573 self.symlink_or_copy(os.path.realpath(data) + ".bai", dest + ".bai") | |
574 else: | |
575 log.warn("Could not find a bam index (.bai file) for %s", data) | |
576 trackDict = { | |
577 "type": "AlignmentsTrack", | |
578 "trackId": tId, | |
579 "name": trackData["name"], | |
580 "assemblyNames": [self.genome_name], | |
581 "adapter": { | |
582 "type": "BamAdapter", | |
583 "bamLocation": bloc, | |
584 "index": { | |
585 "location": { | |
586 "uri": fname + ".bai", | |
587 } | |
588 }, | |
589 }, | |
590 } | |
591 if self.usejson: | |
592 self.tracksToAdd.append(trackDict) | |
593 self.trackIdlist.append(tId) | |
594 else: | |
595 cmd = [ | |
596 "jbrowse", | |
597 "add-track", | |
598 fname, | |
599 "-t", | |
600 "AlignmentsTrack", | |
601 "-l", | |
602 "inPlace", | |
603 "-a", | |
604 self.genome_name, | |
605 "--indexFile", | |
606 fname + ".bai", | |
607 "-n", | |
608 trackData["name"], | |
609 "--target", | |
610 self.outdir, | |
611 ] | |
612 self.subprocess_check_call(cmd) | |
613 | |
614 def add_vcf(self, data, trackData): | |
615 tId = trackData["label"] | |
616 url = "%s/api/datasets/%s/display" % ( | |
617 self.giURL, | |
618 trackData["metadata"]["dataset_id"], | |
619 ) | |
620 url = "%s.vcf.gz" % tId | |
621 dest = os.path.realpath("%s/%s" % (self.outdir, url)) | |
622 cmd = "bgzip -c %s > %s" % (data, dest) | |
623 self.subprocess_popen(cmd) | |
624 cmd = ["tabix", "-p", "vcf", dest] | |
625 self.subprocess_check_call(cmd) | |
626 trackDict = { | |
627 "type": "VariantTrack", | |
628 "trackId": tId, | |
629 "name": trackData["name"], | |
630 "assemblyNames": [self.genome_name], | |
631 "adapter": { | |
632 "type": "VcfTabixAdapter", | |
633 "vcfGzLocation": { | |
634 "uri": url, | |
635 }, | |
636 "index": { | |
637 "location": { | |
638 "uri": url + ".tbi", | |
639 } | |
640 }, | |
641 }, | |
642 "displays": [ | |
643 { | |
644 "type": "LinearVariantDisplay", | |
645 "displayId": "%s-LinearVariantDisplay" % tId, | |
646 }, | |
647 { | |
648 "type": "ChordVariantDisplay", | |
649 "displayId": "%s-ChordVariantDisplay" % tId, | |
650 }, | |
651 { | |
652 "type": "LinearPairedArcDisplay", | |
653 "displayId": "%s-LinearPairedArcDisplay" % tId, | |
654 }, | |
655 ], | |
656 } | |
657 if self.usejson: | |
658 self.tracksToAdd.append(trackDict) | |
659 self.trackIdlist.append(tId) | |
660 else: | |
661 cmd = [ | |
662 "jbrowse", | |
663 "add-track", | |
664 url, | |
665 "-t", | |
666 "VariantTrack", | |
667 "-a", | |
668 self.genome_name, | |
669 "--indexFile", | |
670 url + ".tbi", | |
671 "-n", | |
672 trackData["name"], | |
673 "--load", | |
674 "inPlace", | |
675 "--target", | |
676 self.outdir, | |
677 ] | |
678 self.subprocess_check_call(cmd) | |
679 | |
680 def _sort_gff(self, data, dest): | |
681 # Only index if not already done | |
682 if not os.path.exists(dest + ".gz"): | |
683 cmd = "jbrowse sort-gff %s | bgzip -c > %s.gz" % ( | |
684 data, | |
685 dest, | |
686 ) # "gff3sort.pl --precise '%s' | grep -v \"^$\" > '%s'" | |
687 self.subprocess_popen(cmd) | |
688 self.subprocess_check_call(["tabix", "-f", "-p", "gff", dest + ".gz"]) | |
689 | |
690 def _sort_bed(self, data, dest): | |
691 # Only index if not already done | |
692 if not os.path.exists(dest): | |
693 cmd = "sort -k1,1 -k2,2n %s | bgzip -c > %s" % (data, dest) | |
694 self.subprocess_popen(cmd) | |
695 cmd = ["tabix", "-f", "-p", "bed", dest] | |
696 self.subprocess_check_call(cmd) | |
697 | |
698 def add_gff(self, data, ext, trackData): | |
699 url = "%s.%s" % (trackData["label"], ext) | |
700 dest = os.path.realpath("%s/%s" % (self.outdir, url)) | |
701 self._sort_gff(data, dest) | |
702 url = url + ".gz" | |
703 tId = trackData["label"] | |
704 trackDict = { | |
705 "type": "FeatureTrack", | |
706 "trackId": tId, | |
707 "name": trackData["name"], | |
708 "assemblyNames": [self.genome_name], | |
709 "adapter": { | |
710 "type": "Gff3TabixAdapter", | |
711 "gffGzLocation": { | |
712 "uri": url, | |
713 }, | |
714 "index": { | |
715 "location": { | |
716 "uri": url + ".tbi", | |
717 } | |
718 }, | |
719 }, | |
720 "displays": [ | |
721 { | |
722 "type": "LinearBasicDisplay", | |
723 "displayId": "%s-LinearBasicDisplay" % tId, | |
724 }, | |
725 {"type": "LinearArcDisplay", "displayId": "%s-LinearArcDisplay" % tId}, | |
726 ], | |
727 } | |
728 if self.usejson: | |
729 self.tracksToAdd.append(trackDict) | |
730 self.trackIdlist.append(tId) | |
731 else: | |
732 cmd = [ | |
733 "jbrowse", | |
734 "add-track", | |
735 url, | |
736 "-t", | |
737 "FeatureTrack", | |
738 "-a", | |
739 self.genome_name, | |
740 "-n", | |
741 trackData["name"], | |
742 "--load", | |
743 "inPlace", | |
744 "--target", | |
745 self.outdir, | |
746 ] | |
747 self.subprocess_check_call(cmd) | |
748 | |
749 def add_bed(self, data, ext, trackData): | |
750 url = "%s.%s" % (trackData["label"], ext) | |
751 dest = os.path.realpath("%s/%s.gz" % (self.outdir, url)) | |
752 self._sort_bed(data, dest) | |
753 tId = trackData["label"] | |
754 url = url + ".gz" | |
755 trackDict = { | |
756 "type": "FeatureTrack", | |
757 "trackId": tId, | |
758 "name": trackData["name"], | |
759 "assemblyNames": [self.genome_name], | |
760 "adapter": { | |
761 "type": "BedTabixAdapter", | |
762 "bedGzLocation": { | |
763 "uri": url, | |
764 }, | |
765 "index": { | |
766 "location": { | |
767 "uri": url + ".tbi", | |
768 } | |
769 }, | |
770 }, | |
771 "displays": [ | |
772 { | |
773 "type": "LinearBasicDisplay", | |
774 "displayId": "%s-LinearBasicDisplay" % tId, | |
775 }, | |
776 {"type": "LinearArcDisplay", "displayId": "%s-LinearArcDisplay" % tId}, | |
777 ], | |
778 } | |
779 if self.usejson: | |
780 self.tracksToAdd.append(trackDict) | |
781 self.trackIdlist.append(tId) | |
782 else: | |
783 cmd = [ | |
784 "jbrowse", | |
785 "add-track", | |
786 url, | |
787 "-t", | |
788 "FeatureTrack", | |
789 "-a", | |
790 self.genome_name, | |
791 "--indexFile", | |
792 url + ".tbi", | |
793 "-n", | |
794 trackData["name"], | |
795 "--load", | |
796 "inPlace", | |
797 "--target", | |
798 self.outdir, | |
799 ] | |
800 self.subprocess_check_call(cmd) | |
801 | |
802 def process_annotations(self, track): | |
803 category = track["category"].replace("__pd__date__pd__", TODAY) | |
804 for i, ( | |
805 dataset_path, | |
806 dataset_ext, | |
807 track_human_label, | |
808 extra_metadata, | |
809 ) in enumerate(track["trackfiles"]): | |
810 # Unsanitize labels (element_identifiers are always sanitized by Galaxy) | |
811 for key, value in mapped_chars.items(): | |
812 track_human_label = track_human_label.replace(value, key) | |
813 outputTrackConfig = { | |
814 "category": category, | |
815 } | |
816 if self.debug: | |
817 log.info( | |
818 "Processing category = %s, track_human_label = %s", | |
819 category, | |
820 track_human_label, | |
821 ) | |
822 # We add extra data to hash for the case of REST + SPARQL. | |
823 if ( | |
824 "conf" in track | |
825 and "options" in track["conf"] | |
826 and "url" in track["conf"]["options"] | |
827 ): | |
828 rest_url = track["conf"]["options"]["url"] | |
829 else: | |
830 rest_url = "" | |
831 | |
832 # I chose to use track['category'] instead of 'category' here. This | |
833 # is intentional. This way re-running the tool on a different date | |
834 # will not generate different hashes and make comparison of outputs | |
835 # much simpler. | |
836 hashData = [ | |
837 str(dataset_path), | |
838 track_human_label, | |
839 track["category"], | |
840 rest_url, | |
841 ] | |
842 hashData = "|".join(hashData).encode("utf-8") | |
843 outputTrackConfig["label"] = hashlib.md5(hashData).hexdigest() + "_%s" % i | |
844 outputTrackConfig["metadata"] = extra_metadata | |
845 outputTrackConfig["name"] = track_human_label | |
846 | |
847 if dataset_ext in ("gff", "gff3"): | |
848 self.add_gff( | |
849 dataset_path, | |
850 dataset_ext, | |
851 outputTrackConfig, | |
852 ) | |
853 elif dataset_ext in ("hic",): | |
854 self.add_hic( | |
855 dataset_path, | |
856 outputTrackConfig, | |
857 ) | |
858 elif dataset_ext in ("bed",): | |
859 self.add_bed( | |
860 dataset_path, | |
861 dataset_ext, | |
862 outputTrackConfig, | |
863 ) | |
864 elif dataset_ext in ("maf",): | |
865 self.add_maf( | |
866 dataset_path, | |
867 outputTrackConfig, | |
868 ) | |
869 elif dataset_ext == "bigwig": | |
870 self.add_bigwig( | |
871 dataset_path, | |
872 outputTrackConfig, | |
873 ) | |
874 elif dataset_ext == "bam": | |
875 real_indexes = track["conf"]["options"]["pileup"]["bam_indices"][ | |
876 "bam_index" | |
877 ] | |
878 if not isinstance(real_indexes, list): | |
879 # <bam_indices> | |
880 # <bam_index>/path/to/a.bam.bai</bam_index> | |
881 # </bam_indices> | |
882 # | |
883 # The above will result in the 'bam_index' key containing a | |
884 # string. If there are two or more indices, the container | |
885 # becomes a list. Fun! | |
886 real_indexes = [real_indexes] | |
887 | |
888 self.add_bam( | |
889 dataset_path, | |
890 outputTrackConfig, | |
891 track["conf"]["options"]["pileup"], | |
892 bam_index=real_indexes[i], | |
893 ) | |
894 elif dataset_ext == "blastxml": | |
895 self.add_blastxml( | |
896 dataset_path, outputTrackConfig, track["conf"]["options"]["blast"] | |
897 ) | |
898 elif dataset_ext == "vcf": | |
899 self.add_vcf(dataset_path, outputTrackConfig) | |
900 else: | |
901 log.warn("Do not know how to handle %s", dataset_ext) | |
902 | |
903 def clone_jbrowse(self, jbrowse_dir, destination, minimal=False): | |
904 """Clone a JBrowse directory into a destination directory.""" | |
905 cmd = ["jbrowse", "create", "-f", self.outdir] | |
906 self.subprocess_check_call(cmd) | |
907 for fn in [ | |
908 "asset-manifest.json", | |
909 "favicon.ico", | |
910 "robots.txt", | |
911 "umd_plugin.js", | |
912 "version.txt", | |
913 "test_data", | |
914 ]: | |
915 cmd = ["rm", "-rf", os.path.join(self.outdir, fn)] | |
916 self.subprocess_check_call(cmd) | |
917 | |
918 | |
919 if __name__ == "__main__": | |
920 parser = argparse.ArgumentParser(description="", epilog="") | |
921 parser.add_argument("xml", type=argparse.FileType("r"), help="Track Configuration") | |
922 | |
923 parser.add_argument("--jbrowse", help="Folder containing a jbrowse release") | |
924 parser.add_argument("--outdir", help="Output directory", default="out") | |
925 parser.add_argument( | |
926 "--standalone", | |
927 choices=["complete", "minimal", "data"], | |
928 help="Standalone mode includes a copy of JBrowse", | |
929 ) | |
930 parser.add_argument("--version", "-V", action="version", version="%(prog)s 0.8.0") | |
931 args = parser.parse_args() | |
932 | |
933 tree = ET.parse(args.xml.name) | |
934 root = tree.getroot() | |
935 | |
936 # This should be done ASAP | |
937 GALAXY_INFRASTRUCTURE_URL = root.find("metadata/galaxyUrl").text | |
938 # Sometimes this comes as `localhost` without a protocol | |
939 if not GALAXY_INFRASTRUCTURE_URL.startswith("http"): | |
940 # so we'll prepend `http://` and hope for the best. Requests *should* | |
941 # be GET and not POST so it should redirect OK | |
942 GALAXY_INFRASTRUCTURE_URL = "http://" + GALAXY_INFRASTRUCTURE_URL | |
943 | |
944 jc = JbrowseConnector( | |
945 jbrowse=args.jbrowse, | |
946 outdir=args.outdir, | |
947 genomes=[ | |
948 { | |
949 "path": os.path.realpath(x.attrib["path"]), | |
950 "meta": metadata_from_node(x.find("metadata")), | |
951 } | |
952 for x in root.findall("metadata/genomes/genome") | |
953 ], | |
954 standalone=args.standalone, | |
955 ) | |
956 jc.process_genomes() | |
957 | |
958 for track in root.findall("tracks/track"): | |
959 track_conf = {} | |
960 track_conf["trackfiles"] = [] | |
961 | |
962 is_multi_bigwig = False | |
963 try: | |
964 if track.find("options/wiggle/multibigwig") and ( | |
965 track.find("options/wiggle/multibigwig").text == "True" | |
966 ): | |
967 is_multi_bigwig = True | |
968 multi_bigwig_paths = [] | |
969 except KeyError: | |
970 pass | |
971 | |
972 trackfiles = track.findall("files/trackFile") | |
973 if trackfiles: | |
974 for x in track.findall("files/trackFile"): | |
975 if is_multi_bigwig: | |
976 multi_bigwig_paths.append( | |
977 (x.attrib["label"], os.path.realpath(x.attrib["path"])) | |
978 ) | |
979 else: | |
980 if trackfiles: | |
981 metadata = metadata_from_node(x.find("metadata")) | |
982 track_conf["dataset_id"] = metadata["dataset_id"] | |
983 track_conf["trackfiles"].append( | |
984 ( | |
985 os.path.realpath(x.attrib["path"]), | |
986 x.attrib["ext"], | |
987 x.attrib["label"], | |
988 metadata, | |
989 ) | |
990 ) | |
991 else: | |
992 # For tracks without files (rest, sparql) | |
993 track_conf["trackfiles"].append( | |
994 ( | |
995 "", # N/A, no path for rest or sparql | |
996 track.attrib["format"], | |
997 track.find("options/label").text, | |
998 {}, | |
999 ) | |
1000 ) | |
1001 | |
1002 if is_multi_bigwig: | |
1003 metadata = metadata_from_node(x.find("metadata")) | |
1004 | |
1005 track_conf["trackfiles"].append( | |
1006 ( | |
1007 multi_bigwig_paths, # Passing an array of paths to represent as one track | |
1008 "bigwig_multiple", | |
1009 "MultiBigWig", # Giving an hardcoded name for now | |
1010 {}, # No metadata for multiple bigwig | |
1011 ) | |
1012 ) | |
1013 | |
1014 track_conf["category"] = track.attrib["cat"] | |
1015 track_conf["format"] = track.attrib["format"] | |
1016 try: | |
1017 # Only pertains to gff3 + blastxml. TODO? | |
1018 track_conf["style"] = {t.tag: t.text for t in track.find("options/style")} | |
1019 except TypeError: | |
1020 track_conf["style"] = {} | |
1021 pass | |
1022 track_conf["conf"] = etree_to_dict(track.find("options")) | |
1023 jc.process_annotations(track_conf) | |
1024 print("## processed", str(track_conf), "trackIdlist", jc.trackIdlist) | |
1025 print( | |
1026 "###done processing, trackIdlist=", | |
1027 jc.trackIdlist, | |
1028 "config=", | |
1029 str(jc.config_json), | |
1030 ) | |
1031 jc.config_json["tracks"] = jc.tracksToAdd | |
1032 if jc.usejson: | |
1033 jc.write_config() | |
1034 jc.add_default_view() |