Mercurial > repos > jbrayet > annotatepeaks
comparison createControlPeakSubSet.pl @ 5:1ec8c1d9c3c0 draft
Uploaded
author | jbrayet |
---|---|
date | Thu, 13 Aug 2015 11:02:59 -0400 |
parents | |
children |
comparison
equal
deleted
inserted
replaced
4:82b10f23f1fe | 5:1ec8c1d9c3c0 |
---|---|
1 #:t:::::::::::::::::g@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ | |
2 #:t::::::::::::::;@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ | |
3 #:::::::::::::z;@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ | |
4 #::::::::::::i@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ | |
5 #::::::::::::@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@$@@@@ | |
6 #:::::::::::3@@@@@@@@@@@@@@@@@@@@@@@@@B@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ | |
7 #::::::::::3@@@@@@@@@@@@@@@@@@@@@BEEESSE5EEEEBBM@@@@@@@@@@@@@@@@@@@@@@@@@@ | |
8 #::::::::::3@@@@@@@@@@@@@@@@@@@@BEEEEEE35EE55E2355E5SBMB@@@@@@@@@@@@@@@@@$ | |
9 #::::::::::@@@@@@@@@@@@@@@@@@@EEEE55533t3tttt::::::!!!!7755E755SBBMMM@@@MM | |
10 #::::::::::3@@@@@@@@@@@@@@@@@@EEEE2t3ttttt:::::::::::::::::::::::!7?5225EE | |
11 #::::::::::3@@@@@@@@@@@@@@@@@@EEEEE31t::::::::::::::::::::::::::::::::3E5@ | |
12 #::::::::::3@@@@@@@@@@@@@@@@@@EEEEEEtt:::::::::::::::::::::::::::::::::353 | |
13 #::::::::::3@@@@@@@@@@@@@@@@@@EEEEEE1ttz::::::::::::::::::::::::::::::::35 | |
14 #:::::::::::@@@@@@@@@@@@@@@@@@EEEEEEEtz1::::::::::::::::::::::::::::::::t: | |
15 #:::::::::!3@@@@@@@@@@@@@@@@@@@EEEEEttt::::::::::::::::::::::::::::::::;zz | |
16 #::::::::::@@@@@@@@@@@@@@@@@@@@EEEEEttt:::::z;z:::::::::::::::::::::::::13 | |
17 #::::::::::3B@@@@@@@@@@@@@@@@@@EEEEEEE3tt:czzztti;:::::::::::::::::::::::3 | |
18 #::::ttt::::3@@@@@@@@@@@@@@@@EEEEE5EE25Ezt1EEEz5Etzzz;;;;::::::::::::::::: | |
19 #:::::::::::I9@@@@@@@@@@@@@@@@@@@@@@@@@@EEEEEE@@@@@@@@@@@@@@Ez;::::::::::: | |
20 #:::::::::::::E@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@Ez:::::: | |
21 #::::::::::::::E@@@@@@@@@@@@@@@@@@@@@@@@@@@@@BE5EBB@@@@@@@@@@@@@@@EEE::::: | |
22 #:::::::::::::::@@@@@@@@@@@@@@@@@@@@@@@@@@@@E1::35@@@@@@@@@@ME3MMME2:::::: | |
23 #:::::::::::::::?@@@@@@@@@@@@@@@@@@M@@@@@@@EE:::::3SB@@BBESEEt:::::::::::: | |
24 #::::::::::::::::J$@@@@@@@B@@@@@@@@@@@@@@@@EE:::::::!35E33t::::::::::::::: | |
25 #:::::::::::::::::3@E@@@EE5EESE5EESE@@@@@@@Et::::::::::::tz::::::::::::::: | |
26 #:::::::::::::::::J@E$@EEE5133555SE@@@@@@@@Et::::::::::::::::::::::::::::: | |
27 #::::::::::::::::::E@E@EEEEtt3523EEE@@@@@@@E:::::::::::::::::::::::::::::: | |
28 #:t::::::::::::::::JEE3@@@EEEEEEEEEE@@@@@@@E:::::::::t;::::::::::::::::::: | |
29 #:t:::::::::::::::::!5ES@EEEEEEEEES@@@@@@@@@E;:::;;;:3Ez:::::::::::::::::: | |
30 #:t::::::::::::::::::::JE@@EEEEEEE@@@@@@@@@@@@@@@@ME!:::;::::::::::::::::: | |
31 #:tz::::::::::::::::::::JE@@@EEEE@@@@@@@@@@@@@@EE!:::::::t:::::::::::::::: | |
32 #:t::::::::::::::::::::::3@@@@@@@@@@@@@@@@@@ESBE:::::::::::::::::::::::::: | |
33 #:::::::::::::::::::::::::Q@@@@@@@@@@@@@@@@EE3EE;:::::zzzz:::::::::::::::: | |
34 #:::::::::::::::::::::::::3@@@@@@@@@@@@@@@@@@@@@@NN@@@@@@Ez::::::::::::::: | |
35 #:zt:::::::::::::::::::::::3@@@@EE@@@@@@@@@@EEEEt::;z113E5t::::::::::::::: | |
36 #::tt:::::::::::::::::::::::3@@@E@@@@@@@@@@@@@@@@BEt::::::::::::::::t::::: | |
37 #:tt:t:::::::::::::::::::::::?S@@@@@@@@@@@BBEEE51!::::::::::::::zzzEt::::: | |
38 #::::::::::::::::::::::::::::::3Q@@@@@@@BEEEEEt:::::::::::::;zz@@@EE:::::: | |
39 #::::::::::::::::::::::::::::::::75B@@@@@EEEtt;:::::::::;zz@@@@BEEEtz::::: | |
40 #::::::::::::::::::::::::::::::::::::?9@@@@@@@@@@@E2Ezg@@@@@B@@@EEEE1t:::: | |
41 #:::::::::::::::::::::::::::::::::::::::3@@@@@@@@@@@@@@@@@@@E@EEEEEEEzzz:: | |
42 #::::::::::::::::::::::::::::::::::::;@@@@@@@@@@@@@@@@@@@@@@@EEEEEEE5ttttt | |
43 #:::::::::::::::::::::::::::::::;g@@@@@@@@@@@@@@@@@@@@@@@@@@EEEEEEEEEEEtzt | |
44 #::::::::::::::::::::::::::::;@@@@@@@@@@@@@@@@@@@@@@@@@@E@@EEEEEEEEEEEE@@@ | |
45 #::::::::::::::::::::::::::g@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@EEEE3EEEE@@@@@@@ | |
46 #:::::::::::::::::::::;;g@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@EEEt33@@@@@@@@@@ | |
47 #:::::::::::::::::;g@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@E@@@@@@EEEtg@@@@@@@@@@@@ | |
48 #::::::::::::::;@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@EEEE@@@@@@@@@@@@@@@@@@@@@@@@ | |
49 #:::::::::::::@@@@@@@@@@@@@@@@@$@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ | |
50 #::::::::::;@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ | |
51 # | |
52 # Copyleft ↄ⃝ 2012 Institut Curie | |
53 # Author(s): Valentina Boeva, Alban Lermine (Institut Curie) 2012 | |
54 # Contact: valentina.boeva@curie.fr, alban.lermine@curie.fr | |
55 # This software is distributed under the terms of the GNU General | |
56 # Public License, either Version 2, June 1991 or Version 3, June 2007. | |
57 | |
58 #!/usr/bin/perl | |
59 | |
60 #create a control dataset with the same number of reads as in the SAMPLE (highest peaks) | |
61 | |
62 use strict; | |
63 use warnings; | |
64 use diagnostics; | |
65 | |
66 my $usage = qq{ | |
67 $0 | |
68 | |
69 ----------------------------- | |
70 mandatory parameters: | |
71 | |
72 -f CHiP_file | |
73 -c control_file | |
74 -o output file | |
75 ----------------------------- | |
76 optional parameters: | |
77 | |
78 -n number of files to create | |
79 | |
80 none | |
81 }; | |
82 | |
83 if(scalar(@ARGV) == 0){ | |
84 print $usage; | |
85 exit(0); | |
86 } | |
87 | |
88 ## mandatory arguments | |
89 | |
90 my $filename = ""; | |
91 my $output_fname = ""; | |
92 | |
93 my $controlFilename = ""; | |
94 | |
95 my $nBootstrap = 1; | |
96 | |
97 | |
98 ## optional arguments | |
99 | |
100 ## parse command line arguments | |
101 | |
102 while(scalar(@ARGV) > 0){ | |
103 my $this_arg = shift @ARGV; | |
104 if ( $this_arg eq '-h') {print "$usage\n"; exit; } | |
105 | |
106 elsif ( $this_arg eq '-f') {$filename = shift @ARGV;} | |
107 elsif ( $this_arg eq '-c') {$controlFilename = shift @ARGV;} | |
108 elsif ( $this_arg eq '-o') {$output_fname = shift @ARGV;} | |
109 elsif ( $this_arg eq '-n') {$nBootstrap = shift @ARGV;} | |
110 | |
111 elsif ( $this_arg =~ m/^-/ ) { print "unknown flag: $this_arg\n";} | |
112 } | |
113 | |
114 if ( $filename eq ""){ | |
115 die "you should specify chip file\n"; | |
116 } | |
117 if( $controlFilename eq ""){ | |
118 die "you should specify control file\n"; | |
119 } | |
120 | |
121 if( $output_fname eq ""){ | |
122 die "you should specify output filename\n"; | |
123 } | |
124 | |
125 | |
126 print "\n-----------------\n\n"; | |
127 | |
128 my %hash; | |
129 my $chipCount = 0; | |
130 my @header; | |
131 open FILE, "< $filename " || die "$filename : $!\n"; | |
132 while(<FILE>){ | |
133 $chipCount++; | |
134 } | |
135 close FILE; | |
136 #print "ChIP: $chipCount\n"; | |
137 | |
138 | |
139 my $controlCount = 0; | |
140 | |
141 open FILE, "< $controlFilename " || die "$controlFilename : $!\n"; | |
142 while(<FILE>){ | |
143 next if (/track/); | |
144 my $entry = $_; | |
145 my @fields = split(/\t/,$_); | |
146 $hash{$entry} = $fields[4]; | |
147 $controlCount++; | |
148 } | |
149 #print "controlCount : $controlCount\n"; | |
150 | |
151 close FILE; | |
152 open OUT, "> $output_fname" || die "$output_fname: $!\n"; | |
153 my $count = 0; | |
154 if ($controlCount>$chipCount) { | |
155 for my $entry (sort {$hash{$b}<=>$hash{$a}} keys %hash) { | |
156 print OUT $entry; | |
157 $count++; | |
158 if ($count >=$chipCount) { | |
159 last; | |
160 } | |
161 } | |
162 | |
163 } else { | |
164 for my $entry (keys %hash) { | |
165 print OUT $entry; | |
166 } | |
167 } | |
168 | |
169 close OUT; | |
170 | |
171 | |
172 for my $try (2..$nBootstrap) { | |
173 | |
174 open OUT, "> $output_fname$try" || die "$output_fname$try: $!\n"; | |
175 my $count = 0; | |
176 if ($controlCount>$chipCount) { | |
177 my $prob = $chipCount/$controlCount*1.1; | |
178 for my $entry (sort {$hash{$b}<=>$hash{$a}} keys %hash) { | |
179 my $yes = rand(); | |
180 if ($yes<=$prob) {$yes=1;}else {$yes=0;} | |
181 if ($yes) { | |
182 print OUT $entry ; | |
183 $count++; | |
184 } | |
185 | |
186 if ($count >=$chipCount ) { | |
187 last; | |
188 } | |
189 } | |
190 if ($count <$chipCount) { #do it again! | |
191 | |
192 for my $entry (sort {$hash{$b}<=>$hash{$a}} keys %hash) { | |
193 my $yes = rand(); | |
194 if ($yes<=$prob) {$yes=1;}else {$yes=0;} | |
195 if ($yes) { | |
196 print OUT $entry ; | |
197 $count++; | |
198 } | |
199 | |
200 if ($count >=$chipCount ) { | |
201 last; | |
202 } | |
203 } | |
204 | |
205 } | |
206 | |
207 } else { | |
208 for my $entry (keys %hash) { | |
209 print OUT $entry; | |
210 } | |
211 } | |
212 | |
213 close OUT; | |
214 | |
215 } |