|
1
|
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 }
|
