mummer 3.23+dfsg-3 / usr / bin / nucmer2xfig

#!/usr/bin/perl
# (c) Steven Salzberg 2001
# Make an xfig plot for a comparison of a reference chromosome (or single
# molecule) versus a multifasta file of contigs from another genome.  
# The input file here is a NUCmer coords file such as:
# 2551  2577 | 240 266 | 27  27 | 96.30  | 20302755  1424  | 0.00  1.90 | 2R 1972084
# generated by running 'show-coords -c -l <delta file>'
# For the above example, D. melanogaster chr 2R is the reference and the query
# is an assembly with 1000s of contigs from D. pseudoobscura
# The file needs to be sorted by the smaller contig ids, via:
# tail +6 Dmel2R-vs-Dpseudo.coords | sort -k 19n -k 4n  > Dmel2R-vs-Dpseudo-resort.coords 
#  and remember to get rid of top 5 (header) lines.
# Usage: plot-drosoph-align-xfig.perl Dmel2R-vs-Dpseudo-resort.coords 
unless (open(coordsfile,$ARGV[0]))  {
    die ("can't open file $ARGV[0].\n");
}
$Xscale = 0.005;
$Yscale = 20;
$chrcolor = 4;   # 4 is red, 1 is blue, 2 is green
$green = 2;
$contigcolor = 1; # contigs are blue

# print header info
print "#FIG 3.2\nLandscape\nCenter\nInches\nLetter  \n100.00\nSingle\n-2\n1200 2\n";
$first_time = 1;
while (<coordsfile>) {
    ($s1,$e1,$x1,$s2,$e2,$x2,$l1,$l2,$x3,$percent_id,$x4,$Rlen,$Qlen,$x5,$Rcov,$Qcov,$x6,$Rid,$Qid) = split(" ");
    if ($prevQid eq $Qid) {  # query contig is same as prev line
	$dist = abs($s1 - $prev_s1);
	if ( $dist >  2 * $Qlen)  { # if this match is too far away
	    # print the contig here if the matching bit is > 1000
	    if ($right_end{$Qid} - $left_end{$Qid} > 1000) {
		# print it at y=50 rather than 100 because the scale is 0-50,
		# where 0=50% identical and 50 is 100% identical
		print_xfig_line($left_end{$Qid},50,$right_end{$Qid},50,$contigcolor);
		print_label($left_end{$Qid},50,$right_end{$Qid},$Qid,5);
	    }
	    $left_end{$Qid} = $s1;    # then re-set the start and end of the contig
	    $right_end{$Qid} = $e1;   # and we'll print it again later
	}
	else {   # extend the boundaries of the match
	    if ($s1 < $left_end{$Qid}) { $left_end{$Qid} = $s1; }
	    if ($e1 > $right_end{$Qid}) { $right_end{$Qid} = $e1; }
	}
    }
    else { # this is a different contig, first time seeing it
	$left_end{$Qid} = $s1;
	$right_end{$Qid} = $e1;
	# print the previous contig as a line at y=100 for 100%
	if ($first_time < 1) {
	    print_xfig_line($left_end{$prevQid},50,$right_end{$prevQid},50,$contigcolor);
	    print_label($left_end{$prevQid},50,$right_end{$prevQid},$prevQid,5);
	}
	else { $first_time = 0; }
    }
    $prevQid = $Qid;
    $prev_s1 = $s1;
    $prev_Qlen = $Qlen;
    # next print the matching bit as a separate line, with a separate color,
    # with its height determined by percent match
    $Xleft = int($Xscale * $s1);
    $Xright = int($Xscale * $e1);
    if ($Xleft == $Xright) { $Xright += 1; }
    if ($percent_id < 50) { $percent_id = 50; }
    print_xfig_line($s1,$percent_id - 50,$e1,$percent_id - 50,$green);
}
# print very last contig
$left_end{$Qid} = $s1;
$right_end{$Qid} = $e1;
print_xfig_line($s1,50,$e1,50,$contigcolor);
print_label($s1,50,$e1,$Qid,5);
close(coordsfile);

# now draw the horizontal chr line for the reference
$label_xpos = $Xscale * $Rlen;
print "4 0 0 100 0 0 12 0.0000 4 135 405 ";
printf("%.0f %.0f",$label_xpos, 0);
print " ", $Rid, "\\001\n";
print "2 1 0 2 $chrcolor 7 100 0 -1 0.000 0 0 -1 0 0 2\n";		
printf("\t %.0f %.0f %.0f %.0f\n", 0, 0,
       $Xscale * $Rlen, 0);
# print some X-axis coordinates
$pointsize = 5;
for ($i = 250000; $i < $Rlen - 50000; $i+= 250000) { 
    print "4 0 0 100 0 0 $pointsize 0.0000 4 135 405 ";
    printf("%.0f %.0f",$i * $Xscale + 20, -20 + ($Yscale * -0.5));
    print " $i", "\\001\n";
    #print a vertical tic mark
    print "2 1 0 1 0 7 50 0 -1 0.000 0 0 -1 0 0 2\n"; 
    printf("%.0f %.0f %.0f -50\n",
	   $i * $Xscale, 0, $i * $Xscale);
}
# print tic marks indicating % identity on the y-axis
for ($percent_id = 50; $percent_id < 101; $percent_id += 10) {
    print "4 0 0 100 0 0 $pointsize 0.0000 4 135 405 ";
    printf("-150 %.0f", ($percent_id - 50) * $Yscale + 10); # shift down 10 pixels
    print " $percent_id", "\\001\n";
    # print the tic mark
    print "2 1 0 1 0 7 50 0 -1 0.000 0 0 -1 0 0 2\n"; 
    printf("-50 %.0f 0 %.0f\n",
	   ($percent_id - 50) * $Yscale, ($percent_id - 50) * $Yscale);
}

# print a line in the appropriate color, scaled with Xscale,Yscale
sub print_xfig_line {
    my ($xleft,$yleft,$xright,$yright,$color) = @_;
    # print it at the given coordinates, scaled
    $xleft_scaled = int($Xscale * $xleft);
    $xright_scaled = int($Xscale * $xright);
    # Xfig has a bug: if we re-scale and the resulting X coordinates are equal,
    # then it will print a fixed-size (large) rectangle, rather than a 1-pixel
    # wide one.  So check fo this and correct.
    if ($xleft_scaled == $xright_scaled) { $xright_scaled += 1; }
    # set up and print line in xfig format
    print "2 1 0 2 $color 7 100 0 -1 0.000 0 0 -1 0 0 2\n";		
    printf("\t %.0f %.0f %.0f %.0f\n", 
	   $xleft_scaled, $Yscale * $yleft, $xright_scaled, $Yscale * $yright);
}

# print a label for each contig using its ID
sub print_label {
    my ($xleft,$yleft,$xright,$id,$psize) = @_;
    # print it at the left edge of the contig.  The angle
    # of 4.7124 means text goes down vertically.  5th argument
    # here is pointsize of the label.  
    # bump the label down 20 pixels
    $yposition = $yleft * $Yscale + 20;
    # to keep the display clean, don't print the label unless the
    # contig itself is wider than the width of the text in the label
    $xleft_scaled = $xleft * $Xscale;
    $xright_scaled = $xright * $Xscale;
    # each point of type is 8 pixels
    $textheight = $psize * 8;
    if ($xright_scaled - $xleft_scaled > $textheight) {
	print "4 0 0 50 0 0 $psize 4.7124 4 135 435 ";
	printf("%.0f %.0f",$xleft_scaled, $yposition);
	print " $id", "\\001\n";
    }
}