Mercurial > repos > alenail > map_to_known_genes
comparison map_to_known_genes/map_to_known_genes.py @ 0:3f12a2b0af50 draft default tip
Uploaded
author | alenail |
---|---|
date | Thu, 14 Apr 2016 11:34:47 -0400 |
parents | |
children |
comparison
equal
deleted
inserted
replaced
-1:000000000000 | 0:3f12a2b0af50 |
---|---|
1 #!/usr/local/bin/python | |
2 | |
3 import sys, os | |
4 from optparse import OptionParser | |
5 from collections import defaultdict as dd | |
6 from csv import DictReader, DictWriter | |
7 | |
8 from chipsequtil import MACSFile, BEDFile, KnownGeneFile, parse_number | |
9 from chipsequtil.util import MultiLineHelpFormatter | |
10 | |
11 usage = '%prog [options] <knownGene file> <knownGene xRef file> <peaks file>' | |
12 description = """ | |
13 Map the peaks in <peaks file> to genes in <knownGene file>. <knownGene file> is\ | |
14 format is as specified in http://hgdownload.cse.ucsc.edu/goldenPath/hg19/database/knownGene.sql.\ | |
15 <peaks file> format is as produced by MACS. If *auto* is chosen (default) file extension \ | |
16 is examined for *.xls* for default MACS format or *.bed* for BED format. If the --detail \ | |
17 option is provided, the following extra fields are appended to each row: | |
18 | |
19 peak loc, dist from feature, map type, map subtype | |
20 """ | |
21 epilog = '' | |
22 parser = OptionParser(usage=usage,description=description,epilog=epilog,formatter=MultiLineHelpFormatter()) | |
23 parser.add_option('--upstream-window',dest='upst_win',type='int',default=5500,help='window width in base pairs to consider promoter region [default: %default]') | |
24 parser.add_option('--downstream-window',dest='dnst_win',type='int',default=2500,help='window width in base pairs to consider downstream region [default: %default]') | |
25 parser.add_option('--tss',dest='tss',action='store_true',help='calculate downstream window from transcription start site instead of transcription end site') | |
26 parser.add_option('--map-output',dest='peak_output',default=None,help='filename to output mapped peaks to [default: stdout]') | |
27 parser.add_option('--stats-output',dest='stats_output',default=sys.stderr,help='filename to output summary stats in conversion [default: stderr]') | |
28 parser.add_option('--peaks-format',dest='peaks_fmt',default='auto',type='choice',choices=['auto','MACS','BED'],help='format of peaks input file [default: %default]') | |
29 parser.add_option('--detail',dest='detail',action='store_true',help='add extra fields to output, see description') | |
30 parser.add_option('--intergenic',dest='intergenic',action='store_true',help='write intergenic peaks to the gene file as well with None as gene ID') | |
31 #parser.add_option('--symbol-xref',dest='symbol_xref',default=None,help='use the kgXref table file supplied to find a gene symbol, output as second column') | |
32 | |
33 # TODO - options | |
34 #parser.add_option('--use-cds',dest='use_cds',action='store_true',help='use cdsStart and cdsEnd fields instead of txStart and txEnd to do mapping') | |
35 #parser.add_option('--capture-intergenic'...) | |
36 #parser.add_option('--map-format',dest='peak_format',type='choice',choices=['default','BED'],help='format of peak output [default: %default]') | |
37 #parser.add_option('--stats-format',dest='stats_format',type='choice',choices=['human','python'],help='format of summary stats output [default: %default]') | |
38 | |
39 def parse_gene_ref(ref_gene) : | |
40 reader = KnownGeneFile(ref_gene) | |
41 gene_ref = dd(list) | |
42 for ref_dict in reader : | |
43 gene_ref[ref_dict['chrom']].append(ref_dict) | |
44 | |
45 return gene_ref | |
46 | |
47 | |
48 def parse_gene_ref_line(l) : | |
49 l = map(parse_number, l) # coerce to numbers where possible | |
50 l[9] = map(parse_number, l[9].split(',')) # turn 'x,x,x,...' into list | |
51 l[10] = map(parse_number, l[10].split(',')) | |
52 return l | |
53 | |
54 | |
55 def main(): | |
56 opts, args = parser.parse_args(sys.argv[1:]) | |
57 | |
58 if len(args) < 3 : | |
59 parser.error('Must provide three filename arguments') | |
60 | |
61 gene_ref = parse_gene_ref(args[0]) | |
62 xref_fn = args[1] | |
63 peaks_fn = args[2] | |
64 | |
65 if opts.peaks_fmt == 'MACS' : | |
66 peaks_reader_cls = MACSFile | |
67 chr_field, start_field, end_field = 'chr', 'start', 'end' | |
68 elif opts.peaks_fmt == 'BED' : | |
69 peaks_reader_cls = BEDFile | |
70 chr_field, start_field, end_field = 'chrom', 'chromStart', 'chromEnd' | |
71 else : | |
72 # should never happen | |
73 fieldnames = [] | |
74 | |
75 #peaks_reader = DictReader(open(args[1]),fieldnames=fieldnames,delimiter='\t') | |
76 peaks_reader = peaks_reader_cls(peaks_fn) | |
77 | |
78 # default output format: | |
79 if opts.peak_output : | |
80 peak_output = open(opts.peak_output,'w') | |
81 else : | |
82 peak_output = sys.stdout | |
83 | |
84 fieldnames = peaks_reader.FIELD_NAMES | |
85 if opts.detail : | |
86 fieldnames += ["peak loc","dist from feature","map type","map subtype"]#"score" | |
87 output_fields = ['knownGeneID']+fieldnames | |
88 | |
89 # see if the user wants gene symbols too | |
90 # TODO - actually make this an option, or make it required | |
91 opts.symbol_xref = xref_fn | |
92 if opts.symbol_xref : | |
93 kgXref_fieldnames = ['kgID','mRNA','spID','spDisplayID','geneSymbol','refseq','protAcc','description'] | |
94 symbol_xref_reader = DictReader(open(opts.symbol_xref),fieldnames=kgXref_fieldnames,delimiter='\t') | |
95 symbol_xref_map = {} | |
96 for rec in symbol_xref_reader : | |
97 symbol_xref_map[rec['kgID']] = rec | |
98 output_fields = ['knownGeneID','geneSymbol']+fieldnames | |
99 | |
100 peaks_writer = DictWriter(peak_output,output_fields,delimiter='\t',extrasaction='ignore',lineterminator='\n') | |
101 peaks_writer.writerow(dict([(k,k) for k in output_fields])) | |
102 unique_genes = set() | |
103 map_stats = dd(int) | |
104 for peak in peaks_reader : | |
105 | |
106 # if this is a comment or header line get skip it | |
107 if peak[fieldnames[0]].startswith('#') or \ | |
108 peak[fieldnames[0]] == fieldnames[0] or \ | |
109 peak[fieldnames[0]].startswith('track') : continue | |
110 | |
111 # coerce values to numeric if possible | |
112 for k,v in peak.items() : peak[k] = parse_number(v) | |
113 | |
114 # MACS output gives us summit | |
115 if opts.peaks_fmt == 'MACS' : | |
116 peak_loc = peak[start_field]+peak['summit'] | |
117 else : # peak assumed to be in the middle of the reported peak range | |
118 peak_loc = (peak[start_field]+peak[end_field])/2 | |
119 | |
120 chrom_genes = gene_ref[peak[chr_field]] | |
121 | |
122 if len(chrom_genes) == 0 : | |
123 sys.stdout.write('WARNING: peak chromosome %s not found in gene reference, skipping: %s\n'%(peak[chr_field],peak)) | |
124 continue | |
125 | |
126 mapped = False | |
127 | |
128 # walk through the genes for this chromosome | |
129 for gene in chrom_genes : | |
130 | |
131 # reusable dictionary for output | |
132 out_d = {}.fromkeys(output_fields,0) | |
133 out_d.update(peak) | |
134 out_d['map type'] = '' | |
135 out_d['chromo'] = peak[chr_field] | |
136 out_d['peak loc'] = peak_loc | |
137 | |
138 # determine intervals for promoter, gene, and downstream | |
139 if gene['strand'] == '+' : | |
140 promoter_coords = max(gene['txStart']-1-opts.upst_win,0), gene['txStart']-1 | |
141 if opts.tss : | |
142 gene_coords = gene['txStart'], min(gene['txEnd'],gene['txStart']+opts.dnst_win) | |
143 downstream_coords = gene['txEnd']+1,gene['txStart']+opts.dnst_win | |
144 else : | |
145 gene_coords = gene['txStart'], gene['txEnd'] | |
146 downstream_coords = gene['txEnd']+1, gene['txEnd']+1+opts.dnst_win | |
147 else : | |
148 promoter_coords = gene['txEnd']+1, gene['txEnd']+1+opts.upst_win # +1 because we're using 1 based indexing | |
149 if opts.tss : | |
150 gene_coords = max(gene['txStart'],gene['txEnd']-opts.upst_win), gene['txEnd'] | |
151 downstream_coords = gene['txEnd']-1-opts.dnst_win, gene['txStart']-1 # -1 because we're using 1 based indexing | |
152 else : | |
153 gene_coords = gene['txStart'], gene['txEnd'] | |
154 downstream_coords = gene['txStart']-1-opts.dnst_win, gene['txStart']-1 # -1 because we're using 1 based indexing | |
155 | |
156 # check for promoter | |
157 if peak_loc >= promoter_coords[0] and peak_loc <= promoter_coords[1] : | |
158 out_d['map type'] = 'promoter' | |
159 out_d['dist from feature'] = peak_loc - promoter_coords[1] if gene['strand'] == '+' else promoter_coords[0] - peak_loc | |
160 | |
161 # check for gene | |
162 elif peak_loc >= gene_coords[0] and peak_loc <= gene_coords[1] : | |
163 # check for intron/exon | |
164 exon_coords = zip(gene['exonStarts'],gene['exonEnds']) | |
165 in_exon = False | |
166 for st,en in exon_coords : | |
167 if peak_loc >= st and peak_loc <= en : | |
168 in_exon = True | |
169 break | |
170 out_d['map type'] = 'gene' | |
171 out_d['map subtype'] = 'exon' if in_exon else 'intron' | |
172 | |
173 #Commented out to keep score reported in bed file - AJD 7/29/14 | |
174 # score = (peak-TSS)/(TSE-TSS) - peak distance from TSS as fraction of length of gene | |
175 #gene_len = float(gene_coords[1]-gene_coords[0]) | |
176 #out_d['score'] = (peak_loc-gene_coords[0])/gene_len if gene['strand'] == '+' else (gene_coords[1]-peak_loc)/gene_len | |
177 | |
178 # distance calculated from start of gene | |
179 out_d['dist from feature'] = peak_loc - promoter_coords[1] if gene['strand'] == '+' else promoter_coords[0] - peak_loc | |
180 | |
181 map_stats[out_d['map subtype']] += 1 | |
182 | |
183 # check for downstream | |
184 elif peak_loc >= downstream_coords[0] and peak_loc <= downstream_coords[1] : | |
185 out_d['map type'] = 'after' | |
186 if opts.tss : | |
187 out_d['dist from feature'] = peak_loc - gene_coords[0] if gene['strand'] == '+' else gene_coords[1] - peak_loc | |
188 else : | |
189 out_d['dist from feature'] = peak_loc - downstream_coords[0] if gene['strand'] == '+' else downstream_coords[1] - peak_loc | |
190 | |
191 # does not map to this gene | |
192 else : | |
193 pass | |
194 | |
195 # map type is not blank if we mapped to something | |
196 if out_d['map type'] != '' : | |
197 | |
198 #out_d = {'knownGeneID':gene['name']} | |
199 out_d['knownGeneID'] = gene['name'] | |
200 if opts.symbol_xref : | |
201 out_d['geneSymbol'] = symbol_xref_map[gene['name']]['geneSymbol'] | |
202 peaks_writer.writerow(out_d) | |
203 mapped = True | |
204 | |
205 # reset map_type | |
206 out_d['map type'] = '' | |
207 | |
208 if not mapped : | |
209 if opts.intergenic : | |
210 out_d['knownGeneID'] = 'None' | |
211 out_d['geneSymbol'] = 'None' | |
212 out_d['map type'] = 'intergenic' | |
213 peaks_writer.writerow(out_d) | |
214 map_stats['intergenic'] += 1 | |
215 | |
216 if peak_output != sys.stdout: | |
217 peak_output.close() | |
218 | |
219 #if opts.stats_output != sys.stderr : | |
220 # opts.stats_output = open(opts.stats_output,'w') | |
221 | |
222 #for k,v in map_stats.items() : | |
223 # opts.stats_output.write('%s: %s\n'%(k,v)) | |
224 | |
225 #if opts.stats_output != sys.stderr : | |
226 # opts.stats_output.close() | |
227 | |
228 | |
229 if __name__ == '__main__' : | |
230 main() |