[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.1 - (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 :     my $sims = "$simsD/$fromO-$toO";
73 :    
74 :     &FIG::run("$FIG_Config::bin/sims_between $from $to p 98 0.7 > $sims");
75 :    
76 :     print "orgdir=$orgdir\n";
77 :     open(SETS, "$FIG_Config::bin/connections_between $orgdir $simsD 2 1 |") or
78 :     die "Cannot open pipe from $FIG_Config::bin/connections_between $simsD 2 1: $!";
79 :    
80 :     my %mapping;
81 :    
82 :     while (<SETS>)
83 :     {
84 :     chomp;
85 :     if ($_ =~ /^(fig\|(\d+\.\d+)\.peg\.\d+)\t(fig\|(\d+\.\d+)\.peg\.\d+)$/)
86 :     {
87 :     if (($2 eq $fromO) && ($4 eq $toO))
88 :     {
89 :     $mapping{$1} = $3;
90 :     }
91 :     elsif (($2 eq $toO) && ($4 eq $fromO))
92 :     {
93 :     $mapping{$3} = $1;
94 :     }
95 :     else
96 :     {
97 :     die "INVALID ORGS: $_";
98 :     }
99 :     }
100 :     }
101 :     close(SETS);
102 :    
103 :     my $x;
104 :     foreach $_ (sort { $a =~ /(\d+)$\t/; $x = $1; $b =~ /(\d+)\t/; $x <=> $1 } keys(%mapping))
105 :     {
106 :     print MAPS join("\t",($_,$mapping{$_})),"\n";
107 :     }
108 :     close(MAPS);
109 :     system("rm", "-r", $simsD);
110 :    
111 :     sub find_genome
112 :     {
113 :     my($fa) = @_;
114 :     my $g;
115 :     open(F, "<$fa") or die "Cannot open fasta file $fa: $!";
116 :    
117 :     $_ = <F>;
118 :     if (/^>fig\|(\d+\.\d+)\.peg/)
119 :     {
120 :     $g = $1;
121 :     }
122 :     else
123 :     {
124 :     die "Invalid fasta file $fa";
125 :     }
126 :     while (<F>)
127 :     {
128 :     if (/^>(\S+)/)
129 :     {
130 :     my $id = $1;
131 :     if ($id !~ /^fig\|$g\.peg/)
132 :     {
133 :     die "Invalid id $id found in $fa";
134 :     }
135 :     }
136 :     }
137 :     close(F);
138 :     return $g;
139 :     }

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3