[Bio] / FigKernelScripts / make_peg_maps_from_fasta.pl Repository:
ViewVC logotype

Annotation of /FigKernelScripts/make_peg_maps_from_fasta.pl

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.2 - (view) (download) (as text)

1 : olson 1.1 #
2 :     # Copyright (c) 2003-2006 University of Chicago and Fellowship
3 :     # for Interpretations of Genomes. All Rights Reserved.
4 :     #
5 :     # This file is part of the SEED Toolkit.
6 :     #
7 :     # The SEED Toolkit is free software. You can redistribute
8 :     # it and/or modify it under the terms of the SEED Toolkit
9 :     # Public License.
10 :     #
11 :     # You should have received a copy of the SEED Toolkit Public License
12 :     # along with this program; if not write to the University of Chicago
13 :     # at info@ci.uchicago.edu or the Fellowship for Interpretation of
14 :     # Genomes at veronika@thefig.info or download a copy from
15 :     # http://www.theseed.org/LICENSE.TXT.
16 :     #
17 :    
18 :     #
19 :     # Script derived from make_peg_maps, but optimized for the case of comparing
20 :     # two fasta files between which we wish to map.
21 :     #
22 :    
23 :     use FIGV;
24 :     use FIG;
25 :     use strict;
26 :    
27 :     my $usage = "usage: make_peg_maps fasta-from fasta-to mapfile";
28 :    
29 :     my($from, $to, $maps);
30 :    
31 :     (
32 :     ($from = shift @ARGV) &&
33 :     ($to = shift @ARGV) &&
34 :     ($maps = shift @ARGV)
35 :     )
36 :     || die $usage;
37 :    
38 :     #
39 :     # Determine the genome ID for each, and ensure that all the sequences
40 :     # in each file has the same id.
41 :     #
42 :    
43 :     my($fromO, $toO);
44 :    
45 :     if ($from =~ /(\d+\.\d+)$/ and -f "$from/GENOME")
46 :     {
47 :     $fromO = $1;
48 :     $from = "$from/Features/peg/fasta";
49 :     }
50 :     else
51 :     {
52 :     $fromO = find_genome($from);
53 :     }
54 :    
55 :     my $orgdir;
56 :     if ($to =~ /(\d+\.\d+)$/ and -f "$to/GENOME")
57 :     {
58 :     $toO = $1;
59 :     $orgdir = "-orgdir $to";
60 :     $to = "$to/Features/peg/fasta";
61 :     }
62 :     else
63 :     {
64 :     $toO = find_genome($to);
65 :     }
66 :    
67 :     open(MAPS, ">$maps") or die "Cannot open $maps for writing: $!";
68 :    
69 :     my $simsD = "$FIG_Config::temp/mapsims.$$";
70 :     mkdir($simsD, 0777) || die "could not mkdir $simsD: $!";
71 :    
72 : olson 1.2 my $sims = "$simsD/$toO-$fromO";
73 :    
74 :     &FIG::run("$FIG_Config::bin/sims_between $to $from p 98 0.7 > $sims");
75 :    
76 : olson 1.1 my $sims = "$simsD/$fromO-$toO";
77 :    
78 :     &FIG::run("$FIG_Config::bin/sims_between $from $to p 98 0.7 > $sims");
79 :    
80 :     print "orgdir=$orgdir\n";
81 :     open(SETS, "$FIG_Config::bin/connections_between $orgdir $simsD 2 1 |") or
82 :     die "Cannot open pipe from $FIG_Config::bin/connections_between $simsD 2 1: $!";
83 :    
84 :     my %mapping;
85 :    
86 :     while (<SETS>)
87 :     {
88 :     chomp;
89 :     if ($_ =~ /^(fig\|(\d+\.\d+)\.peg\.\d+)\t(fig\|(\d+\.\d+)\.peg\.\d+)$/)
90 :     {
91 :     if (($2 eq $fromO) && ($4 eq $toO))
92 :     {
93 :     $mapping{$1} = $3;
94 :     }
95 :     elsif (($2 eq $toO) && ($4 eq $fromO))
96 :     {
97 :     $mapping{$3} = $1;
98 :     }
99 :     else
100 :     {
101 :     die "INVALID ORGS: $_";
102 :     }
103 :     }
104 :     }
105 :     close(SETS);
106 :    
107 :     my $x;
108 :     foreach $_ (sort { $a =~ /(\d+)$\t/; $x = $1; $b =~ /(\d+)\t/; $x <=> $1 } keys(%mapping))
109 :     {
110 :     print MAPS join("\t",($_,$mapping{$_})),"\n";
111 :     }
112 :     close(MAPS);
113 :     system("rm", "-r", $simsD);
114 :    
115 :     sub find_genome
116 :     {
117 :     my($fa) = @_;
118 :     my $g;
119 :     open(F, "<$fa") or die "Cannot open fasta file $fa: $!";
120 :    
121 :     $_ = <F>;
122 :     if (/^>fig\|(\d+\.\d+)\.peg/)
123 :     {
124 :     $g = $1;
125 :     }
126 :     else
127 :     {
128 :     die "Invalid fasta file $fa";
129 :     }
130 :     while (<F>)
131 :     {
132 :     if (/^>(\S+)/)
133 :     {
134 :     my $id = $1;
135 :     if ($id !~ /^fig\|$g\.peg/)
136 :     {
137 :     die "Invalid id $id found in $fa";
138 :     }
139 :     }
140 :     }
141 :     close(F);
142 :     return $g;
143 :     }

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3