#!/usr/bin/python2 """Generate JSON datahub file for data from 5x5 project. Originally written by Daofeng (/bar/dli/pyScript/gen_hub_json_twhmm.py) !! IMPORTANT: RUN WITH PYTHON2. !! Example $ cat methylCRF_prehub.txt #sample specimen tissue assay file TW534 TW630 ZYA23-colon1 methylCRF ZYA23-colon1 colon1 methylCRF /bar/twlab-shared/Epigenome_Evolution/human/batch_methylcrf/methylcrf_practice/ZYA23-colon1_mcrf.bed TW535 TW631 ZYA24-colon2 methylCRF ZYA24-colon2 colon2 methylCRF /bar/twlab-shared/Epigenome_Evolution/human/batch_methylcrf/methylcrf_practice/ZYA24-colon2_mcrf.bed $ python2 gen_hub_json_twhmm.py --file methylcrf_prehub --output-directory /bar/nrockweiler/public_html/hmm/methylcrf --assembly hg19 """ # TODO # Correct hierarchy of metadata terms. Use browser metadata term ids. # Make sure umask is set right so final created files are readable by others # Xiaoyu get it from Nicole. don't really need to touch it, just prepare the pre-hub.txt file. from __future__ import print_function import sys # Check if running python < 3 if sys.version_info >= (3,0): print('ERROR: script requires Python 2.x, not Python 3.x. Rerun with \'python2 gen_hub_json_twhmm.py\'\n\n', file=sys.stderr) sys.exit(1) def main(): args = _process_cmd_line() if not args.dry_run: _mk_dir_safe(args.output_directory) (sample_metadata_lookup_table, assay_lookup_table, institution_lookup_table) = _load_metadata(args.sample_metadata_file, args.experimental_metadata_file, args.institution_metadata_file) dat = [] c = 0 # Determine datahub track information for each sample with open(args.file, 'rU') as input_fh: for line in input_fh: line = line.rstrip('\n') # Skip comment and blank rows if line.startswith('#'): continue elif not line: continue cols = line.split('\t') sample_name = cols[0] cell_line = cols[1] developmental_stage = cols[2] strain = cols[3] tissue = cols[4] assay = cols[5] institution = cols[6] track_fn = cols[7] details_source = cols[8] # print("sample_name", sample_name, "cell_line", cell_line, "developmental_stage", developmental_stage, "strain", strain, "tissue", tissue, "assay", assay, "institution", institution, "track_fn", track_fn, sep="|") # Lookup the metadata ids that correspond to the metadata terms cell_line_id = _metadata2id(cell_line, sample_metadata_lookup_table) developmental_stage_id = _metadata2id(developmental_stage, sample_metadata_lookup_table) strain_id = _metadata2id(strain, sample_metadata_lookup_table) tissue_id = _metadata2id(tissue, sample_metadata_lookup_table) assay_id = _metadata2id(assay, assay_lookup_table) institution_id = _metadata2id(institution, institution_lookup_table) if (assay == 'methylcrf'): assay = 'methylCRF' # Convention used by trackDbParser track_bn = os.path.basename(track_fn) (track_bn_no_ext, track_ext) = os.path.splitext(track_bn) track_ext = track_ext[1:].lower() # Remove leading '.' and lowercase # Create output browser-friendly track files # If the track file is medip/mre, file should be a bigwig (binary), softlink file to output directory # If the track file is methylCRF, convert to bedGraph, bgzip, and tabix index in the output directory # If the track file is mnm, convert to bedgraph, bgzip and tabix index in the output directory if (assay == 'medip') | (assay == 'mre'): (final_track_bn, final_track_ext) = _link_file(args.output_directory, track_fn, track_bn, track_ext, args.dry_run) elif ((assay == 'methylmnm') | (assay == 'methylCRF')): if assay == 'methylmnm': # Create bedgraph track_bedgraph_fn = os.path.join(args.output_directory, track_bn_no_ext + '.bedgraph') if not args.dry_run: _methylmnm2bedgraph(track_fn, track_bedgraph_fn, args.overwrite, args.easy_going) else: # Create bedgraph track_bedgraph_fn = os.path.join(args.output_directory, track_bn_no_ext + '.bedgraph') if not args.dry_run: if _can_create_output_f(track_bedgraph_fn, args.overwrite, args.easy_going): methylcrf2bedgraph_cmd = ['cut', '-f', '1-3,5', track_fn] _run_cmd(methylcrf2bedgraph_cmd, track_bedgraph_fn, args.overwrite) # bzip track_bedgraph_bzip_fn = _create_bzip(track_bedgraph_fn, args.overwrite, args.easy_going, args.dry_run) final_track_bn = os.path.basename(track_bedgraph_bzip_fn) final_track_ext = 'bedgraph' # tabix track_bedgraph_tabix_fn = _create_tabix(track_bedgraph_bzip_fn, args.overwrite, args.easy_going, args.dry_run) # Try to guess from the file extension # bigwig elif ((track_ext == 'bw') or (track_ext == 'bigwig')): track_ext = 'bigwig' # Where to softlink the file? I think it should go here: /srv/epgg/data/data/subtleKnife/dm6 (final_track_bn, final_track_ext) = _link_file(args.output_directory, track_fn, track_bn, track_ext, args.dry_run) else: print('ERROR: unrecognized assay \'%s\'\n' % (assay), file=sys.stderr) sys.exit(1) if final_track_bn is None: continue i = {} i['name'] = sample_name i['url'] = '/'.join([args.url, final_track_bn]) i['type'] = final_track_ext i['metadata'] = {} # print("cell_line_id", cell_line_id, "tissue_id", tissue_id, "strain_id", strain_id, "developmental_stage_id", developmental_stage_id, "assay_id", assay_id, "institution_id", institution_id, sep="\t") i['metadata']['md1'] = [cell_line_id, tissue_id, strain_id, developmental_stage_id] i['metadata']['md2'] = [assay_id] i['metadata']['md3'] = [institution_id] i['details'] = {} i['details']['source'] = details_source if c <= args.max_show_tracks: i['mode'] = 'show' else: i['mode'] = 'hide' if assay == 'medip': if track_ext == 'bigwig': i['qtc'] = {'smooth': 7} col = tkColor[assay] i['colorpositive'] = 'rgb({})'.format(col) i['height'] = 30 elif assay == 'methylCRF': i['barplot_bg'] = '#C0C0C0' i['fixedscale'] = {'min':0,'max':1} col = tkColor[assay] i['colorpositive'] = 'rgb({})'.format(col) i['height'] = 30 elif assay == 'methylmnm': i['qtc'] = {'anglescale':1, 'pr':255, 'pg':102, 'pb':51, 'nr':0, 'ng':119, 'nb':158, 'pth':'rgb(178,71,35)', 'nth':'rgb(0,83,110)', 'thtype':2, 'thmin':0, 'thmax':10, 'thpercentile':100, 'height':19, 'summeth':1 } elif track_ext == 'bigwig': i['qtc'] = {'anglescale':1, 'anglescale': 1, 'pr': 0, 'pg': 0, 'pb': 230, 'nr': 255, 'ng': 0, 'nb': 0, 'pth': '#000099', 'nth': '#800000', 'thtype': 0, 'thmin': 0, 'thmax': 10, 'thpercentile': 90, 'height': 10, 'summeth': 1 } else: col = tkColor[assay] i['colorpositive'] = 'rgb({})'.format(col) i['height'] = 30 dat.append(i) c += 1 # Add native tracks nat = {} nat['type'] = 'native_track' natlis = [] nat['list'] = natlis natlis.append({'name':'refGene','mode':'full'}) dat.append(nat) # Add metadata track information cm = {} cm['type'] = 'metadata' cm['vocabulary_set'] = {} cm['vocabulary_set']['md1'] = args.sample_metadata_url cm['vocabulary_set']['md2'] = args.experimental_metadata_url cm['vocabulary_set']['md3'] = args.institution_metadata_url cm['show_terms'] = {} cm['show_terms']['md1'] = ['Cell Line', 'Developmental Stage', 'Strain', 'Tissue'] cm['show_terms']['md2'] = ['Assay'] dat.append(cm) # Create dathub file if args.dry_run: json.dump(dat, sys.stdout, sort_keys=True) print('INFO: dry run datahub file would have been created here \'%s\' with %d tracks.\n' % (args.output_hub, len(dat) - 2), file=sys.stderr) print('INFO: dry run datahub would have been viewable at \'http://epigenomegateway.wustl.edu/browser/?genome=%s&datahub=%s\'\n' % (args.assembly, args.output_hub_url), file=sys.stderr) else: if _can_create_output_f(args.output_hub, args.overwrite, args.easy_going): with open(args.output_hub, 'w') as hub_fh: json.dump(dat, hub_fh, sort_keys=True) print('INFO: created datahub file \'%s\' with %d tracks.\n' % (args.output_hub, len(dat) - 2), file=sys.stderr) print('INFO: view the datahub at \'http://epigenomegateway.wustl.edu/browser/?genome=%s&datahub=%s\'\n' % (args.assembly, args.output_hub_url), file=sys.stderr) def _process_cmd_line(): """Parse the command-line arguments Returns: args (argparse.Namespace): list of arguments """ # Have your cake and eat it too (http://stackoverflow.com/questions/18462610/argumentparser-epilog-and-description-formatting-in-conjunction-with-argumentdef) class CustomFormatter(argparse.ArgumentDefaultsHelpFormatter, argparse.RawDescriptionHelpFormatter): pass parser = argparse.ArgumentParser( formatter_class=CustomFormatter, description=__doc__, ) parser.add_argument( '--file', '-f', required=True, help='path to file containing list of track files to add to the hub. File should be tab-delimited with the following columns: 1) sample name, 2) sample, 3) tissue, 4) assay, 5) data file. (See note at end about file formats.) Blank rows and rows starting with # are ignored. The order of the tracks will be in the same order as this file. For medip and mre tracks, the data file should be in bigwig format; for methylCRF tracks the file should be a bed of the methylation predictions, e.g., ZYA28-testis2_mcrf.bed; for MnM tracks, the data file should be a tsv of the significant DMRs (i.e., from the MnM.selectDMR step), e.g., DMR_e5_testis2_brain2_human.bed.\n' ) parser.add_argument( '--output-directory', '-o', default=os.path.join(os.path.sep, 'bar', getpass.getuser(), 'public_html'), help='path to write output files. Path should be accessible by the browser, i.e., it must start with /bar//public_html. The directory is created if it doesn\'t already exist.', ) parser.add_argument( '--url', '-u', help='URL of output directory. If not specified, the output URL is assumed to be http://wangftp.wustl.edu/~/.', ) parser.add_argument( '--output-hub', '-oh', default='hub', help='basename of output hub file. The hub file is created in the output directory UNLESS --dry-run/-dr is specfied. In this case, the hub file is printed to STDOUT.', ) parser.add_argument( '--assembly', '-a', default='hg19', choices=['hg19', 'mm9', 'rn4', 'rheMac3', 'canFam3', 'galGal4', 'dm6', 'dm3'], help='genome assembly for tracks. Used for generating the browser link. Option is CASE SENSITIVE.', ) parser.add_argument( '--max-show-tracks', '-m', type=int, default=10, help='Maximum numer of tracks to show.', # Should probably make a way to have no limit, or could just use a large number ) parser.add_argument( '--sample-metadata-url', '-s', default='http://vizhub.wustl.edu/metadata/fruitfly/Fly_Sample_metadata', help='URL for sample metadata.', ) parser.add_argument( '--experimental-metadata-url', '-ex', default='http://vizhub.wustl.edu/metadata/Experimental_assays', help='URL experimental assays for metadata.', ) parser.add_argument( '--institution-metadata-url', '-i', default='http://vizhub.wustl.edu/metadata/Institutions', help='URL for institution metadata.', ) parser.add_argument( '--sample-metadata-file', '-sf', default='/ale/expr/nrockweiler/browser/code/meta/fruitfly/Fly_Sample_metadata', help='path for sample metadata.', ) parser.add_argument( '--experimental-metadata-file', '-exf', default='/ale/expr/nrockweiler/browser/code/meta/Experimental_assays', help='path experimental assays for metadata.', ) parser.add_argument( '--institution-metadata-file', '-if', default='/ale/expr/nrockweiler/browser/code/meta/Institutions', help='path for institution metadata.', ) parser.add_argument( '--dry-run', '-dr', action='store_true', default=False, help='specify option to only print hub file to STDOUT. No other files/links/directories will be created. Useful when you want to create a hub file but don\'t have the permission to write in output directory, e.g., when you\'re creating a hub for someone else..', ) # Only --easy-going and --overwrite are mutually exclusive. group = parser.add_mutually_exclusive_group() group.add_argument( '--easy-going', '-e', action='store_true', default=False, help='specify option to only create output files if they don\'t already exist. Useful for when you want to subset or combine existing datahubs. Use with caution: if only final files exist, but initial files don\'t, only the initial files will be made and may no longer be consistent with the final files.', ) group.add_argument( '--overwrite', '-ow', action='store_true', default=False, help='specify option to write overwrite existing output files.', ) args = parser.parse_args() args.output_directory = os.path.realpath(args.output_directory) # Construct URL if args.url is None: FTP_SITE = 'http://wangftp.wustl.edu' m = re.match('/bar/(\w+)/public_html?/', args.output_directory) if m: output_user = m.group(1) args.url = '/'.join([FTP_SITE, '~' + output_user, args.output_directory[m.end(0):]]) else: print('ERROR: could not find username in public html dir\n\n', file=sys.stderr) sys.exit(1) else: args.url = args.url.rstrip('/') args.output_hub_url = '/'.join([args.url, args.output_hub]) args.output_hub = os.path.join(args.output_directory, args.output_hub) if args.dry_run: print('INFO: dry run output directory (containing hub file and tracks) would have been \'%s\'.\n' % (args.output_directory), file=sys.stderr) else: print('INFO: output directory (containing hub file and tracks) is \'%s\'.\n' % (args.output_directory), file=sys.stderr) return(args) def _mk_dir_safe(d): """Create directory if it doesn't already exist. Avoids race condition when checking and creating a directory. See http://stackoverflow.com/questions/273192/in-python-check-if-a-directory-exists-and-create-it-if-necessary. Args: d (str): directory name. """ try: os.makedirs(d) except OSError as exception: if exception.errno != errno.EEXIST: # if error is not because the dir already existed, raise it raise def _can_create_output_f(output_f, overwrite, easy_going): """Check if the output file can be created. Args: output_f (str): output file. overwrite (bool): True if the output file can be overwritten, false otherwise. easy_going (bool): True if the output file should be created only if it doesn't already exist, false otherwise. Returns: create_output_f (bool): True if the output file should be crated, false otherwise. """ create_output_f = True if os.path.exists(output_f): if not overwrite and not easy_going: create_output_f = False print('ERROR: output file \'%s\' already exists and you didn\'t want to overwrite it or use the existing file. Use --overwrite/-ow to overwrite existing output files or --easy-going/-e to use existing files.\n' % (output_f), file=sys.stderr) sys.exit(1) elif easy_going: create_output_f = False print('INFO: output file \'%s\' already exists and you wanted to keep it. Use --overwrite to overwrite the file.\n' % (output_f), file=sys.stderr) return create_output_f def _run_cmd(cmd, output_f, overwrite=False, mode='w'): """Execute external command. Args: cmd (list): sequence of command arguments output_f (str): name of file to write STDOUT to overwrite (flag): True if the output file should be overwritten, false otherwise. mode (str): 'w' to write to output file, 'a' to append to output file. """ if ((os.path.exists(output_f)) and (mode != 'a')): if not overwrite: print('ERROR: output file \'%s\' already exists and you didn\'t want to overwrite it or use the existing file. Use --overwrite/-ow to overwrite existing output files or --easy-going/-e to use existing files.\n' % (output_f), file=sys.stderr) sys.exit(1) # Note: Can't use 'with' and subprocess.Popen in python < v3 (http://bugs.python.org/issue13202) p = subprocess.Popen(cmd, stdout=open(output_f, mode)) p.wait() return_code = p.returncode if return_code: print('ERROR: command \'%s\' had a nonzero (%d) exit status.\n' % (' '.join(cmd), return_code), file=sys.stderr) sys.exit(1) def _methylmnm2bedgraph(mnm_fn, output_bedgraph_fn, overwrite, easy_going): """Create bedgraph of methylMnM data. MethylMnM .tsv file is converted to bedgraph format. The score in the bedgraph is defined as -log10 q-value * sign(DMR). The header is removed (if present). Args: mnm_fn (str): MethylMnM tsv file containing the significant DMRs (i.e., from the MnM.selectDMR step), e.g., DMR_e5_testis2_brain2_human.bed output_bedgraph_fn (str): output file. overwrite (bool): True if the output file can be overwritten, false otherwise. easy_going (bool): True if the output file should be created only if it doesn't already exist, false otherwise. """ if _can_create_output_f(output_bedgraph_fn, overwrite, easy_going): with open(mnm_fn, 'rU') as mnm_fh, open(output_bedgraph_fn, 'w') as output_bedgraph_fh: for line in mnm_fh: # Skip header if line.startswith('chr\tchrSt'): continue line = line.rstrip('\n') cols = line.split('\t') chrom = cols[0] start = cols[1] end = cols[2] test_statistic = float(cols[10]) q_value = float(cols[11]) if q_value == 0: # TODO can the browser display Inf? # Min q-value (excluding 0) for 5x5 MnM data was 4E-58. Not sure what the lower limit is on the statistical test so decided on a value a little lower than that: 1E-60 q_value = 1E-60 score = -1 * math.log(q_value, 10) * math.copysign(1, test_statistic) # score = -log10(q-value) * sign(DMR) print(chrom, start, end, score, sep="\t", file=output_bedgraph_fh) def _create_bzip(fn, overwrite, easy_going, dry_run): """bzip file. Compression is done in place, i.e., the original file will be removed and replaced with the compressed version. Args: fn (str): file to compress. overwrite (bool): True if the output file can be overwritten, false otherwise. easy_going (bool): True if the output file should be created only if it doesn't already exist, false otherwise. Returns: bgzip_output_f (str): compressed file. """ bgzip_output_f = fn + '.gz' if not dry_run: if _can_create_output_f(bgzip_output_f, overwrite, easy_going): bgzip_cmd = ['bgzip', '-f', fn] # Using -f (force) because already checked that output existed subprocess.check_call(bgzip_cmd) return(bgzip_output_f) def _create_tabix(fn, overwrite, easy_going, dry_run): """Create tabix index of a file. Args: fn (str): file to index. overwrite (bool): True if the output file can be overwritten, false otherwise. easy_going (bool): True if the output file should be created only if it doesn't already exist, false otherwise. Returns: tabix_output_f (str): index file. """ tabix_output_f = fn + '.tbi' if not dry_run: if _can_create_output_f(tabix_output_f, overwrite, easy_going): tabix_cmd = ['tabix', '-p', 'bed', '-f', fn] # Using -f (force) because already checked that output existed subprocess.check_call(tabix_cmd) return(tabix_output_f) def _link_file(output_directory, track_fn, track_bn, track_ext, dry_run): target = os.path.join(output_directory, track_bn) if os.path.exists(track_fn): if not dry_run: # The file must be readable by the world to view it on the browser st = os.stat(track_fn) # Inspired by http://stackoverflow.com/questions/12791997/how-do-you-do-a-simple-chmod-x-from-within-python os.chmod(track_fn, st.st_mode | stat.S_IROTH) track_fn = os.path.realpath(track_fn) if os.path.exists(target): if os.path.realpath(target) == os.path.realpath(track_fn): pass # Link is correct else: if args.overwrite: os.unlink(target) os.symlink(track_fn, target) # Meh, a race condition else: print('ERROR: wanted to link source \'%s\' to target \'%s\', but the target already existed and you didn\'t want to overwrite it. Use --overwrite/-ow to overwrite existing output files.\n' % (track_fn, target), file=sys.stderr) sys.exit(1) else: os.symlink(track_fn, target) # Meh, a race condition final_track_bn = track_bn final_track_ext = track_ext else: print('WARNING: track file \'%s\' does not exist and will NOT be included in the hub.\n' % (track_fn), file=sys.stderr) final_track_bn = None final_track_ext = None return(final_track_bn, final_track_ext) def _load_metadata(sample_metadata_fn, assay_metadata_fn, institution_metadata_fn): # print("SAMPLE") with open(sample_metadata_fn) as sample_metadata_fh: sample_metadata_lookup_table = _invert_dict(json.load(sample_metadata_fh)['terms']) # print("ASSAY") with open(assay_metadata_fn) as assay_metadata_fh: assay_lookup_table = _invert_dict(json.load(assay_metadata_fh)['terms']) # print("INSTITUTION") with open(institution_metadata_fn) as institution_metadata_fh: institution_lookup_table = _invert_dict(json.load(institution_metadata_fh)['terms']) return(sample_metadata_lookup_table, assay_lookup_table, institution_lookup_table) def _invert_dict(d): # key = metadata id, value = list of metadata terms d_inverse = {} # print(d.values()) for metadata_id, metadata_list in d.items(): #print(metadata_id, metadata_list, sep="\t") # FIXME I think the second term is always equal to first? Why is this? metadata_list_uniq = set(metadata_list); for metadata in metadata_list_uniq: if metadata in d_inverse.keys(): # print('WARNING: nonunique keys. Got key \'%s\' multiple times. Will use the first occurrence (\'%s\')\n' % (metadata, d_inverse[metadata]), file=sys.stderr) # This is a known issue. Will resolve in the future. TODO pass else: d_inverse[metadata] = metadata_id return(d_inverse) def _metadata2id(metadata, metadata_lookup_table): metadata_id = None #print("\n".join(metadata_lookup_table.keys()), "\n") if metadata in metadata_lookup_table.keys(): metadata_id = metadata_lookup_table[metadata] else: print('ERROR: metadata term \'%s\' did not exist in the metadata lookup table\n\n. Keys:\n%s' % (metadata, "\n".join(metadata_lookup_table.keys())), file=sys.stderr) sys.exit(1) return(metadata_id) if __name__=="__main__": import os import stat import json import urllib2 import argparse import re import getpass import errno import subprocess import math from trackDbParser import tkColor main()