|
0
|
1 #!/usr/bin/perl -w
|
|
|
2 #==========================================================================================
|
|
|
3 # Date:
|
|
|
4 # Title:
|
|
|
5 # Comment: Program to plot gene structure
|
|
|
6 # Input: 1. input file of Gene region annotation which format like GenePred
|
|
|
7 # 2. input file of Transcripts region annotation which format like GenePred
|
|
|
8 # 3. input file of gene snp detail info
|
|
|
9 # Output: output file of gene structure graph by html or svg formt
|
|
|
10 # Test Usage:
|
|
|
11 #========================================================================================
|
|
|
12 #use strict;
|
|
|
13 my $version=1.00;
|
|
|
14 use SVG;
|
|
|
15 use Getopt::Long;
|
|
|
16 my %opt;
|
|
|
17 GetOptions(\%opt,"g=s","l=s","chro=s","mark=s","span=s","te=s","t=s","cen=s","c=s","o=s","out=s","h");
|
|
|
18 if (!(defined $opt{g} and defined $opt{c} and defined $opt{l} and defined $opt{chro} and defined $opt{mark} and defined $opt{o}) || defined $opt{h}) {
|
|
|
19 &usage;
|
|
|
20 }
|
|
|
21 my $span=$opt{span};
|
|
|
22 #===============================Define Attribute==========================================
|
|
|
23 my %attribute=(
|
|
|
24 canvas=>{
|
|
|
25 'width'=>1500,
|
|
|
26 'height'=>1800
|
|
|
27 },
|
|
|
28 text=>{
|
|
|
29 'stroke'=>"#000000",
|
|
|
30 'fill'=>"none",
|
|
|
31 'stroke-width'=>0.5
|
|
|
32 },
|
|
|
33 line=>{
|
|
|
34 'stroke'=>"black",
|
|
|
35 'stroke-width'=>1
|
|
|
36 },
|
|
|
37 csv=>{
|
|
|
38 'stroke'=>"red",
|
|
|
39 'stroke-width'=>0.5
|
|
|
40 },
|
|
|
41 exon=>{
|
|
|
42 'stroke'=>"black",
|
|
|
43 'stroke-width'=>1
|
|
|
44 },
|
|
|
45 intron=>{
|
|
|
46 'stroke'=>"black",
|
|
|
47 'stroke-width'=>1.5
|
|
|
48 },
|
|
|
49 font=>{
|
|
|
50 'fill'=>"#000000",
|
|
|
51 'font-size'=>12,
|
|
|
52 'font-size2'=>10,
|
|
|
53 #'font-weight'=>'bold',
|
|
|
54 'font-family'=>"Arial"
|
|
|
55 #'font-family'=>"ArialNarrow-bold"
|
|
|
56 },
|
|
|
57 rect=>{
|
|
|
58 'fill'=>"lightgreen",
|
|
|
59 'stroke'=>"black",
|
|
|
60 'stroke-width'=>0.5
|
|
|
61 },
|
|
|
62 readwidth=>0.5
|
|
|
63 );
|
|
|
64 #############################s#define start coordinate and scale
|
|
|
65
|
|
|
66 my $length=$opt{"l"};#11
|
|
|
67 #my @target_e_value=qw(1.78 1.83 1.92 2.00 1.92 2.00 1.93 2.11 2.05 2.03 1.92 1.91 1.54 1.72 1.67);
|
|
|
68
|
|
|
69 my $chr=$opt{"chro"};
|
|
|
70 my $start=1;
|
|
|
71 my $XOFFSET=50;
|
|
|
72 my $YOFFSET=60;
|
|
|
73 #my $length=$end-$start+1;
|
|
|
74 my $Xscale=600/$length;#定义X轴比例尺 1:1000 x轴的坐标长度都要按照此比例尺换算
|
|
|
75 #my $high_cov=$high_cov9B1=0.5;#定义峰图最高峰
|
|
|
76 #my $Yscale=1/$high_cov;#定义Y轴比例尺 1:60 y轴的坐标长度都要按照此比例尺换算
|
|
|
77 #========================================New canvas============================
|
|
|
78 #---------------------------------------------------------------
|
|
|
79 if (defined($opt{cen})) {
|
|
|
80 open(CEN,"$opt{cen}")||die"cannot open the file $opt{cen}";
|
|
|
81 my %centromere;
|
|
|
82 while (my $aline=<CEN>) {
|
|
|
83 chomp $aline;
|
|
|
84 next if($aline=~/^\#/);
|
|
|
85 my @temp=split/\t/,$aline;
|
|
|
86 $centromere{$temp[0]}[0]=$temp[1];
|
|
|
87 $centromere{$temp[0]}[1]=$temp[2];
|
|
|
88 }
|
|
|
89 close CEN;
|
|
|
90 }
|
|
|
91
|
|
|
92 #### Starting ####
|
|
|
93 #新建画布
|
|
|
94 my $svg=SVG->new();
|
|
|
95 #画图起始点
|
|
|
96 my $canvas_start_x=$XOFFSET;
|
|
|
97 my $canvas_end_x=$XOFFSET+$length*$Xscale;#按照比例尺 画线
|
|
|
98 my $canvas_start_y=$YOFFSET;
|
|
|
99 my $canvas_end_y=$YOFFSET;
|
|
|
100 #Draw a straight line between two points start(x1,y1) and end(x2,y2).
|
|
|
101 #location attribute
|
|
|
102 $svg->line(id=>'l1',x1=>$canvas_start_x,y1=>$canvas_start_y,x2=>$canvas_end_x,y2=>$canvas_end_y,'stroke',$attribute{line}{'stroke'},'stroke-width',$attribute{line}{'stroke-width'});
|
|
|
103 $long_scale=int ($length/10);#十等分 大刻度
|
|
|
104 #大坐标刻度
|
|
|
105 for ($i=0;$i<=10;$i++) {
|
|
|
106 my $long_x_start=$XOFFSET+$long_scale*$i*$Xscale;
|
|
|
107 my $long_x_end=$long_x_start;
|
|
|
108 my $long_y_start=$YOFFSET;
|
|
|
109 my $long_y_end=$YOFFSET-5;
|
|
|
110 $svg->line('x1',$long_x_start,'y1',$long_y_start,'x2',$long_x_end,'y2',$long_y_end,'stroke',$attribute{line}{'stroke'},'stroke-width',$attribute{line}{'stroke-width'});
|
|
|
111 my $Bscale=$start+$long_scale*$i;
|
|
|
112 my $cdata=int ($Bscale/1000000);
|
|
|
113 $svg->text('x',$long_x_start,'y',$long_y_start-10,'style','fill:black;text-anchor:middle','stroke',$attribute{text}{'stroke'},'stroke-width',$attribute{text}{'stroke-width'},'font-size',12,'font-family',$attribute{font}{'font-family'},'-cdata',$cdata."M");
|
|
|
114 }
|
|
|
115
|
|
|
116 if (defined($opt{cen})) {
|
|
|
117 $svg->rect('x',$XOFFSET+$centromere{$chr2}[0]*$Xscale,'y',$YOFFSET-2,'width',($centromere{$chr2}[1]-$centromere{$chr2}[0]+1)*$Xscale,'height',5,'stroke',"black",'stroke-width',"5",'fill',"black");
|
|
|
118
|
|
|
119 $svg->rect('x',$XOFFSET+$length*$Xscale+15,'y',$YOFFSET+20,'width',10,'height',5,'stroke',"black",'stroke-width',$attribute{intron}{'stroke-width'},'fill',"black");
|
|
|
120 $svg->text('x',$XOFFSET+$length*$Xscale+35,'y',$YOFFSET+20+5,'stroke',"black",'stroke-width',$attribute{text}{'stroke-width'},'font-size',$attribute{font}{'font-size'},'font-family',$attribute{font}{'font-family'},'-cdata',"Centromere");
|
|
|
121 }
|
|
|
122
|
|
|
123 $svg->text('x',$XOFFSET+$length*$Xscale*0.42,'y',$YOFFSET-30,'stroke',"black",'stroke-width',1,'font-size',15,'font-family',$attribute{font}{'font-family'},'-cdata',$chr);
|
|
|
124 #====================================================MAIN PROGRAM========================
|
|
|
125
|
|
|
126 #===========================================gene list data================================
|
|
|
127 #open (TXT,">$opt{out}");
|
|
|
128 #open(TE,"$opt{te}")||die"cannot open the file $opt{te}";
|
|
|
129 #my %te;
|
|
|
130 #while (my $aline=<TE>) {
|
|
|
131 # chomp $aline;
|
|
|
132 # my @temp=split/\t/,$aline;
|
|
|
133 # if ($temp[2] eq "Y") {
|
|
|
134 # $te{$temp[0]}=1;
|
|
|
135 # }
|
|
|
136 #}
|
|
|
137 #close TE;
|
|
|
138 #===================================Exp gene list data =================================
|
|
|
139 #open(TARGET,"$opt{t}")||die"cannot open the file $opt{t}";
|
|
|
140 #my %target_e;
|
|
|
141 #my %target_rpkm;
|
|
|
142 #while (my $aline=<TARGET>) {
|
|
|
143 # chomp $aline;##Gene 19B1 1PA1 1LC1 29B2 2PA2 2LC2 39B3 3PA3 3LC3 93-4C PA-4C LY-4C 9343 PA43 LY43 19B1_VS_1LC1 19B1_VS_1PA1 1PA1_VS_1LC1 29B2_VS_2LC2 29B2_VS_2PA2 2PA2_VS_2LC2 39B3_VS_3LC3 39B3_VS_3PA3 3PA3_VS_3LC3 934C_LY4C 934C_PA4C PA4C_LY4C 9343_VS_LY43 9343_VS_PA43 PA43_VS_LY43 mpv_1 fold_1 mpv_2 fold_2 mpv_3 fold_3 mpv_4_B fold_4_B mpv_leaf fold_lead
|
|
|
144 # next if($aline=~/^\#/);
|
|
|
145 # my @temp=split/\t/,$aline;
|
|
|
146 # $target_rpkm{$temp[0]}="$temp[7]\t$temp[8]\t$temp[9]";
|
|
|
147 # if ($temp[7]>$target_e_value[6]||$temp[8]>$target_e_value[7]||$temp[9]>$target_e_value[8]) {
|
|
|
148 # $target_e{$temp[0]}="$temp[7]\t$temp[8]\t$temp[9]";
|
|
|
149 # }
|
|
|
150 #}
|
|
|
151 #close TARGET;
|
|
|
152 #====================================================================================
|
|
|
153 my @genelist;
|
|
|
154 #my @te_genelist;
|
|
|
155 my @target_e;
|
|
|
156 open(GENE,"$opt{g}")||die"cannot open the file $opt{g}";
|
|
|
157 while (my $aline=<GENE>) {
|
|
|
158 chomp $aline;#LOC_Os01g01280 Chr1 133291 134685 +
|
|
|
159 my @temp=split/\t/,$aline;
|
|
|
160 #my $te;
|
|
|
161 if ($temp[1]=~/^Chr(\d)$/) {
|
|
|
162 $temp[1]="Chr0$1";
|
|
|
163 }
|
|
|
164 next if($temp[1] ne $chr);
|
|
|
165 #push @genelist,[$temp[0],$temp[3],$temp[4]];
|
|
|
166 push @genelist,[$temp[0],$temp[2],$temp[3]];
|
|
|
167 # if ($te{$temp[0]}) {
|
|
|
168 # push @te_genelist,[$temp[0],$temp[3],$temp[4]];
|
|
|
169 # }
|
|
|
170 # else{
|
|
|
171 # push @genelist,[$temp[0],$temp[3],$temp[4]];
|
|
|
172 # }
|
|
|
173 # if ($target_e{$temp[0]}) {
|
|
|
174 # push @target_e,[$temp[0],$temp[3],$temp[4]];
|
|
|
175 # }
|
|
|
176
|
|
|
177 }
|
|
|
178 close GENE;
|
|
|
179 my @gene_desity;
|
|
|
180 #my @region_target_rpkm;
|
|
|
181 @genelist=sort{$a->[1] <=> $b->[1]}@genelist;
|
|
|
182 for (my $i=0;$i<@genelist ;$i++) {
|
|
|
183 # if ($genelist[$i][1]>($num1+1)*$span) {
|
|
|
184 # $num1++;
|
|
|
185 # }
|
|
|
186 # if ($genelist[$i][1]>$num1*$span&&$genelist[$i][1]<=($num1+1)*$span) {
|
|
|
187 # $gene_desity[$num1]++;
|
|
|
188 #
|
|
|
189 # }
|
|
|
190 my $start=int($genelist[$i][1]/$span);
|
|
|
191 my $end=int($genelist[$i][2]/$span);
|
|
|
192 #my @t_rpkm=split/\t/,$target_rpkm{$genelist[$i][0]};
|
|
|
193 if ($start==$end) {
|
|
|
194 $gene_desity[$start]++;
|
|
|
195 #for (my $rs=0;$rs<3 ;$rs++) {
|
|
|
196 #print TXT "$rs\t$i\t$genelist[$i][0]\t$target_rpkm{$genelist[$i][0]}\n";
|
|
|
197 #$region_target_rpkm[$start][$rs]+=$t_rpkm[$rs];
|
|
|
198 #}
|
|
|
199
|
|
|
200 #$target_rpkm_desity[$start][0]+=$temp[0];
|
|
|
201 #$target_rpkm_desity[$start][1]+=$temp[1];
|
|
|
202 #$target_rpkm_desity[$start][2]+=$temp[2];
|
|
|
203 }
|
|
|
204 else{
|
|
|
205 for (my $k=$start;$k<=$end ;$k++) {
|
|
|
206 $gene_desity[$k]++;
|
|
|
207 #for (my $rs=0;$rs<3 ;$rs++) {
|
|
|
208 #$region_target_rpkm[$k][$rs]+=$t_rpkm[$rs];
|
|
|
209 #}
|
|
|
210 #$target_rpkm_desity[$k][0]+=$temp[0];
|
|
|
211 #$target_rpkm_desity[$k][1]+=$temp[1];
|
|
|
212 #$target_rpkm_desity[$k][2]+=$temp[2];
|
|
|
213 }
|
|
|
214 }
|
|
|
215 }
|
|
|
216 #------------------------------------------region_gene_number-------------------------
|
|
|
217 my $max_gene_number=0;
|
|
|
218 my $total=0;
|
|
|
219 for (my $i=0;$i<@gene_desity ;$i++) {
|
|
|
220 if (!(defined($gene_desity[$i]))) {
|
|
|
221 $gene_desity[$i]=0;
|
|
|
222 }
|
|
|
223 if ($gene_desity[$i]>$max_gene_number) {
|
|
|
224 $max_gene_number=$gene_desity[$i];
|
|
|
225 }
|
|
|
226 #print TXT "$i\t$gene_desity[$i]\n";
|
|
|
227 $total+=$gene_desity[$i];
|
|
|
228 }
|
|
|
229 print "Gene max:$max_gene_number\ntotal:$total\n";
|
|
|
230 my $abc=@genelist;
|
|
|
231 #my $cba=@te_genelist;
|
|
|
232 print "aaaaaa:$abc\n";
|
|
|
233 #print "bbbbbb:$cba\n";
|
|
|
234 #---------------------------------------------region_gene_rpkm------------------------
|
|
|
235 #my $max_region_gene_rpkm=0;
|
|
|
236 #my @max_region_gene_rpkm=qw(0 0 0);
|
|
|
237 #my @total_region_gene_rpkm=qw(0 0 0);
|
|
|
238 #for (my $i=0;$i<@region_target_rpkm ;$i++) {
|
|
|
239 # for (my $s=0; $s<3;$s++) {
|
|
|
240 #
|
|
|
241 # if (!(defined($region_target_rpkm[$i][$s]))) {
|
|
|
242 # $region_target_rpkm[$i][$s]=0;
|
|
|
243 # }
|
|
|
244 # $total_region_gene_rpkm[$s]+=$region_target_rpkm[$i][$s];
|
|
|
245 # if ($max_region_gene_rpkm<$region_target_rpkm[$i][$s]) {
|
|
|
246 # $max_region_gene_rpkm=$region_target_rpkm[$i][$s];
|
|
|
247 # }
|
|
|
248 # if ($max_region_gene_rpkm[$s]<$region_target_rpkm[$i][$s]) {
|
|
|
249 # $max_region_gene_rpkm[$s]=$region_target_rpkm[$i][$s];
|
|
|
250 # }
|
|
|
251 # }
|
|
|
252 #}
|
|
|
253 #my @ave_region_gene_rpkm=qw(0 0 0);
|
|
|
254 #my $max_ave_rpkm=0;
|
|
|
255 #for (my $i=0;$i<3;$i++) {
|
|
|
256 # $ave_region_gene_rpkm[$i]=$total_region_gene_rpkm[$i]/($#region_target_rpkm+1);
|
|
|
257 # if ($max_ave_rpkm<$ave_region_gene_rpkm[$i]) {
|
|
|
258 # $max_ave_rpkm=$ave_region_gene_rpkm[$i];
|
|
|
259 # }
|
|
|
260 #}
|
|
|
261 #
|
|
|
262 #print "***max region gene rpkm :$max_region_gene_rpkm\n\n";
|
|
|
263 #print "@max_region_gene_rpkm\n";
|
|
|
264 #if ($max_region_gene_rpkm>10*$max_ave_rpkm) {
|
|
|
265 # $max_region_gene_rpkm=10*$max_ave_rpkm;
|
|
|
266 #}
|
|
|
267 #my @max_region_rpkm;
|
|
|
268 #----------------------------------------------------------------------------------
|
|
|
269 #my @te_gene_desity;
|
|
|
270 #@te_genelist=sort{$a->[1] <=> $b->[1]}@te_genelist;
|
|
|
271 ##my $num2=0;
|
|
|
272 #for (my $i=0;$i<@te_genelist ;$i++) {
|
|
|
273 ## if ($te_genelist[$i][1]>($num2+1)*$span) {
|
|
|
274 ## $num2=int($te_genelist[$i][1]/$span);
|
|
|
275 ## }
|
|
|
276 ## if ($te_genelist[$i][1]>$num2*$span&&$te_genelist[$i][1]<=($num2+1)*$span) {
|
|
|
277 ## $te_gene_desity[$num2]++;
|
|
|
278 ## }
|
|
|
279 # my $start=int($te_genelist[$i][1]/$span);
|
|
|
280 # my $end=int($te_genelist[$i][2]/$span);
|
|
|
281 # if ($start==$end) {
|
|
|
282 # $te_gene_desity[$start]++;
|
|
|
283 # }
|
|
|
284 # else{
|
|
|
285 # for (my $k=$start;$k<=$end ;$k++) {
|
|
|
286 # $te_gene_desity[$k]++;
|
|
|
287 # }
|
|
|
288 # }
|
|
|
289 #}
|
|
|
290 #$max_te_gene_number=0;
|
|
|
291 #$total=0;
|
|
|
292 #for (my $i=0;$i<@te_gene_desity ;$i++) {
|
|
|
293 # if (!(defined($te_gene_desity[$i]))) {
|
|
|
294 # $te_gene_desity[$i]=0;
|
|
|
295 # }
|
|
|
296 # if ($te_gene_desity[$i]>$max_te_gene_number) {
|
|
|
297 # $max_te_gene_number=$te_gene_desity[$i];
|
|
|
298 # }
|
|
|
299 # print TXT "$i\t$te_gene_desity[$i]\n";
|
|
|
300 # $total+=$te_gene_desity[$i];
|
|
|
301 #}
|
|
|
302 #
|
|
|
303 #print "TE gene max:$max_te_gene_number\ntotal:$total\n";
|
|
|
304 #-------------------------------------------------------
|
|
|
305 #my @target_e_desity;
|
|
|
306 #@target_e=sort{$a->[1] <=> $b->[1]}@target_e;
|
|
|
307 #for (my $i=0;$i<@target_e ;$i++) {
|
|
|
308 # my $start=int($target_e[$i][1]/$span);
|
|
|
309 # my $end=int($target_e[$i][2]/$span);
|
|
|
310 # if ($start==$end) {
|
|
|
311 # $target_e_desity[$start]++;
|
|
|
312 # }
|
|
|
313 # else{
|
|
|
314 # for (my $k=$start;$k<=$end ;$k++) {
|
|
|
315 # $target_e_desity[$k]++;
|
|
|
316 # }
|
|
|
317 # }
|
|
|
318 #}
|
|
|
319 #my $max_target_e_number=0;
|
|
|
320 #$total=0;
|
|
|
321 #for (my $i=0;$i<@target_e_desity ;$i++) {
|
|
|
322 # if (!(defined($target_e_desity[$i]))) {
|
|
|
323 # $target_e_desity[$i]=0;
|
|
|
324 # }
|
|
|
325 # if ($target_e_desity[$i]>$max_target_e_number) {
|
|
|
326 # $max_target_e_number=$target_e_desity[$i];
|
|
|
327 # }
|
|
|
328 # print TXT "$i\t$target_e_desity[$i]\n";
|
|
|
329 # $total+=$target_e_desity[$i];
|
|
|
330 #}
|
|
|
331 #
|
|
|
332 #print "Target max:$max_target_e_number\ntotal:$total\n";
|
|
|
333
|
|
|
334 #====================================cluster data=======================================
|
|
|
335 open(CLUSTER,"$opt{c}")||die"cannot open the file $opt{c}";
|
|
|
336 my $mark="$opt{mark}";
|
|
|
337 my @sample=split/\#/,$mark;
|
|
|
338 my $sample_num=@sample;
|
|
|
339 my @cluster;
|
|
|
340 my @cluster_density;
|
|
|
341 my @max_cluster_read=(0)x$sample_num;
|
|
|
342 while (my $aline=<CLUSTER>) {
|
|
|
343 next if($aline=~/^\#/);
|
|
|
344 chomp $aline;##ID chr strand start end 19B1
|
|
|
345 #clusterID major_length percent sample1 sample2 sample3
|
|
|
346 #Chr01:192429-192452 24nt 1.00 0.00 97222.22 0.00
|
|
|
347 my @temp=split/\t/,$aline;
|
|
|
348 my @id=split/\:/,$temp[0];
|
|
|
349 my @posit=split/\-/,$id[1];
|
|
|
350 next if($id[0] ne $chr);#Chr.01 chr04
|
|
|
351 push @cluster,[$id[0],$posit[0],$posit[1],@temp[3..3+$sample_num+1]];
|
|
|
352 #push @cluster,[$temp[0],$temp[3],$temp[4],$temp[11],$temp[12],$temp[13]];#ID start end rpkm(19B1,1PA1,1LC1);
|
|
|
353 for (my $i=0;$i<@sample ;$i++) {
|
|
|
354 if($temp[3+$i]>$max_cluster_read[$i]){
|
|
|
355 $max_cluster_read[$i]=$temp[3+$i];
|
|
|
356 }
|
|
|
357 }
|
|
|
358 }
|
|
|
359 close CLUSTER;
|
|
|
360 @cluster=sort{$a->[1] <=> $b->[1]}@cluster;
|
|
|
361 for(my $i=0;$i<$#cluster;$i++) {
|
|
|
362 my $start=int($cluster[$i][1]/$span);
|
|
|
363 my $end=int($cluster[$i][2]/$span);
|
|
|
364 if ($start==$end) {
|
|
|
365 for (my $j=0;$j<$sample_num ;$j++) {
|
|
|
366 $cluster_density[$start][$j]+=$cluster[$i][$j+$sample_num];
|
|
|
367 }
|
|
|
368
|
|
|
369 }
|
|
|
370 else{
|
|
|
371 for (my $m=$start;$m<=$end ;$m++) {
|
|
|
372 for (my $j=0;$j<$sample_num ;$j++) {
|
|
|
373 $cluster_density[$m][$j]+=$cluster[$i][$j+$sample_num];
|
|
|
374 }
|
|
|
375 }
|
|
|
376 }
|
|
|
377 }
|
|
|
378 my @max_cluster_density=(0)x$sample_num;
|
|
|
379 my @total_cluster_density=(0)x$sample_num;
|
|
|
380 my $max_cluster_density=0;
|
|
|
381 for (my $i=0;$i<@sample;$i++) {
|
|
|
382 for (my $k=0;$k<$#cluster_density ;$k++) {
|
|
|
383
|
|
|
384 if (!(defined($cluster_density[$k][$i]))) {
|
|
|
385 $cluster_density[$k][$i]=0;
|
|
|
386 }
|
|
|
387 $total_cluster_density[$i]+=$cluster_density[$k][$i];
|
|
|
388 if ($cluster_density[$k][$i]>$max_cluster_density[$i]) {
|
|
|
389 $max_cluster_density[$i]=$cluster_density[$k][$i];
|
|
|
390 }
|
|
|
391 if ($cluster_density[$k][$i]>$max_cluster_density) {
|
|
|
392 $max_cluster_density=$cluster_density[$k][$i];
|
|
|
393
|
|
|
394 }
|
|
|
395 }
|
|
|
396 }
|
|
|
397 my @ave_cluster_density=(0)x$sample_num;
|
|
|
398 my $max_ave=0;
|
|
|
399 for (my $i=0;$i<@sample;$i++) {
|
|
|
400 $ave_cluster_density[$i]=$total_cluster_density[$i]/($#cluster_density+1);
|
|
|
401 if ($max_ave<$ave_cluster_density[$i]) {
|
|
|
402 $max_ave=$ave_cluster_density[$i];
|
|
|
403 }
|
|
|
404 }
|
|
|
405
|
|
|
406 print "max cluster reads:@max_cluster_read\n";
|
|
|
407 print "max cluster density:@max_cluster_density\n";
|
|
|
408 if ($max_cluster_density>10*$max_ave) {
|
|
|
409 $max_cluster_density=10*$max_ave;
|
|
|
410 }
|
|
|
411 #===================================== plot gene desity =================================
|
|
|
412 #===row
|
|
|
413 my $Y1scale=5;
|
|
|
414 my $y1_gene_density=$YOFFSET+$max_gene_number*$Y1scale+10;
|
|
|
415 for (my $i=0;$i<@gene_desity-1 ;$i++) {
|
|
|
416 my $density_start_x=$XOFFSET+$i*$span*$Xscale;
|
|
|
417 my $density_start_y=$y1_gene_density-$gene_desity[$i]*$Y1scale;
|
|
|
418 my $density_end_x=$XOFFSET+($i+1)*$span*$Xscale;
|
|
|
419 my $density_end_y=$y1_gene_density-$gene_desity[$i+1]*$Y1scale;
|
|
|
420 my $heigth=$gene_desity[$i]/$max_gene_number;
|
|
|
421 $svg->rect('x',$density_start_x,'y',$density_start_y,'width',$span*$Xscale,'height',$y1_gene_density-$density_start_y,'stroke',"read",'stroke-width',$attribute{intron}{'stroke-width'},'fill',"green");
|
|
|
422 #$svg->line('x1',$density_start_x,'y1',$density_start_y,'x2',$density_end_x,'y2',$density_end_y,'stroke',"red",'stroke-width',$attribute{csv}{'stroke-width'});
|
|
|
423
|
|
|
424 }
|
|
|
425 $svg->line('x1',$XOFFSET,'y1',$y1_gene_density,'x2',$XOFFSET+$length*$Xscale,'y2',$y1_gene_density,'stroke',"black",'stroke-width',"black");
|
|
|
426
|
|
|
427 #====clomun
|
|
|
428 $svg->line('x1',$XOFFSET,'y1',$y1_gene_density,'x2',$XOFFSET,'y2',$y1_gene_density-$max_gene_number*$Y1scale,'stroke',"black",'stroke-width',"black");
|
|
|
429 $svg->text('x',$XOFFSET-15,'y',$y1_gene_density,'stroke',"black",'stroke-width',$attribute{text}{'stroke-width'},'font-size',$attribute{font}{'font-size2'},'font-family',$attribute{font}{'font-family'},'-cdata',"0");#min gene number
|
|
|
430 $svg->text('x',$XOFFSET-20,'y',$y1_gene_density-$max_gene_number*$Y1scale,'stroke',"black",'stroke-width',$attribute{text}{'stroke-width'},'font-size',$attribute{font}{'font-size2'},'font-family',$attribute{font}{'font-family'},'-cdata',$max_gene_number);#max gene number
|
|
|
431 #=========================================================================================
|
|
|
432 #=============================plot TE gene desity=========================================
|
|
|
433 #===row
|
|
|
434 =cut
|
|
|
435 my $Y2scale=$Y1scale;
|
|
|
436 my $y2_te_gene_density=$y1_gene_density;
|
|
|
437 for (my $i=0;$i<@te_gene_desity-1 ;$i++) {
|
|
|
438 my $te_density_start_x=$XOFFSET+$i*$span*$Xscale;
|
|
|
439 my $te_density_start_y=$y2_te_gene_density+$te_gene_desity[$i]*$Y2scale;
|
|
|
440 my $te_density_end_x=$XOFFSET+($i+1)*$span*$Xscale;
|
|
|
441 my $te_density_end_y=$y2_te_gene_density+$te_gene_desity[$i+1]*$Y2scale;
|
|
|
442 #my $te_heigth=$te_gene_desity[$i]/$max_gene_number;
|
|
|
443 $svg->rect('x',$te_density_start_x,'y',$y2_te_gene_density,'width',$span*$Xscale,'height',$te_density_start_y-$y2_te_gene_density,'stroke',"read",'stroke-width',$attribute{intron}{'stroke-width'},'fill',"green");
|
|
|
444 #$svg->line('x1',$density_start_x,'y1',$density_start_y,'x2',$density_end_x,'y2',$density_end_y,'stroke',"red",'stroke-width',$attribute{csv}{'stroke-width'});
|
|
|
445
|
|
|
446 }
|
|
|
447 #column
|
|
|
448 $svg->line('x1',$XOFFSET,'y1',$y2_te_gene_density,'x2',$XOFFSET,'y2',$y2_te_gene_density+$max_te_gene_number*$Y2scale,'stroke',"black",'stroke-width',"black");
|
|
|
449 $svg->text('x',$XOFFSET-20,'y',$y2_te_gene_density+$max_te_gene_number*$Y2scale,'stroke',"black",'stroke-width',$attribute{text}{'stroke-width'},'font-size',$attribute{font}{'font-size2'},'font-family',$attribute{font}{'font-family'},'-cdata',$max_te_gene_number);#min gene number
|
|
|
450 =cut
|
|
|
451 #=======================gene density txt==================================================
|
|
|
452 my $md=$span/1000;
|
|
|
453 #$svg->text('x',$XOFFSET-30,'y',$YOFFSET+10,'width',"698",'height',"298",'stroke',"black",'stroke-width',$attribute{text}{'stroke-width'},'font-size',$attribute{font}{'font-size2'},'font-family',$attribute{font}{'font-family'},'-cdata',"Number per \\30 kb","rotate","90",'writing-mode',"tb-rl");
|
|
|
454
|
|
|
455 $svg->rect('x',$XOFFSET+$length*$Xscale+15,'y',$y1_gene_density-$max_gene_number*$Y1scale*0.7,'width',10,'height',5,'stroke',"green",'stroke-width',$attribute{intron}{'stroke-width'},'fill',"green");
|
|
|
456 $svg->text('x',$XOFFSET+$length*$Xscale+35,'y',$y1_gene_density-$max_gene_number*$Y1scale*0.7+5,'stroke',"green",'stroke-width',$attribute{text}{'stroke-width'},'font-size',$attribute{font}{'font-size'},'font-family',$attribute{font}{'font-family'},'-cdata',"Genes");
|
|
|
457
|
|
|
458 #$svg->rect('x',$XOFFSET+$length*$Xscale+15,'y',$y1_gene_density-$max_gene_number*$Y1scale*0.7+20,'width',10,'height',5,'stroke',"green",'stroke-width',$attribute{intron}{'stroke-width'},'fill',"green");
|
|
|
459 #$svg->text('x',$XOFFSET+$length*$Xscale+35,'y',$y1_gene_density-$max_gene_number*$Y1scale*0.7+20+5,'stroke',"green",'stroke-width',$attribute{text}{'stroke-width'},'font-size',$attribute{font}{'font-size'},'font-family',$attribute{font}{'font-family'},'-cdata',"TE Genes");
|
|
|
460
|
|
|
461 $svg->text('x',$XOFFSET+$length*$Xscale+35,'y',$y1_gene_density-$max_gene_number*$Y1scale*0.7+20,'stroke',"black",'stroke-width',$attribute{text}{'stroke-width'},'font-size',$attribute{font}{'font-size'},'font-family',$attribute{font}{'font-family'},'-cdata',"genes number \/ $md kb");
|
|
|
462 #=========================================================================================
|
|
|
463 #=============================plot exp gene desity=========================================
|
|
|
464 =cut
|
|
|
465 my $y3_target_e_density=$y2_te_gene_density+$max_te_gene_number*$Y2scale+$max_target_e_number*$Y2scale+20;
|
|
|
466 #my $Y3scale=$Y1scale;
|
|
|
467 for (my $i=0;$i<@target_e_desity-1 ;$i++) {
|
|
|
468 my $target_e_density_start_x=$XOFFSET+$i*$span*$Xscale;
|
|
|
469 my $target_e_density_start_y=$y3_target_e_density-$target_e_desity[$i]*$Y2scale;
|
|
|
470 my $target_e_density_end_x=$XOFFSET+($i+1)*$span*$Xscale;
|
|
|
471 my $target_e_density_end_y=$y3_target_e_density-$target_e_desity[$i+1]*$Y2scale;
|
|
|
472 $svg->rect('x',$target_e_density_start_x,'y',$target_e_density_start_y,'width',$span*$Xscale,'height',$y3_target_e_density-$target_e_density_start_y,'stroke',"read",'stroke-width',$attribute{intron}{'stroke-width'},'fill',"red");
|
|
|
473
|
|
|
474 }
|
|
|
475 $svg->line('x1',$XOFFSET,'y1',$y3_target_e_density,'x2',$XOFFSET+$length*$Xscale,'y2',$y3_target_e_density,'stroke',"black",'stroke-width',"black");
|
|
|
476 #column
|
|
|
477 $svg->line('x1',$XOFFSET,'y1',$y3_target_e_density,'x2',$XOFFSET,'y2',$y3_target_e_density-$max_te_gene_number*$Y2scale,'stroke',"black",'stroke-width',"black");
|
|
|
478 $svg->text('x',$XOFFSET-15,'y',$y3_target_e_density,'stroke',"black",'stroke-width',$attribute{text}{'stroke-width'},'font-size',$attribute{font}{'font-size2'},'font-family',$attribute{font}{'font-family'},'-cdata',0);#max gene number
|
|
|
479 $svg->text('x',$XOFFSET-20,'y',$y3_target_e_density-$max_te_gene_number*$Y2scale+10,'stroke',"black",'stroke-width',$attribute{text}{'stroke-width'},'font-size',$attribute{font}{'font-size2'},'font-family',$attribute{font}{'font-family'},'-cdata',$max_target_e_number);#max gene number
|
|
|
480 #=========================================exp gene indensity txt==========================
|
|
|
481 $svg->rect('x',$XOFFSET+$length*$Xscale+15,'y',$y3_target_e_density-$max_target_e_number*$Y2scale*0.7,'width',10,'height',5,'stroke',"red",'stroke-width',$attribute{intron}{'stroke-width'},'fill',"red");
|
|
|
482 $svg->text('x',$XOFFSET+$length*$Xscale+35,'y',$y3_target_e_density-$max_target_e_number*$Y2scale*0.7+5,'stroke',"red",'stroke-width',$attribute{text}{'stroke-width'},'font-size',$attribute{font}{'font-size'},'font-family',$attribute{font}{'font-family'},'-cdata',"Exp Genes");
|
|
|
483 $svg->text('x',$XOFFSET+$length*$Xscale+35,'y',$y3_target_e_density-$max_target_e_number*$Y2scale*0.7+20,'stroke',"black",'stroke-width',$attribute{text}{'stroke-width'},'font-size',$attribute{font}{'font-size'},'font-family',$attribute{font}{'font-family'},'-cdata',"genes number \/ 50kb");
|
|
|
484 #calculate the different
|
|
|
485 =cut
|
|
|
486 #========================================================================================
|
|
|
487 my $Y3scale=0.04/$max_cluster_density*2500;
|
|
|
488 #my $y3_cluster_density=$y2_te_gene_density+$max_te_gene_number*5+$max_cluster_density*$Y3scale+10;
|
|
|
489 my $y3_cluster_density=$y1_gene_density+20;
|
|
|
490 my $y3_total=$y1_gene_density+10;
|
|
|
491 my @cluster_color=qw(fuchsia violet tomato);
|
|
|
492 for (my $s=0;$s<3 ;$s++) {
|
|
|
493 $y3_total=$y3_total+$max_cluster_density*$Y3scale+5;
|
|
|
494 $svg->line('x1',$XOFFSET,'y1',$y3_total,'x2',$XOFFSET+$length*$Xscale,'y2',$y3_total,'stroke',"black",'stroke-width',$attribute{csv}{'stroke-width'});
|
|
|
495 $svg->line('x1',$XOFFSET,'y1',$y3_total,'x2',$XOFFSET,'y2',$y3_total-$max_cluster_density*$Y3scale,'stroke',"black",'stroke-width',$attribute{csv}{'stroke-width'});
|
|
|
496
|
|
|
497 $svg->text('x',$XOFFSET-15,'y',$y3_total,'stroke',"black",'stroke-width',$attribute{text}{'stroke-width'},'font-size',$attribute{font}{'font-size2'},'font-family',$attribute{font}{'font-family'},'-cdata',"0");#min gene number
|
|
|
498
|
|
|
499 if ($max_cluster_density[$s]>$max_cluster_density) {
|
|
|
500 $svg->text('x',$XOFFSET-40,'y',$y3_total-$max_cluster_density*$Y3scale+10,'stroke',"black",'stroke-width',$attribute{text}{'stroke-width'},'font-size',$attribute{font}{'font-size2'},'font-family',$attribute{font}{'font-family'},'-cdata',$max_cluster_density[$s]);#max gene number
|
|
|
501 }
|
|
|
502 else{
|
|
|
503 $svg->text('x',$XOFFSET-40,'y',$y3_total-$max_cluster_density[$s]*$Y3scale+10,'stroke',"black",'stroke-width',$attribute{text}{'stroke-width'},'font-size',$attribute{font}{'font-size2'},'font-family',$attribute{font}{'font-family'},'-cdata',$max_cluster_density[$s]);#max gene number
|
|
|
504 }
|
|
|
505
|
|
|
506 for (my $i=0;$i<$#cluster_density ;$i++) {
|
|
|
507 my $cluster_density_start_x=$XOFFSET+$i*$span*$Xscale;
|
|
|
508 my $cluster_density_end_x=$XOFFSET+($i+1)*$span*$Xscale;
|
|
|
509 if ($cluster_density[$i][$s]>$max_cluster_density) {
|
|
|
510 $cluster_density[$i][$s]=$max_cluster_density;
|
|
|
511 }
|
|
|
512 if ($cluster_density[$i+1][$s]>$max_cluster_density) {
|
|
|
513 $cluster_density[$i+1][$s]=$max_cluster_density;
|
|
|
514 }
|
|
|
515 my $cluster_density_start_y=$y3_total-$cluster_density[$i][$s]*$Y3scale;
|
|
|
516 my $cluster_density_end_y=$y3_total-$cluster_density[$i+1][$s]*$Y3scale;
|
|
|
517 $svg->line('x1',$cluster_density_start_x,'y1',$cluster_density_start_y,'x2',$cluster_density_end_x,'y2',$cluster_density_end_y,'stroke',$cluster_color[$s],'stroke-width',$attribute{csv}{'stroke-width'});
|
|
|
518 }
|
|
|
519 }
|
|
|
520 for (my $s=0;$s<@sample ;$s++) {
|
|
|
521 $svg->line('x1',$XOFFSET+$length*$Xscale+10,'y1',$y3_cluster_density+($y3_total-$y3_cluster_density)*0.3+30*$s,'x2',$XOFFSET+$length*$Xscale+30,'y2',$y3_cluster_density+($y3_total-$y3_cluster_density)*0.3+30*$s,'stroke',$cluster_color[$s],'stroke-width',$attribute{csv}{'stroke-width'});
|
|
|
522 $svg->text('x',$XOFFSET+$length*$Xscale+35,'y',$y3_cluster_density+($y3_total-$y3_cluster_density)*0.3+30*$s+5,'stroke',$cluster_color[$s],'stroke-width',$attribute{text}{'stroke-width'},'font-size',$attribute{font}{'font-size'},'font-family',$attribute{font}{'font-family'},'-cdata',"Small RNAs ".$sample[$s]);
|
|
|
523 }
|
|
|
524
|
|
|
525 $svg->text('x',$XOFFSET+$length*$Xscale+35,'y',$y3_cluster_density+($y3_total-$y3_cluster_density)*0.3+30*3-5,'stroke',"black",'stroke-width',$attribute{text}{'stroke-width'},'font-size',$attribute{font}{'font-size'},'font-family',$attribute{font}{'font-family'},'-cdata',"indensity of sRNA \/ $md kb");
|
|
|
526 #===================================plot region target gene rpkm================
|
|
|
527 =cut
|
|
|
528 my $y4_region_gene_rpkm_1=$y3_total+10;
|
|
|
529 my $y4_region_gene_rpkm=$y4_region_gene_rpkm_1;
|
|
|
530 my $Y4scale=0.1/$max_region_gene_rpkm*1000;
|
|
|
531 my @cluster_color_t=qw(blue slateblue steelblue);
|
|
|
532 for (my $s=0;$s<3 ;$s++) {
|
|
|
533 $y4_region_gene_rpkm+=$max_region_gene_rpkm*$Y4scale+5;
|
|
|
534
|
|
|
535 $svg->line('x1',$XOFFSET,'y1',$y4_region_gene_rpkm,'x2',$XOFFSET+$length*$Xscale,'y2',$y4_region_gene_rpkm,'stroke',"black",'stroke-width',$attribute{csv}{'stroke-width'});
|
|
|
536
|
|
|
537 $svg->line('x1',$XOFFSET,'y1',$y4_region_gene_rpkm,'x2',$XOFFSET,'y2',$y4_region_gene_rpkm-$max_region_gene_rpkm*$Y4scale,'stroke',"black",'stroke-width',$attribute{csv}{'stroke-width'});
|
|
|
538
|
|
|
539 $svg->text('x',$XOFFSET-15,'y',$y4_region_gene_rpkm,'stroke',"black",'stroke-width',$attribute{text}{'stroke-width'},'font-size',$attribute{font}{'font-size2'},'font-family',$attribute{font}{'font-family'},'-cdata',"0");#min gene number
|
|
|
540
|
|
|
541 if ($max_region_gene_rpkm[$s]>$max_region_gene_rpkm) {
|
|
|
542 $svg->text('x',$XOFFSET-40,'y',$y4_region_gene_rpkm-$max_region_gene_rpkm*$Y4scale+10,'stroke',"black",'stroke-width',$attribute{text}{'stroke-width'},'font-size',$attribute{font}{'font-size2'},'font-family',$attribute{font}{'font-family'},'-cdata',$max_region_gene_rpkm[$s]);#max gene number
|
|
|
543 }
|
|
|
544 else{
|
|
|
545 $svg->text('x',$XOFFSET-40,'y',$y4_region_gene_rpkm-$max_region_gene_rpkm[$s]*$Y4scale+10,'stroke',"black",'stroke-width',$attribute{text}{'stroke-width'},'font-size',$attribute{font}{'font-size2'},'font-family',$attribute{font}{'font-family'},'-cdata',$max_region_gene_rpkm[$s]);#max gene number
|
|
|
546 }
|
|
|
547 for (my $i=0;$i<$#region_target_rpkm ;$i++) {
|
|
|
548 my $region_target_rpkm_start_x=$XOFFSET+$i*$span*$Xscale;
|
|
|
549 my $region_target_rpkm_end_x=$XOFFSET+($i+1)*$span*$Xscale;
|
|
|
550 if ($region_target_rpkm[$i][$s]>$max_region_gene_rpkm) {
|
|
|
551 $region_target_rpkm[$i][$s]=$max_region_gene_rpkm;
|
|
|
552 }
|
|
|
553 if ($region_target_rpkm[$i+1][$s]>$max_region_gene_rpkm) {
|
|
|
554 $region_target_rpkm[$i+1][$s]=$max_region_gene_rpkm;
|
|
|
555 }
|
|
|
556 my $region_target_rpkm_start_y=$y4_region_gene_rpkm-$region_target_rpkm[$i][$s]*$Y4scale;
|
|
|
557 my $region_target_rpkm_end_y=$y4_region_gene_rpkm-$region_target_rpkm[$i+1][$s]*$Y4scale;
|
|
|
558 $svg->line('x1',$region_target_rpkm_start_x,'y1',$region_target_rpkm_start_y,'x2',$region_target_rpkm_end_x,'y2',$region_target_rpkm_end_y,'stroke',$cluster_color_t[$s],'stroke-width',$attribute{csv}{'stroke-width'});
|
|
|
559 }
|
|
|
560
|
|
|
561 }
|
|
|
562 for (my $s=0;$s<3 ;$s++) {
|
|
|
563 $svg->line('x1',$XOFFSET+$length*$Xscale+10,'y1',$y4_region_gene_rpkm_1+($y4_region_gene_rpkm-$y4_region_gene_rpkm_1)*0.3+30*$s,'x2',$XOFFSET+$length*$Xscale+30,'y2',$y4_region_gene_rpkm_1+($y4_region_gene_rpkm-$y4_region_gene_rpkm_1)*0.3+30*$s,'stroke',$cluster_color_t[$s],'stroke-width',$attribute{csv}{'stroke-width'});
|
|
|
564 $svg->text('x',$XOFFSET+$length*$Xscale+35,'y',$y4_region_gene_rpkm_1+($y4_region_gene_rpkm-$y4_region_gene_rpkm_1)*0.3+30*$s+5,'stroke',$cluster_color_t[$s],'stroke-width',$attribute{text}{'stroke-width'},'font-size',$attribute{font}{'font-size'},'font-family',$attribute{font}{'font-family'},'-cdata',"Target Genes ".$sample[$s]);
|
|
|
565 }
|
|
|
566 $svg->text('x',$XOFFSET+$length*$Xscale+35,'y',$y4_region_gene_rpkm_1+($y4_region_gene_rpkm-$y4_region_gene_rpkm_1)*0.3+30*3-5,'stroke',"black",'stroke-width',$attribute{text}{'stroke-width'},'font-size',$attribute{font}{'font-size'},'font-family',$attribute{font}{'font-family'},'-cdata',"indensity of genes \/ 50kb");
|
|
|
567 =cut
|
|
|
568 #for (my $s=0;$s<3 ;$s++) {
|
|
|
569 # $svg->line('x1',$XOFFSET+$length*$Xscale+10,'y1',$y4_region_gene_rpkm_1+($y4_region_gene_rpkm-$y3_cluster_density)*0.3+30*$s,'x2',$XOFFSET+$length*$Xscale+30,'y2',$y4_region_gene_rpkm_1+($y4_region_gene_rpkm-$y3_cluster_density)*0.3+30*$s,'stroke',$cluster_color[$s],'stroke-width',$attribute{csv}{'stroke-width'});
|
|
|
570 # $svg->text('x',$XOFFSET+$length*$Xscale+35,'y',$y4_region_gene_rpkm_1+($y4_region_gene_rpkm-$y3_cluster_density)*0.3+30*$s+5,'stroke',$cluster_color[$s],'stroke-width',$attribute{text}{'stroke-width'},'font-size',$attribute{font}{'font-size'},'font-family',$attribute{font}{'font-family'},'-cdata',$sample[$s]);
|
|
|
571 #}
|
|
|
572 #=========================================================================================
|
|
|
573
|
|
|
574 #sub ExpG_number {
|
|
|
575 #
|
|
|
576 #}
|
|
|
577 #sub ExpCluster_number{
|
|
|
578 #
|
|
|
579 #}
|
|
|
580
|
|
|
581
|
|
|
582 #========================================================================================
|
|
|
583 open (OUT,">$opt{o}");
|
|
|
584 print OUT $svg->xmlify();
|
|
|
585
|
|
|
586 sub log2 {
|
|
|
587 my $n = shift;
|
|
|
588 return log($n)/log(2);
|
|
|
589 }
|
|
|
590
|
|
|
591 sub usage{
|
|
|
592 print <<"USAGE";
|
|
|
593 Version $version
|
|
|
594 Usage:
|
|
|
595 $0
|
|
|
596 options:
|
|
|
597 -l centromere length
|
|
|
598 -chro
|
|
|
599 -mark \#
|
|
|
600 -g input file of Gene region annotation which format like GenePred
|
|
|
601 -span
|
|
|
602 -c cluster file input
|
|
|
603 -o svg output
|
|
|
604 -h help
|
|
|
605 USAGE
|
|
|
606 exit(1);
|
|
|
607 } |
