Mercurial > repos > rliterman > csp2

comparison CSP2/CSP2_env/env-d9b9114564458d9d-741b3de822f2aaca6c6caa4325c4afce/opt/mummer-3.23/nucmer2xfig @ 69:33d812a61356

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload commit 2e9511a184a1ca667c7be0c6321a36dc4e3d116d
author jpayne
date Tue, 18 Mar 2025 17:55:14 -0400
parents
children
comparison
equal deleted inserted replaced
67:0e9998148a16 69:33d812a61356
1 #!/usr/bin/env perl
2 # (c) Steven Salzberg 2001
3 # Make an xfig plot for a comparison of a reference chromosome (or single
4 # molecule) versus a multifasta file of contigs from another genome.
5 # The input file here is a NUCmer coords file such as:
6 # 2551 2577 | 240 266 | 27 27 | 96.30 | 20302755 1424 | 0.00 1.90 | 2R 1972084
7 # generated by running 'show-coords -c -l <delta file>'
8 # For the above example, D. melanogaster chr 2R is the reference and the query
9 # is an assembly with 1000s of contigs from D. pseudoobscura
10 # The file needs to be sorted by the smaller contig ids, via:
11 # tail +6 Dmel2R-vs-Dpseudo.coords | sort -k 19n -k 4n > Dmel2R-vs-Dpseudo-resort.coords
12 # and remember to get rid of top 5 (header) lines.
13 # Usage: plot-drosoph-align-xfig.perl Dmel2R-vs-Dpseudo-resort.coords
14 unless (open(coordsfile,$ARGV[0])) {
15 die ("can't open file $ARGV[0].\n");
16 }
17 $Xscale = 0.005;
18 $Yscale = 20;
19 $chrcolor = 4; # 4 is red, 1 is blue, 2 is green
20 $green = 2;
21 $contigcolor = 1; # contigs are blue
22
23 # print header info
24 print "#FIG 3.2\nLandscape\nCenter\nInches\nLetter \n100.00\nSingle\n-2\n1200 2\n";
25 $first_time = 1;
26 while (<coordsfile>) {
27 ($s1,$e1,$x1,$s2,$e2,$x2,$l1,$l2,$x3,$percent_id,$x4,$Rlen,$Qlen,$x5,$Rcov,$Qcov,$x6,$Rid,$Qid) = split(" ");
28 if ($prevQid eq $Qid) { # query contig is same as prev line
29 $dist = abs($s1 - $prev_s1);
30 if ( $dist > 2 * $Qlen) { # if this match is too far away
31 # print the contig here if the matching bit is > 1000
32 if ($right_end{$Qid} - $left_end{$Qid} > 1000) {
33 # print it at y=50 rather than 100 because the scale is 0-50,
34 # where 0=50% identical and 50 is 100% identical
35 print_xfig_line($left_end{$Qid},50,$right_end{$Qid},50,$contigcolor);
36 print_label($left_end{$Qid},50,$right_end{$Qid},$Qid,5);
37 }
38 $left_end{$Qid} = $s1; # then re-set the start and end of the contig
39 $right_end{$Qid} = $e1; # and we'll print it again later
40 }
41 else { # extend the boundaries of the match
42 if ($s1 < $left_end{$Qid}) { $left_end{$Qid} = $s1; }
43 if ($e1 > $right_end{$Qid}) { $right_end{$Qid} = $e1; }
44 }
45 }
46 else { # this is a different contig, first time seeing it
47 $left_end{$Qid} = $s1;
48 $right_end{$Qid} = $e1;
49 # print the previous contig as a line at y=100 for 100%
50 if ($first_time < 1) {
51 print_xfig_line($left_end{$prevQid},50,$right_end{$prevQid},50,$contigcolor);
52 print_label($left_end{$prevQid},50,$right_end{$prevQid},$prevQid,5);
53 }
54 else { $first_time = 0; }
55 }
56 $prevQid = $Qid;
57 $prev_s1 = $s1;
58 $prev_Qlen = $Qlen;
59 # next print the matching bit as a separate line, with a separate color,
60 # with its height determined by percent match
61 $Xleft = int($Xscale * $s1);
62 $Xright = int($Xscale * $e1);
63 if ($Xleft == $Xright) { $Xright += 1; }
64 if ($percent_id < 50) { $percent_id = 50; }
65 print_xfig_line($s1,$percent_id - 50,$e1,$percent_id - 50,$green);
66 }
67 # print very last contig
68 $left_end{$Qid} = $s1;
69 $right_end{$Qid} = $e1;
70 print_xfig_line($s1,50,$e1,50,$contigcolor);
71 print_label($s1,50,$e1,$Qid,5);
72 close(coordsfile);
73
74 # now draw the horizontal chr line for the reference
75 $label_xpos = $Xscale * $Rlen;
76 print "4 0 0 100 0 0 12 0.0000 4 135 405 ";
77 printf("%.0f %.0f",$label_xpos, 0);
78 print " ", $Rid, "\\001\n";
79 print "2 1 0 2 $chrcolor 7 100 0 -1 0.000 0 0 -1 0 0 2\n";
80 printf("\t %.0f %.0f %.0f %.0f\n", 0, 0,
81 $Xscale * $Rlen, 0);
82 # print some X-axis coordinates
83 $pointsize = 5;
84 for ($i = 250000; $i < $Rlen - 50000; $i+= 250000) {
85 print "4 0 0 100 0 0 $pointsize 0.0000 4 135 405 ";
86 printf("%.0f %.0f",$i * $Xscale + 20, -20 + ($Yscale * -0.5));
87 print " $i", "\\001\n";
88 #print a vertical tic mark
89 print "2 1 0 1 0 7 50 0 -1 0.000 0 0 -1 0 0 2\n";
90 printf("%.0f %.0f %.0f -50\n",
91 $i * $Xscale, 0, $i * $Xscale);
92 }
93 # print tic marks indicating % identity on the y-axis
94 for ($percent_id = 50; $percent_id < 101; $percent_id += 10) {
95 print "4 0 0 100 0 0 $pointsize 0.0000 4 135 405 ";
96 printf("-150 %.0f", ($percent_id - 50) * $Yscale + 10); # shift down 10 pixels
97 print " $percent_id", "\\001\n";
98 # print the tic mark
99 print "2 1 0 1 0 7 50 0 -1 0.000 0 0 -1 0 0 2\n";
100 printf("-50 %.0f 0 %.0f\n",
101 ($percent_id - 50) * $Yscale, ($percent_id - 50) * $Yscale);
102 }
103
104 # print a line in the appropriate color, scaled with Xscale,Yscale
105 sub print_xfig_line {
106 my ($xleft,$yleft,$xright,$yright,$color) = @_;
107 # print it at the given coordinates, scaled
108 $xleft_scaled = int($Xscale * $xleft);
109 $xright_scaled = int($Xscale * $xright);
110 # Xfig has a bug: if we re-scale and the resulting X coordinates are equal,
111 # then it will print a fixed-size (large) rectangle, rather than a 1-pixel
112 # wide one. So check fo this and correct.
113 if ($xleft_scaled == $xright_scaled) { $xright_scaled += 1; }
114 # set up and print line in xfig format
115 print "2 1 0 2 $color 7 100 0 -1 0.000 0 0 -1 0 0 2\n";
116 printf("\t %.0f %.0f %.0f %.0f\n",
117 $xleft_scaled, $Yscale * $yleft, $xright_scaled, $Yscale * $yright);
118 }
119
120 # print a label for each contig using its ID
121 sub print_label {
122 my ($xleft,$yleft,$xright,$id,$psize) = @_;
123 # print it at the left edge of the contig. The angle
124 # of 4.7124 means text goes down vertically. 5th argument
125 # here is pointsize of the label.
126 # bump the label down 20 pixels
127 $yposition = $yleft * $Yscale + 20;
128 # to keep the display clean, don't print the label unless the
129 # contig itself is wider than the width of the text in the label
130 $xleft_scaled = $xleft * $Xscale;
131 $xright_scaled = $xright * $Xscale;
132 # each point of type is 8 pixels
133 $textheight = $psize * 8;
134 if ($xright_scaled - $xleft_scaled > $textheight) {
135 print "4 0 0 50 0 0 $psize 4.7124 4 135 435 ";
136 printf("%.0f %.0f",$xleft_scaled, $yposition);
137 print " $id", "\\001\n";
138 }
139 }