# Usage: python make_parameter_bw.py # Edit the variable BW_PATHS, which is a file containing (almost all) paths to in-house bw files. Edit PARENT_DIR, which # will be used to make full path to the bw files. Amend name_map_dict. # Purpose # ------- # Make a parameter file for buildCmatrix.py. For example: # # @C0001100007141 # PROTEIN_ID: # PROTEIN_NAME:HTR13 # DATA_NAME:2_71C_10_IP-HTR13-17c # DATA_FORMAT:bw # DESCRIPTION:in house ChIP data # LOCATION:/media/pw_synology3/PW_HiSeq_data/ChIP-seq/Mapped_data/2_71C/20150603_H3_2nd_rep_mapped_bw/ChIP-seq_h31_h33_col_rep1_20150607_max_fragment_500_bw/2_71C_10_IP-HTR13-17c_raw_trimmo_paired_truseq3-PE-2_2_10_5_1_bowtie2_TAIR10_ensembl_nomixed_sorted_rmdup_picard_genomenorm.bw # NOTE: # # 29 NOV 2016, hui, home import sys, glob, os, operator from geneid2name import make_gene_name_AGI_map_dict BW_PATHS = 'bwfiles_unique.txt' # Ideally, the input file should be sorted by library number, then sample number. PARENT_DIR = '/media/pw_synology3/PW_HiSeq_data/ChIP-seq/Mapped_data' GENE_ID_TO_GENE_NAME = '../Data/information/AGI-to-gene-names_v2.txt' def get_library_id(s): index = s.find('C') # in-house chip library id ends with a letter C if index > 0: return s[:index] else: return s def get_bw_file_name(lst): for x in lst: if '.bw' in x: return x return '' def get_sample_name(s): index = s.find('_raw') # a bw file name contains _raw, get the part before _raw if _raw exists if index > 0: return s[:index] else: return s def get_sample_number(s): index = s.find('_S') # sample number is proceeded by _S if index > 0: if s[index+2].isdigit() and s[index+3].isdigit(): # two digit sample number return s[index+2:index+4] else: # one digit sample number return '0' + s[index+2] return '25' def convert_name(s, d): ''' convert s to gene id ''' if s.upper() in d: return d[s.upper()] else: return ' ' def get_gene_id_and_name(s, d): ''' Return a dictionary value, protein name and protein id given a sample file name.''' for k in sorted(d.keys(),reverse=True): if k.isdigit() and k.lower() in s.lower(): # k is a number, e.g., 833 return (d[k], convert_name(d[k], agi2name_dict)) if k.lower() in s.lower(): return (k, convert_name(k, agi2name_dict)) return ('name_unknown', 'id_unknown') ### # key is a name of sample, value is either empty if it is a protein name, or the protein name if the key is a database number, e.g., 833. If key is protein name, the program tries to find its gene id (AT...). If key is a database number, the program gets its protein name and then tries to find its gene id. If not found, PROTEIN_ID or PROTEIN_NAME will be assigned _unknown. So search _unknown to manually annotate PROTEIN_ID or PROTEIN_NAME. name_map_dict = {'HTA11':'', 'H3':'', 'HTR13':'', 'HSP70':'', 'KIN10':'', 'EPM1':'', 'ELF3':'', 'phyB':'', '833':'KIN10', 'hos':'', 'HOS':'', 'HOS1':'', 'LUX':'', 'YHB':'','HSF1':'','3129':'HSP90','838':'phyB', '1506':'HSF1', 'SUMO':'', 'TIR1':'', 'PIF4':'', 'PIF':'','H2A':'', 'H2AZ':'', 'MNase':'', '544':'ELF4', '834':'PIF4', '745':'ELF3', 'EC1_S1':'', 'EC1_S1':'ELF3', 'EC2_S2':'ELF3', '1166':'LUX', '1167':'LUX', '3239':'REV', '1281':'MPK6', '1278':'SEP3', '1279':'FD', '1280':'FD-like', '1283':'HSP70', '1284':'DELLA', '1762':'FT', 'LFY':'', 'TFL1':''} agi2name_dict = make_gene_name_AGI_map_dict(GENE_ID_TO_GENE_NAME) f = open(BW_PATHS) d_cid = {} for line in f: line = line.strip() lst = line.split('/') bwfile = get_bw_file_name(lst) sample_name = get_sample_name(bwfile) sample_no = get_sample_number(bwfile) library_id = get_library_id(lst[0]) name, gid = get_gene_id_and_name(sample_name, name_map_dict) path = os.path.join(PARENT_DIR, line) if '2_' in library_id: # some in-house library id starts with 2_, indicating that it is a replicate. index = library_id.find('_') library_id = library_id[index+1:] cid = '@C00011000%s%s' % (library_id.zfill(3), sample_no) # the highest digit of the last 9 digits is 1 to indicate that it is a replicate. library number, 3 digits. sample number, 2 digits. else: cid = '@C00010000%s%s' % (library_id.zfill(3), sample_no) # handle duplicate library id if not cid in d_cid: d_cid[cid] = 1 else: # produce a new library id d_cid[cid] += 1 count = 26 head = cid[-2:] # last two digits while True: new_cid = cid[:-2] + '%02d' % (count) if count >= 100: print('Error: count must be less than 100.') sys.exit() count += 1 if not new_cid in d_cid: break cid = new_cid d_cid[cid] = 1 # print contents for the parameter file print(cid) print('PROTEIN_ID:%s' % (gid)) print('PROTEIN_NAME:%s' % (name)) print('DATA_NAME:%s' % (sample_name)) print('DATA_FORMAT:%s' % ('bw')) print('DESCRIPTION:in house ChIP data') print('LOCATION:%s' % (path)) print('NOTE:') print('') f.close()