Mercurial > repos > big-tiandm > sirna_plant

view chr_plot.pl @ 32:399d8dbdedd6 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author big-tiandm
date Sat, 15 Nov 2014 01:33:20 -0500
parents 07745c0958dd
children
line wrap: on
line source

#!/usr/bin/perl -w
#==========================================================================================
# Date: 
# Title: 
# Comment: Program to plot gene structure
# Input: 1. input file of Gene region annotation which format like GenePred
#        2. input file of Transcripts region annotation which format like GenePred
#        3. input file of gene snp detail info
# Output: output file of gene structure graph by html or svg formt
# Test Usage: 
#========================================================================================
#use strict;
my $version=1.00;
use SVG;
use Getopt::Long;
my %opt;
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");
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}) {
&usage;
}
my $span=$opt{span};
#===============================Define Attribute==========================================
my %attribute=(
	canvas=>{
		'width'=>1500,
		'height'=>1800
	},
	text=>{
		'stroke'=>"#000000",
		'fill'=>"none",
		'stroke-width'=>0.5
	},
	line=>{
		'stroke'=>"black",
		'stroke-width'=>1
	},
	csv=>{
		'stroke'=>"red",
		'stroke-width'=>0.5
	},
	exon=>{
		'stroke'=>"black",
		'stroke-width'=>1
	},
	intron=>{
		'stroke'=>"black",
		'stroke-width'=>1.5
	},
	font=>{
		'fill'=>"#000000",
		'font-size'=>12,
		'font-size2'=>10,
		#'font-weight'=>'bold',
		'font-family'=>"Arial"
		#'font-family'=>"ArialNarrow-bold"
	},
	rect=>{
		'fill'=>"lightgreen",
		'stroke'=>"black",
		'stroke-width'=>0.5
	},
	readwidth=>0.5
);
#############################s#define start coordinate and scale

my $length=$opt{"l"};#11
#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);

my $chr=$opt{"chro"};
my $start=1;
my $XOFFSET=50;
my $YOFFSET=60;
#my $length=$end-$start+1;
my $Xscale=600/$length;#定义X轴比例尺 1:1000 x轴的坐标长度都要按照此比例尺换算
#my $high_cov=$high_cov9B1=0.5;#定义峰图最高峰
#my $Yscale=1/$high_cov;#定义Y轴比例尺 1:60 y轴的坐标长度都要按照此比例尺换算
#========================================New canvas============================
#---------------------------------------------------------------
if (defined($opt{cen})) {
	open(CEN,"$opt{cen}")||die"cannot open the file $opt{cen}";
	my %centromere;
	while (my $aline=<CEN>) {
		chomp $aline;
		next if($aline=~/^\#/);
		my @temp=split/\t/,$aline;
		$centromere{$temp[0]}[0]=$temp[1];
		$centromere{$temp[0]}[1]=$temp[2];
	}
	close CEN;
}

####    Starting    ####
#新建画布
my $svg=SVG->new();
#画图起始点
my $canvas_start_x=$XOFFSET;
my $canvas_end_x=$XOFFSET+$length*$Xscale;#按照比例尺 画线
my $canvas_start_y=$YOFFSET;
my $canvas_end_y=$YOFFSET;
#Draw a straight line between two points start(x1,y1) and end(x2,y2).
#location attribute
$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'});
$long_scale=int ($length/10);#十等分 大刻度
#大坐标刻度
for ($i=0;$i<=10;$i++) {
	my $long_x_start=$XOFFSET+$long_scale*$i*$Xscale;
	my $long_x_end=$long_x_start;
	my $long_y_start=$YOFFSET;
	my $long_y_end=$YOFFSET-5;
	$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'});
	my $Bscale=$start+$long_scale*$i;
	my $cdata=int ($Bscale/1000000);
	$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");
}

if (defined($opt{cen})) {
	$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");

	$svg->rect('x',$XOFFSET+$length*$Xscale+15,'y',$YOFFSET+20,'width',10,'height',5,'stroke',"black",'stroke-width',$attribute{intron}{'stroke-width'},'fill',"black");
	$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");
}

$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);
#====================================================MAIN PROGRAM========================

