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 }