Mercurial > repos > rliterman > csp2
diff CSP2/CSP2_env/env-d9b9114564458d9d-741b3de822f2aaca6c6caa4325c4afce/opt/mummer-3.23/nucmer2xfig @ 69:33d812a61356
planemo upload commit 2e9511a184a1ca667c7be0c6321a36dc4e3d116d
author | jpayne |
---|---|
date | Tue, 18 Mar 2025 17:55:14 -0400 |
parents | |
children |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/CSP2/CSP2_env/env-d9b9114564458d9d-741b3de822f2aaca6c6caa4325c4afce/opt/mummer-3.23/nucmer2xfig Tue Mar 18 17:55:14 2025 -0400 @@ -0,0 +1,139 @@ +#!/usr/bin/env 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"; + } +}