Mercurial > repos > big-tiandm > sirna_plant

comparison chr_plot.pl @ 0:07745c0958dd draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author big-tiandm
date Thu, 18 Sep 2014 21:40:25 -0400
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:07745c0958dd
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 }