Mercurial > repos > jbrayet > maketssdist_histones_1_0_docker
comparison createControlPeakSubSet.pl @ 1:635bfae8de52 draft
Uploaded
| author | jbrayet |
|---|---|
| date | Mon, 04 Jan 2016 10:21:59 -0500 |
| parents | |
| children |
comparison
equal
deleted
inserted
replaced
| 0:88a8e2392731 | 1:635bfae8de52 |
|---|---|
| 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 } |
