#!/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 usage $ python2 gen_hub_json_twhmm.py --file methylcrf_prehub.txt --file-format1 --output-directory /bar/nrockweiler/public_html/hmm/methylcrf --assembly hg19 """ # TODO # Correct hierarchy of metadata terms. Use browser metadata term ids. 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() # Determine the species of interest species = {'hg19':'human', 'hg38':'human', 'panTro5':'chimpanzee', 'mm9':'mouse', 'rn4':'rat', 'rheMac3':'monkey', 'canFam3':'dog', 'galGal4':'chicken', 'myoLuc2':'bat' } try: species_of_interest = species[args.assembly] except KeyError as e: print('ERROR: unknown species \'%s\'' % (str(e))) raise _mk_dir_safe(args.output_directory) dat = [] required_metadata_field_names = ['sample', 'assay', 'tissue'] metadata = {field_name: [] for field_name in required_metadata_field_names} if args.file_format1: _parse_header_func_ref = _parse_format1_header _parse_datarow_func_ref = _parse_format1_datarow elif args.file_format2: _parse_header_func_ref = _parse_format2_or_format3_header _parse_datarow_func_ref = _parse_format2_datarow else: _parse_header_func_ref = _parse_format2_or_format3_header _parse_datarow_func_ref = _parse_format3_datarow found_header = False # 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 blank rows if not line: continue cols = line.split('\t') # Parse the header if line.startswith('#'): extra_metadata_field_names = _parse_header_func_ref(cols) num_extra_metadata_cols = len(extra_metadata_field_names) for field_name in extra_metadata_field_names: metadata[field_name] = [] found_header = True # Parse the data row else: if not found_header: sys.exit('ERROR: file containing list of track files (-f/--file \'%s\' does not appear to have a header. The header must be above all data rows and start with #.' % (args.file)) (sample_name, sample, tissue, assay, track_fn, extra_metadata) = _parse_datarow_func_ref(cols, species_of_interest) print(sample_name, sample, tissue, assay, track_fn, extra_metadata) if not sample_name is None: # If not all the extra metadata cols are given, fill with NA # Fill in the missing columns at the end extra_metadata += ['NA'] * (num_extra_metadata_cols - len(extra_metadata)) # From http://stackoverflow.com/questions/3438756/some-built-in-to-pad-a-list-in-python # If at least 1 middle column given, fill in from the beginning for i in xrange(num_extra_metadata_cols): if (extra_metadata[i] == ''): extra_metadata[i] = 'NA' 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 (assay == 'medip'): (final_track_bn, final_track_ext) = _process_medip(args, track_fn, track_bn, track_ext) elif (assay == 'chip' or assay == 'peak' or assay == 'input' or assay.startswith('H3') or assay == 'p300' or assay == 'NR2F1' or assay.startswith('TF')): (final_track_bn, final_track_ext) = _process_chip(args, track_fn, track_bn, track_ext) elif (assay == 'mre'): (final_track_bn, final_track_ext) = _process_mre(args, track_fn, track_bn) elif (assay == 'methylmnm'): (final_track_bn, final_track_ext) = _process_methylmnm(args, track_fn, track_bn_no_ext) elif (assay == 'methylcrf'): (final_track_bn, final_track_ext) = _process_methylcrf(args, track_fn, track_bn_no_ext) else: print('ERROR: unrecognized assay \'%s\'\n' % (assay), file=sys.stderr) sys.exit(1) if (final_track_bn is None): print('WARNING: track file \'%s\' does not exist and will NOT be included in the hub.\n' % (track_fn), file=sys.stderr) else: _add_track_data(args, sample_name, sample, assay, tissue, extra_metadata, final_track_bn, final_track_ext, dat, metadata, extra_metadata_field_names, num_extra_metadata_cols) _complete_hub(dat, metadata, required_metadata_field_names, extra_metadata_field_names, args) 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. Use --file-format1, --file-format2, --file-format3 to specify the format of the file', ) # --file-format1, --file-format2, and --file-format3 are mutually exclusive. group_format = parser.add_mutually_exclusive_group(required=True) group_format.add_argument( '--file-format1', '-f1', action='store_true', default=False, help="""File is tab-delimited with the following columns: 1) sample name, 2) sample, 3) tissue, 4) assay, 5) data file. The file can have optional metadata terms after the last requried column. The metadata terms will be displayed as extra columns in the metadata heatmap. Optional metadata terms that are missing for each sample will be replaced with a NA. 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. Example $ cat methylCRF_prehub.txt #sample specimen tissue assay file [metadata_term_1] ... TW534 TW630 ZYA23-colon1 methylCRF ZYA23-colon1 colon1 methylCRF /bar/twlab-shared/Epigenome_Evolution/human/batch_methylcrf/methylcrf_practice/ZYA23-colon1_mcrf.bed new_filtering TW535 TW631 ZYA24-colon2 methylCRF ZYA24-colon2 colon2 methylCRF /bar/twlab-shared/Epigenome_Evolution/human/batch_methylcrf/methylcrf_practice/ZYA24-colon2_mcrf.bed old_filtering\n""" ) group_format.add_argument( '--file-format2', '-f2', action='store_true', default=False, help="""File is tab-delimited with the following columns: 1) speices, 2) tissue, 3) sample name, 4) assay, 5) library ID, 6) treatment, 7) data path, 8) path to MeDIP .bigWig, 9) path to MRE .bedGraph, 10) batch, 11) exclude (i.e., the format of the master sample spreadsheet). Rows that do not match the species of interest are ignored. Rows that have a 1 in the exclude column are ignored. The order of the tracks will be in the same order as this file. Example #Species Tissue Sample Library type Library ID Treatment Data MeDIP bigwig MRE bedGraph Batch Exclude chicken brain ZYA9 brain1 MeDIP TW519 /bar/dli/HMM/C4-10-MeDIP/trim/combine /bar/dli/HMM/C4-10-MeDIP/trim/combine/TW519_ZYA9-brain1_MeDIP.bigWig NA chicken brain ZYA9 brain1 MeDIP TW726 /bar/dli/HMM/R3-MeDIP/WangT/trim/chicken/combine /bar/dli/HMM/R3-MeDIP/WangT/trim/chicken/combine/TW726_ZYA9-brain1_MeDIP.bigWig NA pool 3 """ ) group_format.add_argument( '--file-format3', '-f3', action='store_true', default=False, help="""File is tab-delimited with the following columns: 1) speices, 2) tissue, 3) sample name, 4) medip library ID, 5) mre library ID, 6) medip extended bed, 7) mre bam 8) methylcrf bed, 9) exclude (i.e., the format of the methylCRF sheet of the master sample spreadsheet). Rows that do not match the species of interest are ignored. Rows that have a 1 in the exclude column are ignored. The order of the tracks will be in the same order as this file. Example #Species Tissue Sample MeDIP library ID MRE library ID MeDIP extended bed MRE bam methylCRF bed Exclude human brain ZYA29-brain1 TW547 TW610 /bar/dli/HMM/CHF-MeDIP/trim/human/combine/TW547_ZYA29-brain1_MeDIP.extended.bed /bar/dli/HMM/RHMDM-MRE/trim/human/combine/TW610_ZYA29-brain1_MRE.bam /bar/twlab-shared/Epigenome_Evolution/human/batch_methylcrf/methylcrf_practice/ZYA29-brain1_mcrf.bed human brain ZYA29-brain1 TW644 TW610 /bar/dli/HMM/K562-ZYA29-MeDIP/trim/combine/TW644_ZYA29-50ng-diagenAgarose_MeDIP.extended.bed /bar/dli/HMM/RHMDM-MRE/trim/human/combine/TW610_ZYA29-brain1_MRE.bam /bar/twlab-shared/Epigenome_Evolution/human/batch_methylcrf/TW644_TW610_brain1_mcrf/ZYA29-brain1_TW644_TW610_mcrf.bed """ ) 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( '--output-hub', '-oh', default='hub', help='basename of output hub file. The hub file is created in the output directory.', ) parser.add_argument( '--assembly', '-a', default='hg19', choices=['hg19', 'hg38', 'panTro5', 'mm9', 'rn4', 'rheMac3', 'canFam3', 'galGal4', 'myoLuc2'], help='genome assembly for tracks. Used for generating the browser link. Option is CASE SENSITIVE.', ) # --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 is 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.', ) parser.add_argument( '--increment-hub', '-i', action='store_true', default=False, help='specify option to increment output hub to the next version, e.g., if hubV1 exists, then hubV2 will be created. NOTE: if the hub already exists --overwrite trumps --increment-hub, and --increment-hub trumps --easy-going.' ) args = parser.parse_args() args.output_directory = os.path.realpath(args.output_directory) # Construct URL 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) args.output_hub = os.path.join(args.output_directory, args.output_hub) 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 _increment_hub(hub_f): """Determine the name of the next hub Args: hub_f (str): hub filename Returns: hub_f (str): next hub filename """ while(os.path.exists(hub_f)): # Figure out which version the hub is at m = re.search('(?PV\d+)?(?P\.json)?$', hub_f) if m: if m.group('version') is None: iteration = 2 else: iteration = int(m.group('version')[1:]) + 1 # Remove the initial V hub_f = '%sV%d' % (hub_f[:m.start()], iteration) if not m.group('extension') is None: hub_f += '.json' else: print('ERROR: could not parse version from hub file \'%s\'\n' % (hub_f), file=sys.stderr) sys.exit(1) return(hub_f) def _can_create_output_f(output_f, overwrite, easy_going, increment_hub=False): """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: if not increment_hub: create_output_f = False print('ERROR: output file \'%s\' already exists and you didn\'t want to overwrite it, use the existing file, or (possibly) increment the hub file. Use --overwrite/-ow to overwrite existing output files, --easy-going/-e to use existing files or possibly --increment-hub/-i to create the next version of the file.\n' % (output_f), file=sys.stderr) sys.exit(1) elif easy_going: if not increment_hub: create_output_f = False print('INFO: output file \'%s\' already exists and you wanted to keep it. Use --overwrite to overwrite the file or possibly --increment-hub/-i to create the next version of 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): """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 _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): """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 _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 _add_track_data(args, sample_name, sample, assay, tissue, extra_metadata, final_track_bn, final_track_ext, dat, metadata, extra_metadata_field_names, num_extra_metadata_cols): """Update the metadata and dat data structures Args: args (ArgumentParser object): command-line arguments sample_name (str): name of sample, e.g., ZYA29-brain1 TW547 MeDIP sample (str): specimen, e.g., ZYA29-brain1. Used as the "sample" metadata term assay (str): assay type, e.g., mre tissue (str): tissue type, e.g., colon extra_metadata (list): list of optional metadata fields. Empty if no optional metadata final_track_bn (str): track basename final_track_ext (str): track extension dat (dict): hub data metadata (dict): hub metadata extra_metadata_field_names (list): list of optional metadata field names (for the hub) num_extra_metadata_cols (int): number of optional metadata field names (for the hub) Returns: Nothing. However, metadata and dat are modified. """ i = {} i['name'] = sample_name i['url'] = '/'.join([args.url, final_track_bn]) i['type'] = final_track_ext i['mode'] = 'show' i['metadata'] = [sample, assay, tissue] + extra_metadata if assay == 'medip': if final_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 } else: try: col = tkColor[assay] except: try: col = tkColor[tissue] except: col = '0,77,0' i['colorpositive'] = 'rgb({})'.format(col) i['height'] = 30 metadata['sample'].append(sample) metadata['assay'].append(assay) metadata['tissue'].append(tissue) for j in xrange(num_extra_metadata_cols): metadata_field_name = extra_metadata_field_names[j] metadata_field_value = extra_metadata[j] metadata[metadata_field_name].append(metadata_field_value) # If the track is of old data, hide track if (metadata_field_name == 'batch'): if (metadata_field_value == 'old_library'): i['mode'] = 'hide' dat.append(i) def _complete_hub(dat, metadata, required_metadata_field_names, extra_metadata_field_names, args): """Finalize hub: add native tracks (e.g., gene track), add metadata, and write the hub to file Args: dat (dict): hub data metadata (dict): hub metadata required_metadata_field_names (list): list of required metadata field names (for the hub) extra_metadata_field_names (list): list of optional metadata field names (for the hub) args (ArgumentParser object): command-line arguments Returns: Nothing. Output hub file is created. """ # 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'] = {} for field_name, field_values in metadata.iteritems(): metadata[field_name] = list(set(field_values)) metadata[field_name].sort() cm['vocabulary'][field_name] = metadata[field_name] cm['show'] = required_metadata_field_names + extra_metadata_field_names # Preserves order dat.append(cm) # Create dathub file if _can_create_output_f(args.output_hub, args.overwrite, args.easy_going, args.increment_hub): if os.path.exists(args.output_hub): if not args.overwrite: if args.increment_hub: args.output_hub = _increment_hub(args.output_hub) with open(args.output_hub, 'w') as hub_fh: json.dump(dat, hub_fh, sort_keys=True) args.output_hub_url = '/'.join([args.url, os.path.basename(args.output_hub)]) 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 \'https://epigenomegateway.wustl.edu/browser/?genome=%s&datahub=%s\'\n' % (args.assembly, args.output_hub_url), file=sys.stderr) def _process_chip(args, track_fn, track_bn, track_ext): """ Same as medip """ final_track_bn, final_track_ext = _process_medip(args, track_fn, track_bn, track_ext) return(final_track_bn, final_track_ext) def _process_medip(args, track_fn, track_bn, track_ext): """Process MeDIP track: file should be a .bigwig (binary). Softlink file to output directory Args: args (ArgumentParser object): command-line arguments track_fn (str): track filename track_bn (str): track basename track_ext (str): track extension Returns: final_track_bn (str): track basename. (None if output file could not be created.) final_track_ext (str): track extension. (None if output file could not be created.) """ final_track_bn = None final_track_ext = None target = os.path.join(args.output_directory, track_bn) if os.path.exists(track_fn): 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 return(final_track_bn, final_track_ext) def _process_mre(args, track_fn, track_bn): """Process MRE track: file should be .bedgraph. bgzip and tabix index the file in the output directory Args: args (ArgumentParser object): command-line arguments track_fn (str): track filename track_bn (str): track basename Returns: final_track_bn (str): track basename. (None if output file could not be created.) final_track_ext (str): track extension. (None if output file could not be created.) """ final_track_bn = None final_track_ext = None target = os.path.join(args.output_directory, track_bn) if os.path.exists(track_fn): 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 # Since the mre files are in d. lis analysis directory, simply bziping that is not possible 1) no write permission 2) this would remove the uncompressed file (unless -c and the output was redirected. A hack around this is to link the file to the working directory, and bgzip it. After compressed, the link will be removed by bzgip. # Link the bedgraph to the working directory # bzip track_bedgraph_bzip_fn = _create_bzip(target, args.overwrite, args.easy_going) 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) return(final_track_bn, final_track_ext) def _process_methylmnm(args, track_fn, track_bn_no_ext): """Process methylMnM track: convert to bedgraph, bgzip, and tabix index in the output directory Args: args (ArgumentParser object): command-line arguments track_fn (str): track filename track_bn_no_ext (str): track basename sans extension Returns: final_track_bn (str): track basename. (None if output file could not be created.) final_track_ext (str): track extension. (None if output file could not be created.) """ # Create bedgraph track_bedgraph_fn = os.path.join(args.output_directory, track_bn_no_ext + '.bedgraph') _methylmnm2bedgraph(track_fn, track_bedgraph_fn, args.overwrite, args.easy_going) # bzip track_bedgraph_bzip_fn = _create_bzip(track_bedgraph_fn, args.overwrite, args.easy_going) 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) return(final_track_bn, final_track_ext) def _process_methylcrf(args, track_fn, track_bn_no_ext): """Process methylCRF track: convert to bedgraph, bgzip, and tabix index in the output directory Args: args (ArgumentParser object): command-line arguments track_fn (str): track filename track_bn_no_ext (str): track basename sans extension Returns: final_track_bn (str): track basename. (None if output file could not be created.) final_track_ext (str): track extension. (None if output file could not be created.) """ # Create bedgraph track_bedgraph_fn = os.path.join(args.output_directory, track_bn_no_ext + '.bedgraph') 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) 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) return(final_track_bn, final_track_ext) def _parse_format1_header(cols): """Determine the optional metadata field names in a format1 file Args: cols (list): list of column headers Returns: extra_metadata_field_names (list): list of optional metadata field names (for the hub). Empty if not optional metadata. """ # Determine the extra metadata column names (if any) extra_metadata_field_names = [] if (len(cols) > 4): extra_metadata_field_names = cols[5:] return(extra_metadata_field_names) def _parse_format1_datarow(cols, species_of_interest): # dummy species_of_interest to keep the two datarow functions have the same number of arguments """Determine the sample information in a format1 datarow Args: cols (list): list of values in a data row species_of_interest (str): species of interest. (Dummy variable to keep the two datarow functions have the same number of arguments.) Returns: sample_name (str): name of sample, e.g., ZYA23-colon1 TW534 MeDIP sample (str): library ID, e.g., ZYA23-colon1 assay (str): assay type, e.g., mre tissue (str): tissue type, e.g., colon extra_metadata (list): list of optional metadata fields. Empty if no optional metadata """ sample_name = cols[0] sample = cols[1] tissue = cols[2] assay = cols[3] track_fn = cols[4] extra_metadata = cols[5:] # Will be [] if no extra metadata terms are given return(sample_name, sample, tissue, assay, track_fn, extra_metadata) def _parse_format2_or_format3_header(cols): """Determine the optional metadata field names in a format2 or format3 file Args: cols (list): list of column headers. (Dummy variable to keep the two header functions have the same number of arguments.) Returns: extra_metadata_field_names (list): list of optional metadata field names (for the hub). Empty if not optional metadata. """ extra_metadata_field_names = ['library_batch'] return(extra_metadata_field_names) def _parse_format2_datarow(cols, species_of_interest): """Determine the sample information in a format2 datarow Args: cols (list): list of values in a data row species_of_interest (str): species of interest. Returns: sample_name (str): name of sample, e.g., ZYA23-colon1 TW534 MeDIP sample (str): library ID, e.g., ZYA23-colon1 assay (str): assay type, e.g., mre tissue (str): tissue type, e.g., colon extra_metadata (list): list of optional metadata fields. Empty if no optional metadata """ sample_name = None sample = None tissue = None assay = None track_fn = None extra_metadata = [] species = cols[0].lower() # e.g., human try: exclude = int(cols[10]) except ValueError as e: exclude = 0 # Skip if the species is not the species of interest # Skip if exclude is true if (species == species_of_interest): if (exclude != 1): tissue = cols[1].lower() # e.g., colon sample = cols[2] # e.g., ZYA23-colon1 assay = cols[3] # e.g., MRE library_id = cols[4] # e.g., TW534 treatment = cols[5] data = cols[6] medip_bigwig = cols[7] mre_bedgraph = cols[8] batch = cols[9] # Rename specimen from 'ZYA25-liver1' to 'ZYA25-liver1 TW534 MRE' sample_name = " ".join([sample, library_id, assay]) assay = assay.lower() # Add new/old library batch information if ((library_id == 'TW644') or (library_id.startswith('TW7'))): extra_metadata = ['new_library'] else: extra_metadata = ['old_library'] # Grab the write path variable if (assay == 'medip'): track_fn = medip_bigwig elif (assay == 'mre'): track_fn = mre_bedgraph else: sys.exit('ERROR: unknown assay \'%s\'.' % (assay)) return(sample_name, sample, tissue, assay, track_fn, extra_metadata) def _parse_format3_datarow(cols, species_of_interest): """Determine the sample information in a format3 datarow Args: cols (list): list of values in a data row species_of_interest (str): species of interest. Returns: sample_name (str): name of sample, e.g., ZYA23-colon1 TW534 MeDIP sample (str): library ID, e.g., ZYA23-colon1 assay (str): assay type, e.g., mre tissue (str): tissue type, e.g., colon extra_metadata (list): list of optional metadata fields. Empty if no optional metadata """ sample_name = None sample = None tissue = None assay = "methylcrf" track_fn = None extra_metadata = [] species = cols[0].lower() # e.g., human try: exclude = int(cols[8]) except ValueError as e: exclude = 0 # Skip if the species is not the species of interest # Skip if exclude is true if (species == species_of_interest): if (exclude != 1): tissue = cols[1].lower() # e.g., colon sample = cols[2] # e.g., ZYA23-colon1 medip_library_id = cols[3] # e.g., TW534 mre_library_id = cols[4] # e.g., TW600 medip_extended_bed = cols[5] mre_bam = cols[6] methylcrf_bed = cols[7] # Rename specimen from 'ZYA25-liver1' to 'ZYA25-liver1 TW534 TW600 methylCRF' sample_name = " ".join([sample, medip_library_id, mre_library_id, 'methylcrf']) # Add new/old library batch information if ((medip_library_id == 'TW644') or (medip_library_id.startswith('TW7'))): extra_metadata = ['new_library'] else: extra_metadata = ['old_library'] track_fn = methylcrf_bed return(sample_name, sample, tissue, assay, track_fn, extra_metadata) if __name__=="__main__": import os import json import urllib2 import argparse import re import getpass import errno import subprocess import math from trackDbParser import tkColor FTP_SITE = 'https://wangftp.wustl.edu' main()