#===========================================gene list data================================
#open (TXT,">$opt{out}");
#open(TE,"$opt{te}")||die"cannot open the file $opt{te}";
#my %te;
#while (my $aline=<TE>) {
#	chomp $aline;
#	my @temp=split/\t/,$aline;
#	if ($temp[2] eq "Y") {
#		$te{$temp[0]}=1;
#	}
#}
#close TE;
#===================================Exp gene list data =================================
#open(TARGET,"$opt{t}")||die"cannot open the file $opt{t}";
#my %target_e;
#my %target_rpkm;
#while (my $aline=<TARGET>) {
#	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
#	next if($aline=~/^\#/);
#	my @temp=split/\t/,$aline;
#	$target_rpkm{$temp[0]}="$temp[7]\t$temp[8]\t$temp[9]";
#	if ($temp[7]>$target_e_value[6]||$temp[8]>$target_e_value[7]||$temp[9]>$target_e_value[8]) {
#		$target_e{$temp[0]}="$temp[7]\t$temp[8]\t$temp[9]";
#	}
#}
#close TARGET;
#====================================================================================
my @genelist;
#my @te_genelist;
my @target_e;
open(GENE,"$opt{g}")||die"cannot open the file $opt{g}";
while (my $aline=<GENE>) {
	chomp $aline;#LOC_Os01g01280  Chr1    133291  134685  +
	my @temp=split/\t/,$aline;
	#my $te;
	if ($temp[1]=~/^Chr(\d)$/) {
		$temp[1]="Chr0$1";
	}
	next if($temp[1] ne $chr);
	#push @genelist,[$temp[0],$temp[3],$temp[4]];
	push @genelist,[$temp[0],$temp[2],$temp[3]];
#	if ($te{$temp[0]}) {
#		push @te_genelist,[$temp[0],$temp[3],$temp[4]];
#	}
#	else{
#		push @genelist,[$temp[0],$temp[3],$temp[4]];
#		}
#	if ($target_e{$temp[0]}) {
#		push @target_e,[$temp[0],$temp[3],$temp[4]];
#	}

}
close GENE;
my @gene_desity;
#my @region_target_rpkm;
@genelist=sort{$a->[1] <=> $b->[1]}@genelist;
for (my $i=0;$i<@genelist ;$i++) {
#	if ($genelist[$i][1]>($num1+1)*$span) {
#			$num1++;
#		}
#	if ($genelist[$i][1]>$num1*$span&&$genelist[$i][1]<=($num1+1)*$span) {
#		$gene_desity[$num1]++;
#		
#	}
	my $start=int($genelist[$i][1]/$span);
	my $end=int($genelist[$i][2]/$span);
	#my @t_rpkm=split/\t/,$target_rpkm{$genelist[$i][0]};
	if ($start==$end) {
		$gene_desity[$start]++;
		#for (my $rs=0;$rs<3 ;$rs++) {
			#print TXT "$rs\t$i\t$genelist[$i][0]\t$target_rpkm{$genelist[$i][0]}\n";
			#$region_target_rpkm[$start][$rs]+=$t_rpkm[$rs];
		#}
		
		#$target_rpkm_desity[$start][0]+=$temp[0];
		#$target_rpkm_desity[$start][1]+=$temp[1];
		#$target_rpkm_desity[$start][2]+=$temp[2];
	}
	else{
		for (my $k=$start;$k<=$end ;$k++) {
			$gene_desity[$k]++;
			#for (my $rs=0;$rs<3 ;$rs++) {
				#$region_target_rpkm[$k][$rs]+=$t_rpkm[$rs];
			#}
			#$target_rpkm_desity[$k][0]+=$temp[0];
			#$target_rpkm_desity[$k][1]+=$temp[1];
			#$target_rpkm_desity[$k][2]+=$temp[2];
		}
	}
}
#------------------------------------------region_gene_number-------------------------
my $max_gene_number=0;
my $total=0;
for (my $i=0;$i<@gene_desity ;$i++) {
	if (!(defined($gene_desity[$i]))) {
		$gene_desity[$i]=0;
	}
	if ($gene_desity[$i]>$max_gene_number) {
		$max_gene_number=$gene_desity[$i];
	}
	#print TXT "$i\t$gene_desity[$i]\n";
	$total+=$gene_desity[$i];
}
print "Gene max:$max_gene_number\ntotal:$total\n";
my $abc=@genelist;
#my $cba=@te_genelist;
print "aaaaaa:$abc\n";
#print "bbbbbb:$cba\n";
#---------------------------------------------region_gene_rpkm------------------------
#my $max_region_gene_rpkm=0;
#my @max_region_gene_rpkm=qw(0 0 0);
#my @total_region_gene_rpkm=qw(0 0 0);
#for (my $i=0;$i<@region_target_rpkm ;$i++) {
#	for (my $s=0; $s<3;$s++) {
#		
#		if (!(defined($region_target_rpkm[$i][$s]))) {
#			$region_target_rpkm[$i][$s]=0;
#		}
#		$total_region_gene_rpkm[$s]+=$region_target_rpkm[$i][$s];
#		if ($max_region_gene_rpkm<$region_target_rpkm[$i][$s]) {
#			$max_region_gene_rpkm=$region_target_rpkm[$i][$s];
#		}
#		if ($max_region_gene_rpkm[$s]<$region_target_rpkm[$i][$s]) {
#			$max_region_gene_rpkm[$s]=$region_target_rpkm[$i][$s];
#		}
#	}
#}
#my @ave_region_gene_rpkm=qw(0 0 0);
#my $max_ave_rpkm=0;
#for (my $i=0;$i<3;$i++) {
#	$ave_region_gene_rpkm[$i]=$total_region_gene_rpkm[$i]/($#region_target_rpkm+1);
#	if ($max_ave_rpkm<$ave_region_gene_rpkm[$i]) {
#		$max_ave_rpkm=$ave_region_gene_rpkm[$i];
#	}
#}
#
#print "***max region gene rpkm :$max_region_gene_rpkm\n\n";
#print "@max_region_gene_rpkm\n";
#if ($max_region_gene_rpkm>10*$max_ave_rpkm) {
#	$max_region_gene_rpkm=10*$max_ave_rpkm;
#}
#my @max_region_rpkm;
#----------------------------------------------------------------------------------
#my @te_gene_desity;
#@te_genelist=sort{$a->[1] <=> $b->[1]}@te_genelist;
##my $num2=0;
#for (my $i=0;$i<@te_genelist ;$i++) {
##	if ($te_genelist[$i][1]>($num2+1)*$span) {
##		$num2=int($te_genelist[$i][1]/$span);
##	}
##	if ($te_genelist[$i][1]>$num2*$span&&$te_genelist[$i][1]<=($num2+1)*$span) {
##		$te_gene_desity[$num2]++;
##	}	
#	my $start=int($te_genelist[$i][1]/$span);
#	my $end=int($te_genelist[$i][2]/$span);
#	if ($start==$end) {
#		$te_gene_desity[$start]++;
#	}
#	else{
#		for (my $k=$start;$k<=$end ;$k++) {
#			$te_gene_desity[$k]++;
#		}
#	}
#}
#$max_te_gene_number=0;
#$total=0;
#for (my $i=0;$i<@te_gene_desity ;$i++) {
#	if (!(defined($te_gene_desity[$i]))) {
#		$te_gene_desity[$i]=0;
#	}
#	if ($te_gene_desity[$i]>$max_te_gene_number) {
#		$max_te_gene_number=$te_gene_desity[$i];
#	}
#	print TXT "$i\t$te_gene_desity[$i]\n";
#	$total+=$te_gene_desity[$i];
#}
#
#print "TE gene max:$max_te_gene_number\ntotal:$total\n";
#-------------------------------------------------------
#my @target_e_desity;
#@target_e=sort{$a->[1] <=> $b->[1]}@target_e;
#for (my $i=0;$i<@target_e ;$i++) {	
#	my $start=int($target_e[$i][1]/$span);
#	my $end=int($target_e[$i][2]/$span);
#	if ($start==$end) {
#		$target_e_desity[$start]++;
#	}
#	else{
#		for (my $k=$start;$k<=$end ;$k++) {
#			$target_e_desity[$k]++;
#		}
#	}
#}
#my $max_target_e_number=0;
#$total=0;
#for (my $i=0;$i<@target_e_desity ;$i++) {
#	if (!(defined($target_e_desity[$i]))) {
#		$target_e_desity[$i]=0;
#	}
#	if ($target_e_desity[$i]>$max_target_e_number) {
#		$max_target_e_number=$target_e_desity[$i];
#	}
#	print TXT "$i\t$target_e_desity[$i]\n";
#	$total+=$target_e_desity[$i];
#}
#
#print "Target max:$max_target_e_number\ntotal:$total\n";

