|
0
|
1 '''
|
|
|
2 NAME
|
|
|
3
|
|
|
4 pieplots_macs.py
|
|
|
5
|
|
|
6 SYNOPSIS
|
|
|
7
|
|
|
8 python pieplots_macs.py --genefile MACSoutfile_genes.txt --peakfile MACSoutfile_peaks.bed --outfile MACSdirectory/piechart.pdf
|
|
|
9
|
|
|
10
|
|
|
11 DESCRIPTION
|
|
|
12
|
|
|
13 Peaks are assigned to the closest gene and then categorized according to their location at different genomic regions (promoter, intron, exon, or after the gene). Sites >10kb away from any gene are considered intergenic. (from Pamela)
|
|
|
14
|
|
|
15 '''
|
|
|
16
|
|
|
17 __author__='Renan Escalante'
|
|
|
18 __email__='renanec@mit.edu'
|
|
|
19
|
|
2
|
20 import sys
|
|
0
|
21 import matplotlib
|
|
|
22 matplotlib.use('pdf')
|
|
|
23 from matplotlib import pyplot as plt
|
|
|
24 matplotlib.rcParams['pdf.fonttype']=42
|
|
|
25 matplotlib.rcParams['font.size']=14
|
|
5
|
26 # import pandas as pd
|
|
0
|
27 from argparse import ArgumentParser
|
|
|
28
|
|
|
29 def map_peaks(gene,peak,outfile,macsFlag):
|
|
|
30 genefile = open(gene, 'r')
|
|
|
31 peakfile = open(peak, 'r')
|
|
|
32
|
|
|
33 types = {'promoter':0, 'after':0, 'intron':0, 'exon': 0}
|
|
|
34
|
|
|
35 #read mapped gene file, store closest map for each peak
|
|
|
36 peaks={} #{chrom:{peakStart:[dist, type]}}
|
|
|
37 for line in genefile:
|
|
|
38 words = line.strip().split('\t')
|
|
|
39 #ignore first line
|
|
|
40 if words[0] == 'knownGeneID': continue
|
|
|
41 chrom = words[2]
|
|
|
42
|
|
6
|
43 if macsFlag:
|
|
7
|
44 print "macsflag on"
|
|
6
|
45 start = int(words[3])-1
|
|
|
46 dist = abs(int(words[12]))
|
|
7
|
47 print dist
|
|
6
|
48 maptype = words[14]
|
|
|
49 if maptype == 'gene':
|
|
|
50 maptype = words[15]
|
|
|
51
|
|
|
52 else:
|
|
7
|
53 print "macsflag off"
|
|
0
|
54 try:
|
|
8
|
55 print words
|
|
0
|
56 start = int(words[3])
|
|
|
57 dist = abs(int(words[15]))
|
|
7
|
58 print dist
|
|
0
|
59 maptype = words[16]
|
|
|
60 if maptype == 'gene':
|
|
|
61 maptype = words[17]
|
|
|
62 except:
|
|
7
|
63 print 'trigger except'
|
|
0
|
64 pass
|
|
|
65
|
|
6
|
66 print "dist:"
|
|
|
67 print dist
|
|
7
|
68 print "that was dist"
|
|
0
|
69
|
|
|
70 if chrom not in peaks:
|
|
|
71 #new chrom
|
|
|
72 peaks[chrom] = {start:[dist,maptype]}
|
|
|
73 else:
|
|
|
74 if start in peaks[chrom].keys():
|
|
|
75 #account for duplicate entries - choose closest gene and store type
|
|
|
76 if dist < peaks[chrom][start][0]:
|
|
|
77 #closer gene
|
|
|
78 peaks[chrom][start] = [dist, maptype]
|
|
|
79 else: peaks[chrom][start] = [dist, maptype]
|
|
|
80
|
|
|
81 #count types - 1 per peak in peak file
|
|
|
82 types = {'promoter':0, 'after':0, 'intron':0, 'exon': 0, 'inter': 0}
|
|
|
83 totalpks = 0
|
|
|
84 #Read peak file in bed format
|
|
|
85 for line in peakfile:
|
|
|
86 words = line.strip().split('\t')
|
|
|
87 chrom = words[0]
|
|
|
88 start = int(words[1])
|
|
|
89 if chrom in peaks:
|
|
|
90 if start in peaks[chrom]:
|
|
|
91 types[peaks[chrom][start][1]] += 1
|
|
|
92 else:
|
|
|
93 types['inter'] += 1
|
|
|
94 else:
|
|
|
95 types['inter'] += 1
|
|
|
96 totalpks += 1
|
|
|
97
|
|
|
98
|
|
|
99 #--------------------------------------------
|
|
|
100 # make a square figure and axes
|
|
|
101 #--------------------------------------------
|
|
|
102
|
|
|
103 fig = plt.figure(figsize=(6,6))
|
|
|
104 pie_ax = fig.add_axes((0.3,0.2,0.4,0.4))
|
|
|
105
|
|
|
106 # The slices will be ordered and plotted counter-clockwise.
|
|
|
107 labels = ['exon: %i'%types['exon'],'intron: %i'%types['intron'],'promoter: %i'%types['promoter'],'intergenic: %i'%types['inter'], 'after: %i'%types['after']]
|
|
|
108 fracs = [types['exon'], types['intron'], types['promoter'], types['inter'], types['after']]
|
|
|
109
|
|
|
110 plt.pie(fracs, labels=labels) #, autopct='%1.1f%%')
|
|
|
111
|
|
|
112 #Export data frame with all the counts
|
|
|
113 indexDataFrame = ['exon','intron','promoter','intergenic','after']
|
|
|
114 df = pd.DataFrame(data=fracs, index=indexDataFrame)
|
|
|
115 dfFileName = outfile.replace("pdf","csv")
|
|
|
116 df.to_csv(dfFileName, sep=',')
|
|
|
117 #plt.title('MACS peaks in %s'%(name))
|
|
|
118 plt.figtext(.5, .1, 'Total: %i'%totalpks, ha='center')
|
|
|
119 fig.savefig(outfile)
|
|
|
120
|
|
|
121 def main():
|
|
|
122 usage = "usage: %prog --genefile MACSoutfile_genes.txt --peakfile MACSoutfile_peaks.bed --outfile MACSdirectory/piechart.pdf"
|
|
|
123 parser = ArgumentParser(usage)
|
|
|
124 parser.add_argument("--genefile", dest="genefile", help="Path to file MACS_mfold10,30_pval1e-5_genes.txt")
|
|
|
125 parser.add_argument("--peakfile", dest="peakfile", help="Path to file MACS_mfold10,30_pval1e-5_peaks.bed")
|
|
|
126 parser.add_argument("--outfile", dest="outfile", default="MACS_piechart.pdf", help="Path to pdf file where you want to store the piechart")
|
|
|
127 parser.add_argument('--MACS',action='store_true',default=False,help='Set this value if you have MACS peaks')
|
|
|
128
|
|
|
129 args=parser.parse_args()
|
|
|
130
|
|
6
|
131 print args.MACS
|
|
|
132
|
|
0
|
133 map_peaks(args.genefile, args.peakfile, args.outfile, args.MACS)
|
|
|
134
|
|
|
135
|
|
|
136 if __name__=='__main__':
|
|
|
137 main()
|
