|
0
|
1 #:::::::::::::::::::::::::g@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@j::::
|
|
|
2 #::::::::::::::::::::::;g@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@E:::
|
|
|
3 #:ttt:::::::::::::::::@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@p@;
|
|
|
4 #:t:::::::::::::::::g@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
|
|
|
5 #:t::::::::::::::;@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
|
|
|
6 #:::::::::::::z;@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
|
|
|
7 #::::::::::::i@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
|
|
|
8 #::::::::::::@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@$@@@@
|
|
|
9 #:::::::::::3@@@@@@@@@@@@@@@@@@@@@@@@@B@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
|
|
|
10 #::::::::::3@@@@@@@@@@@@@@@@@@@@@BEEESSE5EEEEBBM@@@@@@@@@@@@@@@@@@@@@@@@@@
|
|
|
11 #::::::::::3@@@@@@@@@@@@@@@@@@@@BEEEEEE35EE55E2355E5SBMB@@@@@@@@@@@@@@@@@$
|
|
|
12 #::::::::::@@@@@@@@@@@@@@@@@@@EEEE55533t3tttt::::::!!!!7755E755SBBMMM@@@MM
|
|
|
13 #::::::::::3@@@@@@@@@@@@@@@@@@EEEE2t3ttttt:::::::::::::::::::::::!7?5225EE
|
|
|
14 #::::::::::3@@@@@@@@@@@@@@@@@@EEEEE31t::::::::::::::::::::::::::::::::3E5@
|
|
|
15 #::::::::::3@@@@@@@@@@@@@@@@@@EEEEEEtt:::::::::::::::::::::::::::::::::353
|
|
|
16 #::::::::::3@@@@@@@@@@@@@@@@@@EEEEEE1ttz::::::::::::::::::::::::::::::::35
|
|
|
17 #:::::::::::@@@@@@@@@@@@@@@@@@EEEEEEEtz1::::::::::::::::::::::::::::::::t:
|
|
|
18 #:::::::::!3@@@@@@@@@@@@@@@@@@@EEEEEttt::::::::::::::::::::::::::::::::;zz
|
|
|
19 #::::::::::@@@@@@@@@@@@@@@@@@@@EEEEEttt:::::z;z:::::::::::::::::::::::::13
|
|
|
20 #::::::::::3B@@@@@@@@@@@@@@@@@@EEEEEEE3tt:czzztti;:::::::::::::::::::::::3
|
|
|
21 #::::ttt::::3@@@@@@@@@@@@@@@@EEEEE5EE25Ezt1EEEz5Etzzz;;;;:::::::::::::::::
|
|
|
22 #:::::::::::I9@@@@@@@@@@@@@@@@@@@@@@@@@@EEEEEE@@@@@@@@@@@@@@Ez;:::::::::::
|
|
|
23 #:::::::::::::E@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@Ez::::::
|
|
|
24 #::::::::::::::E@@@@@@@@@@@@@@@@@@@@@@@@@@@@@BE5EBB@@@@@@@@@@@@@@@EEE:::::
|
|
|
25 #:::::::::::::::@@@@@@@@@@@@@@@@@@@@@@@@@@@@E1::35@@@@@@@@@@ME3MMME2::::::
|
|
|
26 #:::::::::::::::?@@@@@@@@@@@@@@@@@@M@@@@@@@EE:::::3SB@@BBESEEt::::::::::::
|
|
|
27 #::::::::::::::::J$@@@@@@@B@@@@@@@@@@@@@@@@EE:::::::!35E33t:::::::::::::::
|
|
|
28 #:::::::::::::::::3@E@@@EE5EESE5EESE@@@@@@@Et::::::::::::tz:::::::::::::::
|
|
|
29 #:::::::::::::::::J@E$@EEE5133555SE@@@@@@@@Et:::::::::::::::::::::::::::::
|
|
|
30 #::::::::::::::::::E@E@EEEEtt3523EEE@@@@@@@E::::::::::::::::::::::::::::::
|
|
|
31 #:t::::::::::::::::JEE3@@@EEEEEEEEEE@@@@@@@E:::::::::t;:::::::::::::::::::
|
|
|
32 #:t:::::::::::::::::!5ES@EEEEEEEEES@@@@@@@@@E;:::;;;:3Ez::::::::::::::::::
|
|
|
33 #:t::::::::::::::::::::JE@@EEEEEEE@@@@@@@@@@@@@@@@ME!:::;:::::::::::::::::
|
|
|
34 #:tz::::::::::::::::::::JE@@@EEEE@@@@@@@@@@@@@@EE!:::::::t::::::::::::::::
|
|
|
35 #:t::::::::::::::::::::::3@@@@@@@@@@@@@@@@@@ESBE::::::::::::::::::::::::::
|
|
|
36 #:::::::::::::::::::::::::Q@@@@@@@@@@@@@@@@EE3EE;:::::zzzz::::::::::::::::
|
|
|
37 #:::::::::::::::::::::::::3@@@@@@@@@@@@@@@@@@@@@@NN@@@@@@Ez:::::::::::::::
|
|
|
38 #:zt:::::::::::::::::::::::3@@@@EE@@@@@@@@@@EEEEt::;z113E5t:::::::::::::::
|
|
|
39 #::tt:::::::::::::::::::::::3@@@E@@@@@@@@@@@@@@@@BEt::::::::::::::::t:::::
|
|
|
40 #:tt:t:::::::::::::::::::::::?S@@@@@@@@@@@BBEEE51!::::::::::::::zzzEt:::::
|
|
|
41 #::::::::::::::::::::::::::::::3Q@@@@@@@BEEEEEt:::::::::::::;zz@@@EE::::::
|
|
|
42 #::::::::::::::::::::::::::::::::75B@@@@@EEEtt;:::::::::;zz@@@@BEEEtz:::::
|
|
|
43 #::::::::::::::::::::::::::::::::::::?9@@@@@@@@@@@E2Ezg@@@@@B@@@EEEE1t::::
|
|
|
44 #:::::::::::::::::::::::::::::::::::::::3@@@@@@@@@@@@@@@@@@@E@EEEEEEEzzz::
|
|
|
45 #::::::::::::::::::::::::::::::::::::;@@@@@@@@@@@@@@@@@@@@@@@EEEEEEE5ttttt
|
|
|
46 #:::::::::::::::::::::::::::::::;g@@@@@@@@@@@@@@@@@@@@@@@@@@EEEEEEEEEEEtzt
|
|
|
47 #::::::::::::::::::::::::::::;@@@@@@@@@@@@@@@@@@@@@@@@@@E@@EEEEEEEEEEEE@@@
|
|
|
48 #::::::::::::::::::::::::::g@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@EEEE3EEEE@@@@@@@
|
|
|
49 #:::::::::::::::::::::;;g@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@EEEt33@@@@@@@@@@
|
|
|
50 #:::::::::::::::::;g@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@E@@@@@@EEEtg@@@@@@@@@@@@
|
|
|
51 #::::::::::::::;@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@EEEE@@@@@@@@@@@@@@@@@@@@@@@@
|
|
|
52 #:::::::::::::@@@@@@@@@@@@@@@@@$@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
|
|
|
53 #::::::::::;@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
|
|
|
54 #
|
|
|
55 # Copyleft ↄ⃝ 2012 Institut Curie
|
|
|
56 # Author(s): Valentina Boeva, Alban Lermine (Institut Curie) 2012
|
|
|
57 # Contact: valentina.boeva@curie.fr, alban.lermine@curie.fr
|
|
|
58 # This software is distributed under the terms of the GNU General
|
|
|
59 # Public License, either Version 2, June 1991 or Version 3, June 2007.
|
|
|
60
|
|
|
61 #!/usr/bin/perl -w
|
|
|
62 #outputs statistics for all genes in the list
|
|
|
63
|
|
|
64 #different boundaries
|
|
|
65 #no motif p-value for binding sites
|
|
|
66 #read directly No-resp/Up/down category
|
|
|
67
|
|
|
68 #all isoforms from the file with genes
|
|
|
69
|
|
|
70 #RNApolII sites on junctions
|
|
|
71
|
|
|
72 use strict;
|
|
|
73 use POSIX;
|
|
|
74
|
|
|
75 my $usage = qq{
|
|
|
76 $0
|
|
|
77
|
|
|
78 -----------------------------
|
|
|
79 mandatory parameters:
|
|
|
80
|
|
|
81 -g filename file with all genes
|
|
|
82 -rp filename file with sites of RNApolII
|
|
|
83 -k36 filename file with sites of K36me3
|
|
|
84 -tf filename file with sites of TF1
|
|
|
85
|
|
|
86 -i filename file with a table where to add colomnes
|
|
|
87 -add values which colomns to add
|
|
|
88
|
|
|
89 -----------------------------
|
|
|
90 optional parameters:
|
|
|
91
|
|
|
92 -o filename output filename (defaut "genes.output.txt")
|
|
|
93 -v verbose mode
|
|
|
94 -mir filename file with positions of miRNA
|
|
|
95 -k9 filename ile with sites of K9me3
|
|
|
96
|
|
|
97 -c_rp value cutoff for -rp
|
|
|
98 -c_k9 value cutoff for -k9
|
|
|
99 -c_k36 value cutoff for -k36
|
|
|
100
|
|
|
101 -selG filename selected genes (up-down-regulated)
|
|
|
102 -fluo filename file with fluorescence
|
|
|
103 -gc filename file with gc-islands
|
|
|
104
|
|
|
105 -long for each gene take the longest isoform
|
|
|
106
|
|
|
107
|
|
|
108 };
|
|
|
109
|
|
|
110 if(scalar(@ARGV) == 0){
|
|
|
111 print $usage;
|
|
|
112 exit(0);
|
|
|
113 }
|
|
|
114
|
|
|
115 ## mandatory arguments
|
|
|
116
|
|
|
117 my $RNApolFilename = "";
|
|
|
118 my $H3K36Me3polFilename = "";
|
|
|
119 my $H3K9Me3polFilename = "";
|
|
|
120 my $TF1Filename = "";
|
|
|
121 my $TF2Filename = "";
|
|
|
122 my $GenesFilename = "";
|
|
|
123 my $MirFilename = "";
|
|
|
124 my $TF1FilenameALL = "";
|
|
|
125 my $TF2FilenameALL = "";
|
|
|
126 my $SelectedGenesFilename = "";
|
|
|
127 my $fluoFile = "";
|
|
|
128 my $initialTable = "";
|
|
|
129 my $colomnesToAdd = "";
|
|
|
130
|
|
|
131 my $enhLeft = -30000;
|
|
|
132 my $longEnhLeft = -60000;
|
|
|
133 my $enhRight = -1500;
|
|
|
134 my $immediateDownstream = 2000;
|
|
|
135 my $K9dist = 5000;
|
|
|
136 my $kb5 = 5000;
|
|
|
137 my $INFINITY = 10000000000;
|
|
|
138 my $jonctionSize = 50;
|
|
|
139 ## optional arguments
|
|
|
140 my $outname = "genes.output.txt";
|
|
|
141 my $verbose = 0;
|
|
|
142 my $GCislands = "";
|
|
|
143
|
|
|
144 my $longest = 0;
|
|
|
145
|
|
|
146 #my $cutoff_tf1 = 0;
|
|
|
147 #my $cutoff_tf2 = 0;
|
|
|
148 my $cutoff_tf1All = 0;
|
|
|
149 my $cutoff_tf2All = 0;
|
|
|
150 my $cutoff_rp = 0;
|
|
|
151 my $cutoff_k9 = 0;
|
|
|
152 my $cutoff_k36 = 0;
|
|
|
153 my $ifTFcoord = 0;
|
|
|
154
|
|
|
155 ## parse command line arguments
|
|
|
156
|
|
|
157 while(scalar(@ARGV) > 0){
|
|
|
158 my $this_arg = shift @ARGV;
|
|
|
159 if ( $this_arg eq '-h') {print "$usage\n"; exit; }
|
|
|
160
|
|
|
161 elsif ( $this_arg eq '-selG') {$SelectedGenesFilename = shift @ARGV;}
|
|
|
162
|
|
|
163 elsif ( $this_arg eq '-g') {$GenesFilename = shift @ARGV;}
|
|
|
164 elsif ( $this_arg eq '-rp') {$RNApolFilename = shift @ARGV;}
|
|
|
165 elsif ( $this_arg eq '-k36') {$H3K36Me3polFilename = shift @ARGV;}
|
|
|
166 elsif ( $this_arg eq '-k9') {$H3K9Me3polFilename = shift @ARGV;}
|
|
|
167
|
|
|
168 elsif ( $this_arg eq '-tf') {$TF1Filename = shift @ARGV;}
|
|
|
169
|
|
|
170 elsif ( $this_arg eq '-v') {$verbose = 1;}
|
|
|
171
|
|
|
172 elsif ( $this_arg eq '-long') {$longest = 1;}
|
|
|
173
|
|
|
174
|
|
|
175 elsif ( $this_arg eq '-o') {$outname = shift @ARGV;}
|
|
|
176 elsif ( $this_arg eq '-mir') {$MirFilename = shift @ARGV;}
|
|
|
177
|
|
|
178 elsif ( $this_arg eq '-c_rp') {$cutoff_rp = shift @ARGV;}
|
|
|
179 elsif ( $this_arg eq '-c_k9') {$cutoff_k9 = shift @ARGV;}
|
|
|
180 elsif ( $this_arg eq '-c_k36') {$cutoff_k36 = shift @ARGV;}
|
|
|
181
|
|
|
182 elsif ( $this_arg eq '-fluo') {$fluoFile = shift @ARGV;}
|
|
|
183 elsif ( $this_arg eq '-gc') {$GCislands = shift @ARGV;}
|
|
|
184
|
|
|
185 elsif ( $this_arg eq '-i') {$initialTable = shift @ARGV;}
|
|
|
186 elsif ( $this_arg eq '-add') {$colomnesToAdd = shift @ARGV;}
|
|
|
187
|
|
|
188 elsif ( $this_arg eq '-lp') {$enhRight = shift @ARGV;}
|
|
|
189 elsif ( $this_arg eq '-rightp') {$immediateDownstream = shift @ARGV;}
|
|
|
190 elsif ( $this_arg eq '-enh') {$enhLeft = shift @ARGV;}
|
|
|
191 elsif ( $this_arg eq '-dg') {$kb5 = shift @ARGV;}
|
|
|
192
|
|
|
193 elsif ( $this_arg =~ m/^-/ ) { print "unknown flag: $this_arg\n";}
|
|
|
194 }
|
|
|
195
|
|
|
196
|
|
|
197 if ( $GenesFilename eq "" ){
|
|
|
198 die "you should specify a file with genes \n";
|
|
|
199 }
|
|
|
200 if(( $RNApolFilename eq "")&&($H3K36Me3polFilename eq "")&&($TF1Filename eq "")&&($H3K9Me3polFilename eq "")){
|
|
|
201 die "you should specify at least one file with peaks\n";
|
|
|
202 }
|
|
|
203
|
|
|
204
|
|
|
205 #-----------read selected genes----------------
|
|
|
206 my %selectedGenes;
|
|
|
207 my %selectedGenesFoldChange;
|
|
|
208 if ( $SelectedGenesFilename ne "") {
|
|
|
209 open (FILE, "<$SelectedGenesFilename ") or die "Cannot open file $SelectedGenesFilename !!!!: $!";
|
|
|
210 while (<FILE>) {
|
|
|
211 chomp;
|
|
|
212 my @a = split/\t/;
|
|
|
213 $selectedGenes{$a[1]} = $a[3];
|
|
|
214 $selectedGenesFoldChange{$a[1]} = $a[2];
|
|
|
215 #print "gene:$a[1],reg:$selectedGenes{$a[1]},FC:$selectedGenesFoldChange{$a[1]}\n";
|
|
|
216 }
|
|
|
217
|
|
|
218 close FILE;
|
|
|
219 print "\t\t$SelectedGenesFilename is read!\n" if ($verbose);
|
|
|
220 }
|
|
|
221
|
|
|
222 #-----------read genes with fluorescence---------
|
|
|
223 my %fluoGenes;
|
|
|
224 if ( $fluoFile ne "") {
|
|
|
225 open (FILE, "<$fluoFile ") or die "Cannot open file $fluoFile !!!!: $!";
|
|
|
226 my $gene = "";
|
|
|
227 my $med = 0;
|
|
|
228 my %h;
|
|
|
229 while (<FILE>) {
|
|
|
230 chomp;
|
|
|
231 my @a = split/\t/;
|
|
|
232
|
|
|
233 next if (scalar(@a)<5);
|
|
|
234 next if ($a[0] eq "probesets");
|
|
|
235 next unless ($a[0] =~m/\S/);
|
|
|
236 next unless ($a[4] =~m/\S/);
|
|
|
237 if ($gene ne "" && $gene ne $a[2]) {
|
|
|
238 #calcMed
|
|
|
239 $med = med(keys %h);
|
|
|
240 $fluoGenes{$gene} = $med;
|
|
|
241 $med=0;
|
|
|
242 delete @h{keys %h};
|
|
|
243 } else {
|
|
|
244 #$h{$a[4]} = 1;
|
|
|
245 }
|
|
|
246 $gene = $a[2];
|
|
|
247 $h{$a[4]} = 1;
|
|
|
248 #print "keys ", scalar(keys %h),"\t",keys %h,"\n";
|
|
|
249 }
|
|
|
250 if ($gene ne "") {
|
|
|
251 $med = med(keys %h);
|
|
|
252 $fluoGenes{$gene} = $med;
|
|
|
253 }
|
|
|
254 close FILE;
|
|
|
255 print "\t\t$fluoFile is read!\n" if ($verbose);;
|
|
|
256 }
|
|
|
257 #-----------read GC-islands----------------
|
|
|
258 my %GCislands;
|
|
|
259 if ($GCislands ne "") {
|
|
|
260 open (FILE, "<$GCislands ") or die "Cannot open file $GCislands !!!!: $!";
|
|
|
261
|
|
|
262 while (<FILE>) {
|
|
|
263 chomp;
|
|
|
264 my @a = split/\t/;
|
|
|
265 #bin chrom chromStart chromEnd name length cpgNum gcNum perCpg perGc obsExp
|
|
|
266 #107 chr1 36568608 36569851 CpG: 128 1243 128 766 20.6 61.6 1.09
|
|
|
267 my $chr = $a[1];
|
|
|
268 my $start = $a[2];
|
|
|
269 my $end = $a[3];
|
|
|
270 $GCislands{$chr}->{$start}=$end;
|
|
|
271 }
|
|
|
272 close FILE;
|
|
|
273 if ($verbose) {
|
|
|
274 print "$GCislands is read\n";
|
|
|
275 }
|
|
|
276 } elsif ($verbose) {
|
|
|
277 print "you did not specify a file with GC-islands\n";
|
|
|
278 }
|
|
|
279
|
|
|
280 #-----------read genes----------------
|
|
|
281
|
|
|
282 my %genes;
|
|
|
283
|
|
|
284 my $count = 0;
|
|
|
285
|
|
|
286 open (GENES, "<$GenesFilename") or die "Cannot open file $GenesFilename!!!!: $!";
|
|
|
287 <GENES>;
|
|
|
288 while (<GENES>) {
|
|
|
289 chomp;
|
|
|
290 if (/(chr.*)\s([+-])\s(\d+)\s(\d+)\s(\d+)\s(\d+)\s(\d+)\s(\S+)\s(\S+)\s\S+\s(\S+)/){
|
|
|
291 my $name = $10;
|
|
|
292 my $chr = $1;
|
|
|
293 my $strand = $2;
|
|
|
294 if ($strand eq '+') {
|
|
|
295 $strand = 1;
|
|
|
296 }
|
|
|
297 else {
|
|
|
298 $strand = -1;
|
|
|
299 }
|
|
|
300 my $leftPos = $3;
|
|
|
301 my $rightPos = $4;
|
|
|
302 my $cdsStart= $5;
|
|
|
303 my $cdsEnd= $6;
|
|
|
304 my $exonCount= $7;
|
|
|
305 my $exonStarts= $8;
|
|
|
306 my $exonEnds= $9;
|
|
|
307 my $ID = "$name\t$chr:$leftPos-$rightPos\t$count";
|
|
|
308 my $foldChange = 1;
|
|
|
309 my $reg = "NA";
|
|
|
310 my $fluo = "NA";
|
|
|
311 if ( $SelectedGenesFilename ne "") {
|
|
|
312 #print "$name\t";
|
|
|
313 if (exists($selectedGenes{$name})) {
|
|
|
314 $reg = $selectedGenes{$name};
|
|
|
315 $foldChange = $selectedGenesFoldChange{$name};
|
|
|
316 }
|
|
|
317 }
|
|
|
318 if ( $fluoFile ne "") {
|
|
|
319 if (exists($fluoGenes{$name})) {
|
|
|
320 $fluo = $fluoGenes{$name};
|
|
|
321 }
|
|
|
322 }
|
|
|
323 unless (exists($genes{$chr})) {
|
|
|
324 my %h;
|
|
|
325 $genes{$chr} = \%h;
|
|
|
326 }
|
|
|
327
|
|
|
328 my $RNAlength = 0;
|
|
|
329 my $skip = 0;
|
|
|
330
|
|
|
331 #print "$ID\n";
|
|
|
332 if($longest) {
|
|
|
333 $RNAlength = getRNAlength($exonStarts,$exonEnds);
|
|
|
334 for my $IDgene (keys %{$genes{$chr}}) {
|
|
|
335 my $nameGene= (split('\t', $IDgene))[0];
|
|
|
336 if ($nameGene eq $name && $RNAlength > $genes{$chr}->{$IDgene}{'RNAlength'}) {
|
|
|
337 #print "found longer isofome: $ID longer than $IDgene\n";
|
|
|
338 # print "$RNAlength > ".$genes{$chr}->{$IDgene}{'RNAlength'}."\n";
|
|
|
339
|
|
|
340 $ID=$IDgene;
|
|
|
341 } elsif ($nameGene eq $name && $RNAlength <= $genes{$chr}->{$IDgene}{'RNAlength'}) {
|
|
|
342 #print "found shorter isofome: $ID shorted than $IDgene\nwill skip it\n";
|
|
|
343 #print "$RNAlength <= ".$genes{$chr}->{$IDgene}{'RNAlength'}."\n";
|
|
|
344 $skip = 1;
|
|
|
345 }
|
|
|
346 }
|
|
|
347 }
|
|
|
348
|
|
|
349
|
|
|
350 unless ($skip) {
|
|
|
351
|
|
|
352 unless (exists($genes{$chr}->{$ID})) {
|
|
|
353 my %h1;
|
|
|
354 $genes{$chr}->{$ID} = \%h1;
|
|
|
355 $count++;
|
|
|
356 }
|
|
|
357
|
|
|
358 $genes{$chr}->{$ID}{'name'} = $name ;
|
|
|
359 $genes{$chr}->{$ID}{'left'} = $leftPos ;
|
|
|
360 $genes{$chr}->{$ID}{'right'} = $rightPos ;
|
|
|
361 $genes{$chr}->{$ID}{'cdsStart'} = $cdsStart;
|
|
|
362 $genes{$chr}->{$ID}{'cdsEnd'} = $cdsEnd;
|
|
|
363 $genes{$chr}->{$ID}{'strand'} = $strand;
|
|
|
364 $genes{$chr}->{$ID}{'exonCount'} = $exonCount;
|
|
|
365 $genes{$chr}->{$ID}{'exonStarts'} = $exonStarts;
|
|
|
366 $genes{$chr}->{$ID}{'exonEnds'} = $exonEnds;
|
|
|
367 $genes{$chr}->{$ID}{'TSS'} = ($strand == 1) ? $leftPos :$rightPos ;
|
|
|
368 $genes{$chr}->{$ID}{'TE'} = ($strand == -1) ? $leftPos :$rightPos ;
|
|
|
369 $genes{$chr}->{$ID}{'reg'} = $reg;
|
|
|
370 $genes{$chr}->{$ID}{'foldChange'} = $foldChange;
|
|
|
371 $genes{$chr}->{$ID}{'length'} = abs ($leftPos-$rightPos);
|
|
|
372 $genes{$chr}->{$ID}{'RNAlength'} = $RNAlength ;
|
|
|
373
|
|
|
374 $genes{$chr}->{$ID}{'fluo'} = $fluo;
|
|
|
375
|
|
|
376 $genes{$chr}->{$ID}{'RNApolScore'} = 0;
|
|
|
377 $genes{$chr}->{$ID}{'RNApolDist'} = $INFINITY;
|
|
|
378 $genes{$chr}->{$ID}{'RNApol_junctionScore'} = 0;
|
|
|
379 $genes{$chr}->{$ID}{'RNApol_junctionDist'} = $INFINITY;
|
|
|
380
|
|
|
381 $genes{$chr}->{$ID}{'K36score'} = 0;
|
|
|
382 $genes{$chr}->{$ID}{'K9promScore'} = 0;
|
|
|
383 $genes{$chr}->{$ID}{'K9promDist'} = $INFINITY;
|
|
|
384 $genes{$chr}->{$ID}{'K9largeScore'} = 0;
|
|
|
385 $genes{$chr}->{$ID}{'K9largeDist'} = $INFINITY;
|
|
|
386 $genes{$chr}->{$ID}{'TFpromScore'} = 0;
|
|
|
387 $genes{$chr}->{$ID}{'TFpromDist'} = $INFINITY;
|
|
|
388 $genes{$chr}->{$ID}{'TFenhScore'} = 0;
|
|
|
389 $genes{$chr}->{$ID}{'TFenhDist'} = $INFINITY;
|
|
|
390 $genes{$chr}->{$ID}{'TFintraScore'} = 0;
|
|
|
391 $genes{$chr}->{$ID}{'TFintraDist'} = $INFINITY;
|
|
|
392 $genes{$chr}->{$ID}{'TFallScore'} = 0;
|
|
|
393 $genes{$chr}->{$ID}{'TFallDist'} = $INFINITY;
|
|
|
394
|
|
|
395
|
|
|
396 $genes{$chr}->{$ID}{'TFFirstIntronAndIntraScore'} = 0;
|
|
|
397 $genes{$chr}->{$ID}{'TFFirstIntronAndIntraDist'} = $INFINITY;
|
|
|
398 $genes{$chr}->{$ID}{'TFFirstIntronScore'} = 0;
|
|
|
399 $genes{$chr}->{$ID}{'TFFirstIntronDist'} = $INFINITY;
|
|
|
400 $genes{$chr}->{$ID}{'TFintraMinusFirstIntronScore'} = 0;
|
|
|
401 $genes{$chr}->{$ID}{'TFintraMinusFirstIntronDist'} = $INFINITY;
|
|
|
402 $genes{$chr}->{$ID}{'TFImmDownScore'} = 0;
|
|
|
403 $genes{$chr}->{$ID}{'TFImmDownDist'} = $INFINITY;
|
|
|
404 $genes{$chr}->{$ID}{'TFpromSimpleScore'} = 0;
|
|
|
405 $genes{$chr}->{$ID}{'TFpromSimpleDist'} = $INFINITY;
|
|
|
406 $genes{$chr}->{$ID}{'TFenh60kbScore'} = 0;
|
|
|
407 $genes{$chr}->{$ID}{'TFenh60kbDist'} = $INFINITY;
|
|
|
408
|
|
|
409 $genes{$chr}->{$ID}{'TF_FirstExonScore'} = 0;
|
|
|
410 $genes{$chr}->{$ID}{'TF_FirstExonDist'} = $INFINITY;
|
|
|
411 $genes{$chr}->{$ID}{'TF_junctionScore'} = 0;
|
|
|
412 $genes{$chr}->{$ID}{'TF_junctionDist'} = $INFINITY;
|
|
|
413
|
|
|
414 $genes{$chr}->{$ID}{'TF_junctionAndIntraScore'} = 0;
|
|
|
415 $genes{$chr}->{$ID}{'TF_junctionAndIntraDist'} = $INFINITY;
|
|
|
416 $genes{$chr}->{$ID}{'TF_FirstExonAndIntraScore'} = 0;
|
|
|
417 $genes{$chr}->{$ID}{'TF_FirstExonAndIntraDist'} = $INFINITY;
|
|
|
418
|
|
|
419 $genes{$chr}->{$ID}{'TF_OtherExonsScore'} = 0;
|
|
|
420 $genes{$chr}->{$ID}{'TF_OtherExonsDist'} = $INFINITY;
|
|
|
421 $genes{$chr}->{$ID}{'TF_OtherExonsAndIntraScore'} = 0;
|
|
|
422 $genes{$chr}->{$ID}{'TF_OtherExonsAndIntraDist'} = $INFINITY;
|
|
|
423
|
|
|
424
|
|
|
425 $genes{$chr}->{$ID}{'TF_OtherIntronsScore'} = 0;
|
|
|
426 $genes{$chr}->{$ID}{'TF_OtherIntronsDist'} = $INFINITY;
|
|
|
427 $genes{$chr}->{$ID}{'TF_OtherIntronsAndIntraScore'} = 0;
|
|
|
428 $genes{$chr}->{$ID}{'TF_OtherIntronsAndIntraDist'} = $INFINITY;
|
|
|
429
|
|
|
430 $genes{$chr}->{$ID}{'TF5kbDownScore'} = 0;
|
|
|
431 $genes{$chr}->{$ID}{'TF5kbDownDist'} = $INFINITY;
|
|
|
432 $genes{$chr}->{$ID}{'K9enhScore'} = 0;
|
|
|
433 $genes{$chr}->{$ID}{'K9enhDist'} = $INFINITY;
|
|
|
434
|
|
|
435 ($genes{$chr}->{$ID}{'firstIntronStart'},$genes{$chr}->{$ID}{'firstIntronEnd'}) = getFirstIntron ($exonCount,$exonStarts,$exonEnds,$strand);
|
|
|
436 ($genes{$chr}->{$ID}{'firstExonStart'},$genes{$chr}->{$ID}{'firstExonEnd'}) = getFirstExon ($exonCount,$exonStarts,$exonEnds,$strand);
|
|
|
437
|
|
|
438 $genes{$chr}->{$ID}{'exonCount'} = $exonCount;
|
|
|
439 $genes{$chr}->{$ID}{'exonStarts'} = $exonStarts;
|
|
|
440 $genes{$chr}->{$ID}{'exonEnds'} = $exonEnds;
|
|
|
441
|
|
|
442
|
|
|
443 $genes{$chr}->{$ID}{'GCisland'} = 0;
|
|
|
444 if ($GCislands ne "") {
|
|
|
445 $genes{$chr}->{$ID}{'GCisland'} = checkIfGC ($genes{$chr}->{$ID}{'TSS'},$strand,2000,$GCislands{$chr});
|
|
|
446 }
|
|
|
447 }
|
|
|
448 }
|
|
|
449 }
|
|
|
450
|
|
|
451 print "Total genes (including isoforms) : $count\n" ;
|
|
|
452 close GENES;
|
|
|
453 print "\t\t$GenesFilename is read!\n" if ($verbose);;
|
|
|
454 #for my $gName (sort keys %{$genes{'chr18'}}) {
|
|
|
455
|
|
|
456 # print "$gName\t$genes{'chr18'}->{$gName}{'TSS'}\n";
|
|
|
457 #}
|
|
|
458
|
|
|
459 #-----------read file with sites miRNA, store as genes-----
|
|
|
460
|
|
|
461 if ( $MirFilename eq ""){
|
|
|
462 print "you did not specify file with miRNA\n" if ($verbose);;
|
|
|
463 }
|
|
|
464 else {
|
|
|
465 $count = 0;
|
|
|
466 open (MIR, "<$MirFilename ") or die "Cannot open file $MirFilename !!!!: $!";
|
|
|
467 #chr1 20669090 20669163 mmu-mir-206 960 +
|
|
|
468 while (<MIR>) {
|
|
|
469 chomp;
|
|
|
470 my ($name, $chr, $leftPos, $rightPos, $strand );
|
|
|
471 #1 . miRNA 20669091 20669163 . + . ACC="MI0000249"; ID="mmu-mir-206";
|
|
|
472 if (/([0-9XYM]+)\s.\smiRNA\s(\d+)\s(\d+)\s.\s([+-])\s.\sACC=.*ID=\"(.*)\"/) {
|
|
|
473 $name = $5;
|
|
|
474 $chr = $1;
|
|
|
475 $leftPos = $2;
|
|
|
476 $rightPos = $3;
|
|
|
477 $strand = $4;
|
|
|
478 }
|
|
|
479 elsif (/(.*)\s(\d+)\s(\d+)\s(.*)\s(.*)\s(.*)/){
|
|
|
480 $name = $4;
|
|
|
481 $chr = $1;
|
|
|
482 $leftPos = $2;
|
|
|
483 $rightPos = $3;
|
|
|
484 $strand = $6;
|
|
|
485 } else {
|
|
|
486 next;
|
|
|
487 }
|
|
|
488
|
|
|
489 unless ($chr =~ m/chr/) {
|
|
|
490 $chr = "chr".$chr;
|
|
|
491 }
|
|
|
492 my $ID = "$name\t$chr:$leftPos-$rightPos\t$count";
|
|
|
493
|
|
|
494 if ($strand eq '+') {
|
|
|
495 $strand = 1;
|
|
|
496 }
|
|
|
497 else {
|
|
|
498 $strand = -1;
|
|
|
499 }
|
|
|
500
|
|
|
501 unless (exists($genes{$chr})) {
|
|
|
502 my %h;
|
|
|
503 $genes{$chr} = \%h;
|
|
|
504 }
|
|
|
505 unless (exists($genes{$chr}->{$ID})) {
|
|
|
506 my %h1;
|
|
|
507 $genes{$chr}->{$ID} = \%h1;
|
|
|
508 $count++;
|
|
|
509 }
|
|
|
510 $genes{$chr}->{$ID}{'name'} = $name ;
|
|
|
511 $genes{$chr}->{$ID}{'left'} = $leftPos ;
|
|
|
512 $genes{$chr}->{$ID}{'right'} = $rightPos ;
|
|
|
513 $genes{$chr}->{$ID}{'cdsStart'} = $leftPos ;
|
|
|
514 $genes{$chr}->{$ID}{'cdsEnd'} = $rightPos ;
|
|
|
515 $genes{$chr}->{$ID}{'strand'} = $strand;
|
|
|
516 $genes{$chr}->{$ID}{'length'} = abs ($leftPos-$rightPos);
|
|
|
517 $genes{$chr}->{$ID}{'exonCount'} = 1;
|
|
|
518 $genes{$chr}->{$ID}{'exonStarts'} = $leftPos ;
|
|
|
519 $genes{$chr}->{$ID}{'exonEnds'} = $rightPos ;
|
|
|
520 $genes{$chr}->{$ID}{'TSS'} = ($strand == 1) ? $leftPos :$rightPos ;
|
|
|
521 $genes{$chr}->{$ID}{'TE'} = ($strand == -1) ? $leftPos :$rightPos ;
|
|
|
522 $genes{$chr}->{$ID}{'reg'} = "miRNA";
|
|
|
523 $genes{$chr}->{$ID}{'foldChange'} = 1;
|
|
|
524
|
|
|
525 $genes{$chr}->{$ID}{'fluo'} = "N/A";
|
|
|
526
|
|
|
527
|
|
|
528 $genes{$chr}->{$ID}{'RNApolScore'} = 0;
|
|
|
529 $genes{$chr}->{$ID}{'RNApolDist'} = $INFINITY;
|
|
|
530 $genes{$chr}->{$ID}{'RNApol_junctionScore'} = 0;
|
|
|
531 $genes{$chr}->{$ID}{'RNApol_junctionDist'} = $INFINITY;
|
|
|
532 $genes{$chr}->{$ID}{'K36score'} = 0;
|
|
|
533 $genes{$chr}->{$ID}{'K9promScore'} = 0;
|
|
|
534 $genes{$chr}->{$ID}{'K9promDist'} = $INFINITY;
|
|
|
535 $genes{$chr}->{$ID}{'K9largeScore'} = 0;
|
|
|
536 $genes{$chr}->{$ID}{'K9largeDist'} = $INFINITY;
|
|
|
537 $genes{$chr}->{$ID}{'TFpromScore'} = 0;
|
|
|
538 $genes{$chr}->{$ID}{'TFpromDist'} = $INFINITY;
|
|
|
539 $genes{$chr}->{$ID}{'TFenhScore'} = 0;
|
|
|
540 $genes{$chr}->{$ID}{'TFenhDist'} = $INFINITY;
|
|
|
541 $genes{$chr}->{$ID}{'TFintraScore'} = 0;
|
|
|
542 $genes{$chr}->{$ID}{'TFintraDist'} = $INFINITY;
|
|
|
543 $genes{$chr}->{$ID}{'TFallScore'} = 0;
|
|
|
544 $genes{$chr}->{$ID}{'TFallDist'} = $INFINITY;
|
|
|
545
|
|
|
546 $genes{$chr}->{$ID}{'TFFirstIntronAndIntraScore'} = 0;
|
|
|
547 $genes{$chr}->{$ID}{'TFFirstIntronAndIntraDist'} = $INFINITY;
|
|
|
548 $genes{$chr}->{$ID}{'TFFirstIntronScore'} = 0;
|
|
|
549 $genes{$chr}->{$ID}{'TFFirstIntronDist'} = $INFINITY;
|
|
|
550 $genes{$chr}->{$ID}{'TFintraMinusFirstIntronScore'} = 0;
|
|
|
551 $genes{$chr}->{$ID}{'TFintraMinusFirstIntronDist'} = $INFINITY;
|
|
|
552 $genes{$chr}->{$ID}{'TFImmDownScore'} = 0;
|
|
|
553 $genes{$chr}->{$ID}{'TFImmDownDist'} = $INFINITY;
|
|
|
554 $genes{$chr}->{$ID}{'TFpromSimpleScore'} = 0;
|
|
|
555 $genes{$chr}->{$ID}{'TFpromSimpleDist'} = $INFINITY;
|
|
|
556
|
|
|
557 $genes{$chr}->{$ID}{'TF_FirstExonScore'} = 0;
|
|
|
558 $genes{$chr}->{$ID}{'TF_FirstExonDist'} = $INFINITY;
|
|
|
559 $genes{$chr}->{$ID}{'TF_OtherExonsScore'} = 0;
|
|
|
560 $genes{$chr}->{$ID}{'TF_OtherExonsDist'} = $INFINITY;
|
|
|
561 $genes{$chr}->{$ID}{'TF_OtherExonsAndIntraScore'} = 0;
|
|
|
562 $genes{$chr}->{$ID}{'TF_OtherExonsAndIntraDist'} = $INFINITY;
|
|
|
563 $genes{$chr}->{$ID}{'TF_junctionScore'} = 0;
|
|
|
564 $genes{$chr}->{$ID}{'TF_junctionDist'} = $INFINITY;
|
|
|
565 $genes{$chr}->{$ID}{'TF_OtherIntronsScore'} = 0;
|
|
|
566 $genes{$chr}->{$ID}{'TF_OtherIntronsDist'} = $INFINITY;
|
|
|
567 $genes{$chr}->{$ID}{'TF_OtherIntronsAndIntraScore'} = 0;
|
|
|
568 $genes{$chr}->{$ID}{'TF_OtherIntronsAndIntraDist'} = $INFINITY;
|
|
|
569
|
|
|
570 $genes{$chr}->{$ID}{'TF_junctionAndIntraScore'} = 0;
|
|
|
571 $genes{$chr}->{$ID}{'TF_junctionAndIntraDist'} = $INFINITY;
|
|
|
572 $genes{$chr}->{$ID}{'TF_FirstExonAndIntraScore'} = 0;
|
|
|
573 $genes{$chr}->{$ID}{'TF_FirstExonAndIntraDist'} = $INFINITY;
|
|
|
574
|
|
|
575 $genes{$chr}->{$ID}{'TFenh60kbScore'} = 0;
|
|
|
576 $genes{$chr}->{$ID}{'TFenh60kbDist'} = $INFINITY;
|
|
|
577
|
|
|
578 $genes{$chr}->{$ID}{'TF5kbDownScore'} = 0;
|
|
|
579 $genes{$chr}->{$ID}{'TF5kbDownDist'} = $INFINITY;
|
|
|
580 $genes{$chr}->{$ID}{'K9enhScore'} = 0;
|
|
|
581 $genes{$chr}->{$ID}{'K9enhDist'} = $INFINITY;
|
|
|
582
|
|
|
583 ($genes{$chr}->{$ID}{'firstIntronStart'},$genes{$chr}->{$ID}{'firstIntronEnd'}) = (0,0);
|
|
|
584 ($genes{$chr}->{$ID}{'firstExonStart'},$genes{$chr}->{$ID}{'firstExonEnd'}) = (0,0);
|
|
|
585
|
|
|
586 $genes{$chr}->{$ID}{'GCisland'} = 0;
|
|
|
587
|
|
|
588 $genes{$chr}->{$ID}{'exonCount'} = 1;
|
|
|
589 $genes{$chr}->{$ID}{'exonStarts'} = $leftPos ;
|
|
|
590 $genes{$chr}->{$ID}{'exonEnds'} = $rightPos ;
|
|
|
591
|
|
|
592
|
|
|
593 }
|
|
|
594
|
|
|
595
|
|
|
596 close MIR;
|
|
|
597 print "\t\t$MirFilename is read!\n" if ($verbose) ;
|
|
|
598 print "$count miRNA\n" if ($verbose);;
|
|
|
599 }
|
|
|
600
|
|
|
601 #-----------read file with sites of TF1-----
|
|
|
602 my $numberOfAllSites = 0;
|
|
|
603
|
|
|
604 if ($TF1Filename eq "") {
|
|
|
605 print "No file with peaks of TF1!\n" if ($verbose) ;
|
|
|
606 } else {
|
|
|
607 open (FILE, "<$TF1Filename ") or die "Cannot open file $TF1Filename !!!!: $!";
|
|
|
608 $_ = <FILE>;
|
|
|
609 my $correction = 0;
|
|
|
610 my @a = split /\t/;
|
|
|
611 if ( $a[1] =~ m/chr/ ) {
|
|
|
612 $correction = 1;
|
|
|
613 }
|
|
|
614
|
|
|
615 while (<FILE>) {
|
|
|
616 chomp;
|
|
|
617
|
|
|
618 my @a = split /\t/;
|
|
|
619
|
|
|
620 my $chr = $a[0+$correction];
|
|
|
621 my $firstPos = $a[1+$correction];
|
|
|
622 my $LastPos = $a[2+$correction];
|
|
|
623 my $maxPos = $a[3+$correction];
|
|
|
624 if ($maxPos=~/\D/) {
|
|
|
625 $maxPos = int(($firstPos+$LastPos)/2);
|
|
|
626 }
|
|
|
627 my $score = $a[4+$correction];
|
|
|
628
|
|
|
629 for my $ID (keys %{$genes{$chr}}) {
|
|
|
630
|
|
|
631 my $distTSS = ($maxPos - $genes{$chr}->{$ID}{'TSS'})*$genes{$chr}->{$ID}{'strand'};
|
|
|
632 my $distTE = ($maxPos - $genes{$chr}->{$ID}{'TE'})*$genes{$chr}->{$ID}{'strand'};
|
|
|
633
|
|
|
634 if (($distTSS>= $enhLeft)&&($distTSS<$enhRight)) {
|
|
|
635 if ($genes{$chr}->{$ID}{'TFenhScore'}<$score) {
|
|
|
636 $genes{$chr}->{$ID}{'TFenhScore'}=$score;
|
|
|
637 $genes{$chr}->{$ID}{'TFenhDist'} = $distTSS;
|
|
|
638 }
|
|
|
639 } elsif (($distTSS>= $enhRight)&&($distTSS<=$immediateDownstream)) {
|
|
|
640 if ($genes{$chr}->{$ID}{'TFpromScore'}<$score) {
|
|
|
641 $genes{$chr}->{$ID}{'TFpromScore'}=$score;
|
|
|
642 $genes{$chr}->{$ID}{'TFpromDist'} = $distTSS;
|
|
|
643 }
|
|
|
644 } elsif (($distTSS >= $immediateDownstream)&&($distTE<=0)) {
|
|
|
645 if ($genes{$chr}->{$ID}{'TFintraScore'}<$score) {
|
|
|
646 $genes{$chr}->{$ID}{'TFintraScore'}=$score;
|
|
|
647 $genes{$chr}->{$ID}{'TFintraDist'} = $distTSS;
|
|
|
648 }
|
|
|
649 }
|
|
|
650 if (($distTSS>= $enhLeft)&&($distTE<=$kb5)) {
|
|
|
651 if ($genes{$chr}->{$ID}{'TFallScore'}<$score) {
|
|
|
652 $genes{$chr}->{$ID}{'TFallScore'}=$score;
|
|
|
653 $genes{$chr}->{$ID}{'TFallDist'} = $distTSS;
|
|
|
654 }
|
|
|
655 }
|
|
|
656
|
|
|
657 if (($distTSS>= 0)&&($distTSS<=$immediateDownstream)) {
|
|
|
658 if ($genes{$chr}->{$ID}{'TFImmDownScore'}<$score) {
|
|
|
659 $genes{$chr}->{$ID}{'TFImmDownScore'}=$score;
|
|
|
660 $genes{$chr}->{$ID}{'TFImmDownDist'} = $distTSS;
|
|
|
661 }
|
|
|
662 }
|
|
|
663 if (($distTSS<= 0)&&($distTSS>=$enhRight)) {
|
|
|
664 if ($genes{$chr}->{$ID}{'TFpromSimpleScore'}<$score) {
|
|
|
665 $genes{$chr}->{$ID}{'TFpromSimpleScore'}=$score;
|
|
|
666 $genes{$chr}->{$ID}{'TFpromSimpleDist'} = $distTSS;
|
|
|
667 }
|
|
|
668 }
|
|
|
669
|
|
|
670 my ($firstIntronStart,$firstIntronEnd)=($genes{$chr}->{$ID}{'firstIntronStart'},$genes{$chr}->{$ID}{'firstIntronEnd'});
|
|
|
671 ($firstIntronStart,$firstIntronEnd)= ($firstIntronEnd,$firstIntronStart) if ($firstIntronStart>$firstIntronEnd) ;
|
|
|
672
|
|
|
673 my ($firstExonStart,$firstExonEnd) = ($genes{$chr}->{$ID}{'firstExonStart'},$genes{$chr}->{$ID}{'firstExonEnd'}) ;
|
|
|
674 ($firstExonStart,$firstExonEnd)= ($firstExonEnd,$firstExonStart) if ($firstExonStart>$firstExonEnd) ;
|
|
|
675
|
|
|
676 if ($maxPos>=$firstIntronStart && $maxPos <= $firstIntronEnd) {
|
|
|
677 if ($genes{$chr}->{$ID}{'TFFirstIntronScore'}<$score) {
|
|
|
678 $genes{$chr}->{$ID}{'TFFirstIntronScore'}=$score;
|
|
|
679 $genes{$chr}->{$ID}{'TFFirstIntronDist'} = $distTSS;
|
|
|
680 }
|
|
|
681
|
|
|
682 if (($distTSS >= $immediateDownstream)&&($distTE<=0)) {
|
|
|
683 if ($genes{$chr}->{$ID}{'TFFirstIntronAndIntraScore'}<$score) {
|
|
|
684 $genes{$chr}->{$ID}{'TFFirstIntronAndIntraScore'}=$score;
|
|
|
685 $genes{$chr}->{$ID}{'TFFirstIntronAndIntraDist'} = $distTSS;
|
|
|
686 }
|
|
|
687 }
|
|
|
688 } elsif (($distTSS >= $immediateDownstream)&&($distTE<=0)) {
|
|
|
689 if ($genes{$chr}->{$ID}{'TFintraMinusFirstIntronScore'}<$score) {
|
|
|
690 $genes{$chr}->{$ID}{'TFintraMinusFirstIntronScore'}=$score;
|
|
|
691 $genes{$chr}->{$ID}{'TFintraMinusFirstIntronDist'} = $distTSS;
|
|
|
692 }
|
|
|
693 }
|
|
|
694 if (($distTSS>= $longEnhLeft)&&($distTSS<$enhRight)) {
|
|
|
695 if ($genes{$chr}->{$ID}{'TFenh60kbScore'}<$score) {
|
|
|
696 $genes{$chr}->{$ID}{'TFenh60kbScore'}=$score;
|
|
|
697 $genes{$chr}->{$ID}{'TFenh60kbDist'} = $distTSS;
|
|
|
698 }
|
|
|
699 }
|
|
|
700 if (($distTE>=0)&&($distTE<=$kb5)) {
|
|
|
701 if ($genes{$chr}->{$ID}{'TF5kbDownScore'}<$score) {
|
|
|
702 $genes{$chr}->{$ID}{'TF5kbDownScore'}=$score;
|
|
|
703 $genes{$chr}->{$ID}{'TF5kbDownDist'} = $distTSS;
|
|
|
704 }
|
|
|
705 }
|
|
|
706 if ($distTSS>=0 && $distTE<=0) {
|
|
|
707 my $typeIntra = &getTypeIntra($genes{$chr}->{$ID},$maxPos);
|
|
|
708 if ($typeIntra eq "f_exon") {
|
|
|
709 if ($genes{$chr}->{$ID}{'TF_FirstExonScore'}<$score) {
|
|
|
710 $genes{$chr}->{$ID}{'TF_FirstExonScore'}=$score;
|
|
|
711 $genes{$chr}->{$ID}{'TF_FirstExonDist'} = $distTSS;
|
|
|
712 }
|
|
|
713
|
|
|
714 if ($genes{$chr}->{$ID}{'TF_FirstExonAndIntraScore'}<$score && ($distTSS >= $immediateDownstream)&&($distTE<=0)) {
|
|
|
715 $genes{$chr}->{$ID}{'TF_FirstExonAndIntraScore'}=$score;
|
|
|
716 $genes{$chr}->{$ID}{'TF_FirstExonAndIntraDist'} = $distTSS;
|
|
|
717 }
|
|
|
718
|
|
|
719 } else {
|
|
|
720 if ($typeIntra eq "exon") {
|
|
|
721 if ($genes{$chr}->{$ID}{'TF_OtherExonsScore'}<$score) {
|
|
|
722 $genes{$chr}->{$ID}{'TF_OtherExonsScore'}=$score;
|
|
|
723 $genes{$chr}->{$ID}{'TF_OtherExonsDist'} = $distTSS;
|
|
|
724 }
|
|
|
725
|
|
|
726 if (($distTSS >= $immediateDownstream)&&($distTE<=0)) {
|
|
|
727 if ($genes{$chr}->{$ID}{'TF_OtherExonsAndIntraScore'}<$score) {
|
|
|
728 $genes{$chr}->{$ID}{'TF_OtherExonsAndIntraScore'}=$score;
|
|
|
729 $genes{$chr}->{$ID}{'TF_OtherExonsAndIntraDist'} = $distTSS;
|
|
|
730 }
|
|
|
731 }
|
|
|
732
|
|
|
733 } elsif ($typeIntra eq "intron") {
|
|
|
734 if ($genes{$chr}->{$ID}{'TF_OtherIntronsScore'}<$score) {
|
|
|
735 $genes{$chr}->{$ID}{'TF_OtherIntronsScore'}=$score;
|
|
|
736 $genes{$chr}->{$ID}{'TF_OtherIntronsDist'} = $distTSS;
|
|
|
737 }
|
|
|
738
|
|
|
739 if (($distTSS >= $immediateDownstream)&&($distTE<=0)) {
|
|
|
740 if ($genes{$chr}->{$ID}{'TF_OtherIntronsAndIntraScore'}<$score) {
|
|
|
741 $genes{$chr}->{$ID}{'TF_OtherIntronsAndIntraScore'}=$score;
|
|
|
742 $genes{$chr}->{$ID}{'TF_OtherIntronsAndIntraDist'} = $distTSS;
|
|
|
743 }
|
|
|
744 }
|
|
|
745 } elsif ($typeIntra eq "jonction") {
|
|
|
746 if ($genes{$chr}->{$ID}{'TF_junctionScore'}<$score) {
|
|
|
747 $genes{$chr}->{$ID}{'TF_junctionScore'}=$score;
|
|
|
748 $genes{$chr}->{$ID}{'TF_junctionDist'} = $distTSS;
|
|
|
749 }
|
|
|
750
|
|
|
751 if ($genes{$chr}->{$ID}{'TF_junctionAndIntraScore'}<$score && ($distTSS >= $immediateDownstream)&&($distTE<=0)) {
|
|
|
752 $genes{$chr}->{$ID}{'TF_junctionAndIntraScore'}=$score;
|
|
|
753 $genes{$chr}->{$ID}{'TF_junctionAndIntraDist'} = $distTSS;
|
|
|
754 }
|
|
|
755
|
|
|
756 }
|
|
|
757
|
|
|
758 }
|
|
|
759 }
|
|
|
760
|
|
|
761
|
|
|
762 }
|
|
|
763 $numberOfAllSites++;
|
|
|
764 }
|
|
|
765
|
|
|
766 close FILE;
|
|
|
767 print "\t$TF1Filename is read!\n" if ($verbose) ;
|
|
|
768 print "$numberOfAllSites sites\n" ;
|
|
|
769 }
|
|
|
770
|
|
|
771 #-----------read file with sites RNApolII-----
|
|
|
772 $numberOfAllSites = 0;
|
|
|
773
|
|
|
774 if ($RNApolFilename eq "") {
|
|
|
775 print "No file with peaks of RNA pol II!\n" if ($verbose) ;
|
|
|
776 } else {
|
|
|
777 open (FILE, "<$RNApolFilename ") or die "Cannot open file $RNApolFilename !!!!: $!";
|
|
|
778 $_ = <FILE>;
|
|
|
779 my @a = split /\t/;
|
|
|
780 my $correction = 0;
|
|
|
781
|
|
|
782 if ($a[0]=~m/chr/) {
|
|
|
783 $correction = -1;
|
|
|
784 }
|
|
|
785 #seek (FILE, 0, 0);
|
|
|
786 while (<FILE>) {
|
|
|
787 chomp;
|
|
|
788
|
|
|
789 my @a = split /\t/;
|
|
|
790
|
|
|
791 my $chr = $a[1+$correction];
|
|
|
792 my $firstPos = $a[2+$correction];
|
|
|
793 my $LastPos = $a[3+$correction];
|
|
|
794 my $maxPos = $a[4+$correction];
|
|
|
795 my $score = $a[5+$correction];
|
|
|
796 #print "$numberOfAllSites: $score\n";
|
|
|
797 next if ($score < $cutoff_rp);
|
|
|
798
|
|
|
799 for my $ID (keys %{$genes{$chr}}) {
|
|
|
800
|
|
|
801 my $distTSS = ($maxPos - $genes{$chr}->{$ID}{'TSS'})*$genes{$chr}->{$ID}{'strand'};
|
|
|
802 my $distTE = ($maxPos - $genes{$chr}->{$ID}{'TE'})*$genes{$chr}->{$ID}{'strand'};
|
|
|
803
|
|
|
804 if (($distTSS>= $enhRight)&&($distTSS<=$immediateDownstream)) {
|
|
|
805 if ($genes{$chr}->{$ID}{'RNApolScore'}<$score) {
|
|
|
806 $genes{$chr}->{$ID}{'RNApolScore'}=$score;
|
|
|
807 $genes{$chr}->{$ID}{'RNApolDist'} = $distTSS;
|
|
|
808 }
|
|
|
809 }
|
|
|
810 if ($distTSS>=0 && $distTE<=0) {
|
|
|
811 my $typeIntra = &getTypeIntra($genes{$chr}->{$ID},$maxPos);
|
|
|
812 if ($typeIntra eq "jonction") {
|
|
|
813 if ($genes{$chr}->{$ID}{'RNApol_junctionScore'}<$score) {
|
|
|
814 $genes{$chr}->{$ID}{'RNApol_junctionScore'}=$score;
|
|
|
815 $genes{$chr}->{$ID}{'RNApol_junctionDist'} = $distTSS;
|
|
|
816 }
|
|
|
817 }
|
|
|
818 }
|
|
|
819 }
|
|
|
820
|
|
|
821 $numberOfAllSites++;
|
|
|
822 }
|
|
|
823 close FILE;
|
|
|
824 print "\t$RNApolFilename is read!\n$numberOfAllSites sites\n" if ($verbose) ;
|
|
|
825
|
|
|
826 }
|
|
|
827
|
|
|
828 #-----------read file with sites K36me3-----
|
|
|
829 $numberOfAllSites = 0;
|
|
|
830 my @K36Score;
|
|
|
831
|
|
|
832 if ($H3K36Me3polFilename eq "") {
|
|
|
833 print "No file with peaks of H3K36me3!\n" if ($verbose) ;
|
|
|
834 } else {
|
|
|
835 open (FILE, "<$H3K36Me3polFilename ") or die "Cannot open file $H3K36Me3polFilename !!!!: $!";
|
|
|
836
|
|
|
837 $_ = <FILE>;
|
|
|
838 my @a = split /\t/;
|
|
|
839 my $correction = 0;
|
|
|
840
|
|
|
841 if ($a[0]=~m/chr/) {
|
|
|
842 $correction = -1;
|
|
|
843 }
|
|
|
844
|
|
|
845 while (<FILE>) {
|
|
|
846 chomp;
|
|
|
847
|
|
|
848 my @a = split /\t/;
|
|
|
849
|
|
|
850 my $chr = $a[1+$correction];
|
|
|
851 my $firstPos = $a[2+$correction];
|
|
|
852 my $lastPos = $a[3+$correction];
|
|
|
853 my $maxPos = $a[4+$correction];
|
|
|
854 my $score= $a[5+$correction];
|
|
|
855
|
|
|
856 for my $ID (keys %{$genes{$chr}}) {
|
|
|
857
|
|
|
858 my $distTSS = ($maxPos - $genes{$chr}->{$ID}{'TSS'})*$genes{$chr}->{$ID}{'strand'};
|
|
|
859 my $distTE = ($maxPos - $genes{$chr}->{$ID}{'TE'})*$genes{$chr}->{$ID}{'strand'};
|
|
|
860
|
|
|
861 my $scoreToadd = 0;
|
|
|
862
|
|
|
863 if (($firstPos>=$genes{$chr}->{$ID}{'left'})&&($lastPos<=$genes{$chr}->{$ID}{'right'})) {
|
|
|
864 $scoreToadd = $score/2.*($lastPos-$firstPos+1)/($genes{$chr}->{$ID}{'right'}-$genes{$chr}->{$ID}{'left'}+1);
|
|
|
865 }
|
|
|
866 if (($firstPos>=$genes{$chr}->{$ID}{'left'})&&($firstPos<$genes{$chr}->{$ID}{'right'})&&($lastPos>$genes{$chr}->{$ID}{'right'})) {
|
|
|
867 $scoreToadd = $score/2.*($genes{$chr}->{$ID}{'right'}-$firstPos+1)/($genes{$chr}->{$ID}{'right'}-$genes{$chr}->{$ID}{'left'}+1);
|
|
|
868 }
|
|
|
869 if (($firstPos<$genes{$chr}->{$ID}{'left'})&&($lastPos>$genes{$chr}->{$ID}{'left'})&&($lastPos<=$genes{$chr}->{$ID}{'right'})) {
|
|
|
870 $scoreToadd = $score/2.*($lastPos-$genes{$chr}->{$ID}{'left'}+1)/($genes{$chr}->{$ID}{'right'}-$genes{$chr}->{$ID}{'left'}+1);
|
|
|
871 }
|
|
|
872 if (($firstPos<$genes{$chr}->{$ID}{'left'})&&($lastPos>$genes{$chr}->{$ID}{'right'})) {
|
|
|
873 my $scoreToadd = $score/2.*($lastPos-$firstPos+1);
|
|
|
874 }
|
|
|
875 $genes{$chr}->{$ID}{'K36score'} += $scoreToadd;
|
|
|
876 }
|
|
|
877 $numberOfAllSites++;
|
|
|
878 }
|
|
|
879 close FILE;
|
|
|
880 print "\t$H3K36Me3polFilename is read!\n$numberOfAllSites sites\n" if ($verbose) ;
|
|
|
881 }
|
|
|
882
|
|
|
883 #-----------read file with sites H3K9me3-----
|
|
|
884 $numberOfAllSites = 0;
|
|
|
885
|
|
|
886
|
|
|
887 if ($H3K9Me3polFilename eq "") {
|
|
|
888 print "No file with peaks of H3K9me3!\n" if ($verbose) ;
|
|
|
889 } else {
|
|
|
890 open (FILE, "<$H3K9Me3polFilename ") or die "Cannot open file $H3K9Me3polFilename !!!!: $!";
|
|
|
891
|
|
|
892 $_ = <FILE>;
|
|
|
893 my @a = split /\t/;
|
|
|
894 my $correction = 0;
|
|
|
895
|
|
|
896 if ($a[0]=~m/chr/) {
|
|
|
897 $correction = -1;
|
|
|
898 }
|
|
|
899
|
|
|
900 while (<FILE>) {
|
|
|
901 chomp;
|
|
|
902 my @a = split /\t/;
|
|
|
903 my $chr = $a[1+$correction];
|
|
|
904 my $firstPos = $a[2+$correction];
|
|
|
905 my $LastPos = $a[3+$correction];
|
|
|
906 my $maxPos = $a[4+$correction];
|
|
|
907 my $score = $a[5+$correction];
|
|
|
908
|
|
|
909 next if ($score < $cutoff_k9);
|
|
|
910
|
|
|
911 for my $ID (keys %{$genes{$chr}}) {
|
|
|
912
|
|
|
913 my $distTSS = ($maxPos - $genes{$chr}->{$ID}{'TSS'})*$genes{$chr}->{$ID}{'strand'};
|
|
|
914 my $distTE = ($maxPos - $genes{$chr}->{$ID}{'TE'})*$genes{$chr}->{$ID}{'strand'};
|
|
|
915
|
|
|
916 if (($distTSS>= $enhRight)&&($distTSS<=$immediateDownstream)) {
|
|
|
917 if ($genes{$chr}->{$ID}{'K9promScore'}<$score) {
|
|
|
918 $genes{$chr}->{$ID}{'K9promScore'}=$score;
|
|
|
919 $genes{$chr}->{$ID}{'K9promDist'} = $distTSS;
|
|
|
920 }
|
|
|
921 } elsif (($distTSS >= -$K9dist)&&($distTSS<=$K9dist)) {
|
|
|
922 if ($genes{$chr}->{$ID}{'K9largeScore'}<$score) {
|
|
|
923 $genes{$chr}->{$ID}{'K9largeScore'}=$score;
|
|
|
924 $genes{$chr}->{$ID}{'K9largeDist'} = $distTSS;
|
|
|
925 }
|
|
|
926 }
|
|
|
927 if (($distTSS>= $enhLeft)&&($distTSS<$enhRight)) {
|
|
|
928 if ($genes{$chr}->{$ID}{'K9enhScore'}<$score) {
|
|
|
929 $genes{$chr}->{$ID}{'K9enhScore'}=$score;
|
|
|
930 $genes{$chr}->{$ID}{'K9enhDist'} = $distTSS;
|
|
|
931 }
|
|
|
932 }
|
|
|
933 }
|
|
|
934 $numberOfAllSites++;
|
|
|
935 }
|
|
|
936 close FILE;
|
|
|
937 print "\t$H3K9Me3polFilename is read!\n$numberOfAllSites sites\n" if ($verbose) ;
|
|
|
938 }
|
|
|
939
|
|
|
940 #-----------output all-----
|
|
|
941 #unless($initialTable eq "") {}
|
|
|
942
|
|
|
943
|
|
|
944 open (OUT , ">$outname") or die "Cannot open file $outname!!!!: $!";
|
|
|
945
|
|
|
946 print OUT "name\tchr\tstart\tend\tstrand\tReg\tfoldChange\t";
|
|
|
947
|
|
|
948 if ($GCislands ne "") {
|
|
|
949 print OUT "GC-island\t";
|
|
|
950 }
|
|
|
951
|
|
|
952 if ( $fluoFile ne "") {
|
|
|
953 print OUT "fluorescence\t";
|
|
|
954 }
|
|
|
955
|
|
|
956 if ($RNApolFilename ne "") {
|
|
|
957 print OUT "RNApolII_score\tRNApolII_distTSS\tRNApol_junctionScore\tRNApol_junctionDist\t";
|
|
|
958 }
|
|
|
959 if ($H3K36Me3polFilename ne "") {
|
|
|
960 print OUT "H3K36me3_score\t";
|
|
|
961 }
|
|
|
962 if ($H3K9Me3polFilename ne "") {
|
|
|
963 print OUT "H3K9me3_score_prom\tH3K9me3_distTSS_prom\tH3K9me3_score_large\tH3K9me3_distTSS_large\tH3K9me3_score_enh\tH3K9me3_distTSS_enh\t";
|
|
|
964 }
|
|
|
965 if ($TF1Filename ne "") {
|
|
|
966 print OUT "TF_score_Gene\tTF_distTSS_Gene\t";
|
|
|
967 print OUT "TF_score_Promoter\tTF_distTSS_Promoter\t";
|
|
|
968 print OUT "TF_score_ImmDown\tTF_distTSS_ImmDown\t";
|
|
|
969 print OUT "TF_score_PromoterORImmDown\tTF_distTSS_PromoterORImmDown\t";
|
|
|
970 print OUT "TF_score_Enhancer\tTF_distTSS_Enhancer\t";
|
|
|
971 print OUT "TF_score_Intragenic\tTF_distTSS_Intragenic\t";
|
|
|
972 print OUT "TF_score_GeneDownstream\tTF_distTSS_GeneDownstream\t";
|
|
|
973 print OUT "TF_score_FirstExon\tTF_distTSS_FirstExon\t";
|
|
|
974 print OUT "TF_score_FisrtIntron\tTF_distTSS_FisrtIntron\t";
|
|
|
975 print OUT "TF_score_FirstExonAND>$immediateDownstream\tTF_distTSS_FirstExonAND>$immediateDownstream\t";
|
|
|
976 print OUT "TF_score_FisrtIntronAND>$immediateDownstream\tTF_distTSS_FisrtIntronAND>$immediateDownstream\t";
|
|
|
977 #print OUT "TF_score_IntraMinusFisrtIntron\tTF_distTSS_IntraMinusFisrtIntron\t";
|
|
|
978
|
|
|
979 #print OUT "TF_score_enh60kb\tTF_distTSS_enh60kb\t";
|
|
|
980
|
|
|
981 print OUT "TF_score_Exons2,3,4,etc\tTF_distTSS_Exons2,3,4,etc\t";
|
|
|
982 print OUT "TF_score_Exons2,3,4,etcAND>$immediateDownstream\tTF_distTSS_Exons2,3,4,etcAND>$immediateDownstream\t";
|
|
|
983
|
|
|
984 print OUT "TF_score_Introns2,3,4,etc\tTF_distTSS_Introns2,3,4,etc\t";
|
|
|
985 print OUT "TF_score_Introns2,3,4,etcAND>$immediateDownstream\tTF_distTSS_Introns2,3,4,etcAND>$immediateDownstream\t";
|
|
|
986
|
|
|
987 print OUT "TF_score_EIjunction\tTF_distTSS_EIjunction\t";
|
|
|
988 print OUT "TF_score_EIjunctionAND>$immediateDownstream\tTF_distTSS_EIjunctionAND>$immediateDownstream";
|
|
|
989 }
|
|
|
990
|
|
|
991 print OUT "\n";
|
|
|
992
|
|
|
993 for my $chr (keys %genes) {
|
|
|
994 for my $ID (keys %{$genes{$chr}}) {
|
|
|
995 print OUT "$genes{$chr}->{$ID}{'name'}\t$chr\t$genes{$chr}->{$ID}{'left'}\t$genes{$chr}->{$ID}{'right'}\t$genes{$chr}->{$ID}{'strand'}\t$genes{$chr}->{$ID}{'reg'}\t$genes{$chr}->{$ID}{'foldChange'}\t";
|
|
|
996
|
|
|
997 if ($GCislands ne "") {
|
|
|
998 print OUT "$genes{$chr}->{$ID}{'GCisland'}\t";
|
|
|
999 }
|
|
|
1000
|
|
|
1001 if ( $fluoFile ne "") {
|
|
|
1002 print OUT "$genes{$chr}->{$ID}{'fluo'}\t";
|
|
|
1003 }
|
|
|
1004
|
|
|
1005 if ($RNApolFilename ne "") {
|
|
|
1006 print OUT "$genes{$chr}->{$ID}{'RNApolScore'}\t$genes{$chr}->{$ID}{'RNApolDist'}\t$genes{$chr}->{$ID}{'RNApol_junctionScore'}\t$genes{$chr}->{$ID}{'RNApol_junctionDist'}\t";
|
|
|
1007 }
|
|
|
1008 if ($H3K36Me3polFilename ne "") {
|
|
|
1009 print OUT "$genes{$chr}->{$ID}{'K36score'}\t";
|
|
|
1010 }
|
|
|
1011 if ($H3K9Me3polFilename ne "") {
|
|
|
1012 print OUT "$genes{$chr}->{$ID}{'K9promScore'}\t$genes{$chr}->{$ID}{'K9promDist'}\t$genes{$chr}->{$ID}{'K9largeScore'}\t$genes{$chr}->{$ID}{'K9largeDist'}\t$genes{$chr}->{$ID}{'K9enhScore'}\t$genes{$chr}->{$ID}{'K9enhDist'}\t";
|
|
|
1013 }
|
|
|
1014 if ($TF1Filename ne "") {
|
|
|
1015 print OUT "$genes{$chr}->{$ID}{'TFallScore'}\t$genes{$chr}->{$ID}{'TFallDist'}\t";
|
|
|
1016 print OUT "$genes{$chr}->{$ID}{'TFpromSimpleScore'}\t$genes{$chr}->{$ID}{'TFpromSimpleDist'}\t";
|
|
|
1017 print OUT "$genes{$chr}->{$ID}{'TFImmDownScore'}\t$genes{$chr}->{$ID}{'TFImmDownDist'}\t";
|
|
|
1018 print OUT "$genes{$chr}->{$ID}{'TFpromScore'}\t$genes{$chr}->{$ID}{'TFpromDist'}\t";
|
|
|
1019 print OUT "$genes{$chr}->{$ID}{'TFenhScore'}\t$genes{$chr}->{$ID}{'TFenhDist'}\t";
|
|
|
1020 print OUT "$genes{$chr}->{$ID}{'TFintraScore'}\t$genes{$chr}->{$ID}{'TFintraDist'}\t";
|
|
|
1021 print OUT "$genes{$chr}->{$ID}{'TF5kbDownScore'}\t$genes{$chr}->{$ID}{'TF5kbDownDist'}\t";
|
|
|
1022 print OUT "$genes{$chr}->{$ID}{'TF_FirstExonScore'}\t$genes{$chr}->{$ID}{'TF_FirstExonDist'}\t";
|
|
|
1023 print OUT "$genes{$chr}->{$ID}{'TF_FirstExonAndIntraScore'}\t$genes{$chr}->{$ID}{'TF_FirstExonAndIntraDist'}\t";
|
|
|
1024 print OUT "$genes{$chr}->{$ID}{'TFFirstIntronScore'}\t$genes{$chr}->{$ID}{'TFFirstIntronDist'}\t";
|
|
|
1025 print OUT "$genes{$chr}->{$ID}{'TFFirstIntronAndIntraScore'}\t$genes{$chr}->{$ID}{'TFFirstIntronAndIntraDist'}\t";
|
|
|
1026 #print OUT "$genes{$chr}->{$ID}{'TFintraMinusFirstIntronScore'}\t$genes{$chr}->{$ID}{'TFintraMinusFirstIntronDist'}\t";
|
|
|
1027
|
|
|
1028 #print OUT "$genes{$chr}->{$ID}{'TFenh60kbScore'}\t$genes{$chr}->{$ID}{'TFenh60kbDist'}\t";
|
|
|
1029
|
|
|
1030 print OUT "$genes{$chr}->{$ID}{'TF_OtherExonsScore'}\t$genes{$chr}->{$ID}{'TF_OtherExonsDist'}\t";
|
|
|
1031 print OUT "$genes{$chr}->{$ID}{'TF_OtherExonsAndIntraScore'}\t$genes{$chr}->{$ID}{'TF_OtherExonsAndIntraDist'}\t";
|
|
|
1032
|
|
|
1033 print OUT "$genes{$chr}->{$ID}{'TF_OtherIntronsScore'}\t$genes{$chr}->{$ID}{'TF_OtherIntronsDist'}\t";
|
|
|
1034 print OUT "$genes{$chr}->{$ID}{'TF_OtherIntronsAndIntraScore'}\t$genes{$chr}->{$ID}{'TF_OtherIntronsAndIntraDist'}\t";
|
|
|
1035
|
|
|
1036 print OUT "$genes{$chr}->{$ID}{'TF_junctionScore'}\t$genes{$chr}->{$ID}{'TF_junctionDist'}\t";
|
|
|
1037 print OUT "$genes{$chr}->{$ID}{'TF_junctionAndIntraScore'}\t$genes{$chr}->{$ID}{'TF_junctionAndIntraDist'}";
|
|
|
1038 }
|
|
|
1039 print OUT "\n";
|
|
|
1040 }
|
|
|
1041 }
|
|
|
1042
|
|
|
1043 close OUT;
|
|
|
1044
|
|
|
1045 ###################################
|
|
|
1046 sub med {
|
|
|
1047 my @arr = @_;
|
|
|
1048 my $med = 0;
|
|
|
1049 @arr = sort {$a <=> $b} @arr;
|
|
|
1050 if ((scalar(@arr)/2) =~ m/[\.\,]5/) {
|
|
|
1051 return $arr[floor(scalar(@arr)/2)];
|
|
|
1052 } else {
|
|
|
1053 return ($arr[scalar(@arr)/2]+$arr[scalar(@arr)/2-1])/2;
|
|
|
1054 }
|
|
|
1055 $med;
|
|
|
1056 }
|
|
|
1057
|
|
|
1058 sub checkIfGC {
|
|
|
1059 my ($TSS,$strand,$dist,$GCislandsChr)=@_;
|
|
|
1060 my $ifGC = 0;
|
|
|
1061 my $leftProm=$TSS-$dist;
|
|
|
1062 my $rightProm = $TSS;
|
|
|
1063 if ($strand== -1) {
|
|
|
1064 my $leftProm=$TSS;
|
|
|
1065 my $rightProm = $TSS+$dist;
|
|
|
1066 } #print "$leftProm\t"; print "$rightProm\n";
|
|
|
1067 for my $leftGC (keys %{$GCislandsChr}) {
|
|
|
1068 my $rightGC = $GCislandsChr->{$leftGC};
|
|
|
1069 if ($leftGC>=$leftProm&&$leftGC<=$rightProm || $rightGC>=$leftProm&&$rightGC<=$rightProm) {
|
|
|
1070 return "GC-island";
|
|
|
1071 }
|
|
|
1072 }
|
|
|
1073 return $ifGC ;
|
|
|
1074 }
|
|
|
1075
|
|
|
1076 sub getFirstIntron {
|
|
|
1077 my ($exonCount,$exonStarts,$exonEnds,$strand) = @_;
|
|
|
1078 my ($left,$right);
|
|
|
1079 if ($exonCount == 1) {
|
|
|
1080 return (0,0);
|
|
|
1081 }
|
|
|
1082 if ($strand == 1) {
|
|
|
1083 $left = (split ",", $exonEnds)[0];
|
|
|
1084 $right = (split (",", $exonStarts))[1];
|
|
|
1085 } else {
|
|
|
1086 $left = (split (",", $exonEnds))[$exonCount-2];
|
|
|
1087 $right = (split (",", $exonStarts))[$exonCount-1];
|
|
|
1088 }
|
|
|
1089 ($left,$right);
|
|
|
1090 }
|
|
|
1091
|
|
|
1092 sub getFirstExon {
|
|
|
1093 my ($exonCount,$exonStarts,$exonEnds,$strand) = @_;
|
|
|
1094 my ($left,$right);
|
|
|
1095 if ($exonCount == 1) {
|
|
|
1096 return (0,0);
|
|
|
1097 }
|
|
|
1098 if ($strand == 1) {
|
|
|
1099 $left = (split ",", $exonStarts)[0];
|
|
|
1100 $right = (split (",", $exonEnds))[0]-$jonctionSize;
|
|
|
1101 } else {
|
|
|
1102 $left = (split (",", $exonStarts))[$exonCount-1]+$jonctionSize;
|
|
|
1103 $right = (split (",", $exonEnds))[$exonCount-1];
|
|
|
1104 }
|
|
|
1105 ($left,$right);
|
|
|
1106 }
|
|
|
1107
|
|
|
1108
|
|
|
1109 sub getIntronExon {
|
|
|
1110 my ($pos,$exonCount,$exonStarts,$exonEnds,$strand) = @_;
|
|
|
1111 my (@left,@right);
|
|
|
1112 @left = (split ",", $exonStarts);
|
|
|
1113 @right = (split (",", $exonEnds));
|
|
|
1114
|
|
|
1115 for (my $i = 0; $i<$exonCount;$i++) {
|
|
|
1116 #print "$left[$i] <= $pos ? $pos <= $right[$i]\n";
|
|
|
1117 if (($left[$i]+$jonctionSize < $pos) && ($pos < $right[$i]-$jonctionSize)) {
|
|
|
1118 #print "URA!\n";
|
|
|
1119 return "exon";
|
|
|
1120 } elsif (($i+1<$exonCount)&&($right[$i]+$jonctionSize < $pos) && ($pos < $left[$i+1]-$jonctionSize)) {
|
|
|
1121 return "intron";
|
|
|
1122 }
|
|
|
1123 }
|
|
|
1124 return "jonction";
|
|
|
1125 }
|
|
|
1126
|
|
|
1127
|
|
|
1128 sub getTypeIntra {
|
|
|
1129
|
|
|
1130 my ($geneEntry, $pos) = @_;
|
|
|
1131 my $type;
|
|
|
1132
|
|
|
1133 if (($pos >= $geneEntry->{'firstIntronStart'})&&($pos <=$geneEntry->{'firstIntronEnd'})) {
|
|
|
1134 return "f_intron";
|
|
|
1135 }
|
|
|
1136 if (($pos >= $geneEntry->{'firstExonStart'})&&($pos <=$geneEntry->{'firstExonEnd'})) {
|
|
|
1137 return "f_exon";
|
|
|
1138 }
|
|
|
1139 $type = getIntronExon ($pos, $geneEntry->{'exonCount'},$geneEntry->{'exonStarts'},$geneEntry->{'exonEnds'},$geneEntry->{'strand'});
|
|
|
1140 return $type;
|
|
|
1141 }
|
|
|
1142
|
|
|
1143 sub getRNAlength {
|
|
|
1144 my ($exonStarts,$exonEnds) = @_;
|
|
|
1145 my (@left,@right);
|
|
|
1146 @left = (split ",", $exonStarts);
|
|
|
1147 @right = (split (",", $exonEnds));
|
|
|
1148 my $length = 0;
|
|
|
1149 for (my $i = 0; $i<scalar(@right);$i++) {
|
|
|
1150 $length+=$right[$i]-$left[$i];
|
|
|
1151 }
|
|
|
1152 #print STDERR "length = $length\n";
|
|
|
1153 return $length;
|
|
|
1154 }
|