#====================================cluster data=======================================
open(CLUSTER,"$opt{c}")||die"cannot open the file $opt{c}";
my $mark="$opt{mark}";
my @sample=split/\#/,$mark;
my $sample_num=@sample;
my @cluster;
my @cluster_density;
my @max_cluster_read=(0)x$sample_num;
while (my $aline=<CLUSTER>) {
	next if($aline=~/^\#/);
	chomp $aline;##ID	chr	strand	start	end	19B1
	#clusterID	major_length	 percent	sample1	sample2	sample3
	#Chr01:192429-192452	24nt	1.00	0.00	97222.22	0.00
	my @temp=split/\t/,$aline;
	my @id=split/\:/,$temp[0];
	my @posit=split/\-/,$id[1];
	next if($id[0] ne $chr);#Chr.01 chr04
	push @cluster,[$id[0],$posit[0],$posit[1],@temp[3..3+$sample_num+1]];
	#push @cluster,[$temp[0],$temp[3],$temp[4],$temp[11],$temp[12],$temp[13]];#ID	start	end	rpkm(19B1,1PA1,1LC1);
	for (my $i=0;$i<@sample ;$i++) {
		if($temp[3+$i]>$max_cluster_read[$i]){
			$max_cluster_read[$i]=$temp[3+$i];
		}
	}
}
close CLUSTER;
@cluster=sort{$a->[1] <=> $b->[1]}@cluster;
for(my $i=0;$i<$#cluster;$i++) {
	my $start=int($cluster[$i][1]/$span);
	my $end=int($cluster[$i][2]/$span);
	if ($start==$end) {
		for (my $j=0;$j<$sample_num ;$j++) {
			$cluster_density[$start][$j]+=$cluster[$i][$j+$sample_num];
		}
		
	}
	else{
		for (my $m=$start;$m<=$end ;$m++) {
			for (my $j=0;$j<$sample_num ;$j++) {
				$cluster_density[$m][$j]+=$cluster[$i][$j+$sample_num];
			}
		}
	}
}
my @max_cluster_density=(0)x$sample_num;
my @total_cluster_density=(0)x$sample_num;
my $max_cluster_density=0;
for (my $i=0;$i<@sample;$i++) {
	for (my $k=0;$k<$#cluster_density ;$k++) {
		
		if (!(defined($cluster_density[$k][$i]))) {
			$cluster_density[$k][$i]=0;
		}
		$total_cluster_density[$i]+=$cluster_density[$k][$i];
		if ($cluster_density[$k][$i]>$max_cluster_density[$i]) {
			$max_cluster_density[$i]=$cluster_density[$k][$i];
		}
		if ($cluster_density[$k][$i]>$max_cluster_density) {
			$max_cluster_density=$cluster_density[$k][$i];

		}
	}	
}
my @ave_cluster_density=(0)x$sample_num;
my $max_ave=0;
for (my $i=0;$i<@sample;$i++) {
	$ave_cluster_density[$i]=$total_cluster_density[$i]/($#cluster_density+1);
	if ($max_ave<$ave_cluster_density[$i]) {
		$max_ave=$ave_cluster_density[$i];
	}
}

print "max cluster reads:@max_cluster_read\n";
print "max cluster density:@max_cluster_density\n";
if ($max_cluster_density>10*$max_ave) {
	$max_cluster_density=10*$max_ave;
}
#===================================== plot gene desity =================================
#===row
my $Y1scale=5;
my $y1_gene_density=$YOFFSET+$max_gene_number*$Y1scale+10;
for (my $i=0;$i<@gene_desity-1 ;$i++) {
	my $density_start_x=$XOFFSET+$i*$span*$Xscale;
	my $density_start_y=$y1_gene_density-$gene_desity[$i]*$Y1scale;
	my $density_end_x=$XOFFSET+($i+1)*$span*$Xscale;
	my $density_end_y=$y1_gene_density-$gene_desity[$i+1]*$Y1scale;
	my $heigth=$gene_desity[$i]/$max_gene_number;
	$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");
	#$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'});
	
}
$svg->line('x1',$XOFFSET,'y1',$y1_gene_density,'x2',$XOFFSET+$length*$Xscale,'y2',$y1_gene_density,'stroke',"black",'stroke-width',"black");

#====clomun
$svg->line('x1',$XOFFSET,'y1',$y1_gene_density,'x2',$XOFFSET,'y2',$y1_gene_density-$max_gene_number*$Y1scale,'stroke',"black",'stroke-width',"black");
$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
$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
#=========================================================================================
#=============================plot TE gene desity=========================================
#===row
=cut
my $Y2scale=$Y1scale;
my $y2_te_gene_density=$y1_gene_density;
for (my $i=0;$i<@te_gene_desity-1 ;$i++) {
	my $te_density_start_x=$XOFFSET+$i*$span*$Xscale;
	my $te_density_start_y=$y2_te_gene_density+$te_gene_desity[$i]*$Y2scale;
	my $te_density_end_x=$XOFFSET+($i+1)*$span*$Xscale;
	my $te_density_end_y=$y2_te_gene_density+$te_gene_desity[$i+1]*$Y2scale;
	#my $te_heigth=$te_gene_desity[$i]/$max_gene_number;
	$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");
	#$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'});
	
}
#column
$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");
$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
=cut
#=======================gene density txt==================================================
my $md=$span/1000;
#$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");

$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");
$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");

#$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");
#$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");

$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");
#=========================================================================================
#=============================plot exp gene desity=========================================
=cut
my $y3_target_e_density=$y2_te_gene_density+$max_te_gene_number*$Y2scale+$max_target_e_number*$Y2scale+20;
#my $Y3scale=$Y1scale;
for (my $i=0;$i<@target_e_desity-1 ;$i++) {
	my $target_e_density_start_x=$XOFFSET+$i*$span*$Xscale;
	my $target_e_density_start_y=$y3_target_e_density-$target_e_desity[$i]*$Y2scale;
	my $target_e_density_end_x=$XOFFSET+($i+1)*$span*$Xscale;
	my $target_e_density_end_y=$y3_target_e_density-$target_e_desity[$i+1]*$Y2scale;
	$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");

}
$svg->line('x1',$XOFFSET,'y1',$y3_target_e_density,'x2',$XOFFSET+$length*$Xscale,'y2',$y3_target_e_density,'stroke',"black",'stroke-width',"black");
#column
$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");
$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
$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
#=========================================exp gene indensity txt==========================
$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");
$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");
$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");
#calculate the different
=cut
#========================================================================================
my $Y3scale=0.04/$max_cluster_density*2500;
#my $y3_cluster_density=$y2_te_gene_density+$max_te_gene_number*5+$max_cluster_density*$Y3scale+10;
my $y3_cluster_density=$y1_gene_density+20;
my $y3_total=$y1_gene_density+10;
my @cluster_color=qw(fuchsia violet tomato);
for (my $s=0;$s<3 ;$s++) {
	$y3_total=$y3_total+$max_cluster_density*$Y3scale+5;
	$svg->line('x1',$XOFFSET,'y1',$y3_total,'x2',$XOFFSET+$length*$Xscale,'y2',$y3_total,'stroke',"black",'stroke-width',$attribute{csv}{'stroke-width'});
	$svg->line('x1',$XOFFSET,'y1',$y3_total,'x2',$XOFFSET,'y2',$y3_total-$max_cluster_density*$Y3scale,'stroke',"black",'stroke-width',$attribute{csv}{'stroke-width'});

	$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

	if ($max_cluster_density[$s]>$max_cluster_density) {
		$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
	}
	else{
		$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
	}
	
	for (my $i=0;$i<$#cluster_density ;$i++) {
		my $cluster_density_start_x=$XOFFSET+$i*$span*$Xscale;
		my $cluster_density_end_x=$XOFFSET+($i+1)*$span*$Xscale;
		if ($cluster_density[$i][$s]>$max_cluster_density) {
			$cluster_density[$i][$s]=$max_cluster_density;
		}
		if ($cluster_density[$i+1][$s]>$max_cluster_density) {
			$cluster_density[$i+1][$s]=$max_cluster_density;
		}
		my $cluster_density_start_y=$y3_total-$cluster_density[$i][$s]*$Y3scale;
		my $cluster_density_end_y=$y3_total-$cluster_density[$i+1][$s]*$Y3scale;
		$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'});
	}
}
for (my $s=0;$s<@sample ;$s++) {
	$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'});
	$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]);
}

$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");
#===================================plot region target gene rpkm================
=cut
my $y4_region_gene_rpkm_1=$y3_total+10;
my $y4_region_gene_rpkm=$y4_region_gene_rpkm_1;
my $Y4scale=0.1/$max_region_gene_rpkm*1000;
my @cluster_color_t=qw(blue slateblue steelblue);
for (my $s=0;$s<3 ;$s++) {
	$y4_region_gene_rpkm+=$max_region_gene_rpkm*$Y4scale+5;

	$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'});
	
	$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'});

	$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

	if ($max_region_gene_rpkm[$s]>$max_region_gene_rpkm) {
		$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
	}
	else{
		$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
	}
	for (my $i=0;$i<$#region_target_rpkm ;$i++) {
		my $region_target_rpkm_start_x=$XOFFSET+$i*$span*$Xscale;
		my $region_target_rpkm_end_x=$XOFFSET+($i+1)*$span*$Xscale;
		if ($region_target_rpkm[$i][$s]>$max_region_gene_rpkm) {
			$region_target_rpkm[$i][$s]=$max_region_gene_rpkm;
		}
		if ($region_target_rpkm[$i+1][$s]>$max_region_gene_rpkm) {
			$region_target_rpkm[$i+1][$s]=$max_region_gene_rpkm;
		}
		my $region_target_rpkm_start_y=$y4_region_gene_rpkm-$region_target_rpkm[$i][$s]*$Y4scale;
		my $region_target_rpkm_end_y=$y4_region_gene_rpkm-$region_target_rpkm[$i+1][$s]*$Y4scale;
		$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'});
	}

}
for (my $s=0;$s<3 ;$s++) {
	$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'});
	$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]);
}
$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");
=cut
#for (my $s=0;$s<3 ;$s++) {
#	$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'});
#	$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]);
#}
#=========================================================================================

#sub ExpG_number {
#
#}
#sub ExpCluster_number{
#
#}


#========================================================================================
open (OUT,">$opt{o}");
print OUT $svg->xmlify();

sub log2 {
	my $n = shift;
	return log($n)/log(2);
}

sub usage{
print <<"USAGE";
Version $version
Usage:
$0 
options:
-l centromere length
-chro
-mark \#
-g input file of Gene region annotation which format like GenePred
-span
-c cluster file input
-o svg output
-h help
USAGE
exit(1);
}