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

Annotation of /FigKernelScripts/svr_link_to_compare_regions.pl

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : overbeek 1.1 use strict;
2 :     use Data::Dumper;
3 :     use Carp;
4 :    
5 :     #
6 :     # This is a SAS Component
7 :     #
8 :    
9 :     =head1 svr_link_to_compare_regions
10 :    
11 :     Get link to compare regions organized to show co-occuring genes
12 :    
13 :     ------
14 :    
15 :     Example:
16 :    
17 :     echo 'fig|83333.1.peg.4' | svr_link_to_compare_regions > 2-col.txt
18 :    
19 :     would produce a 2-column table (in this case containing a single row).
20 :     The first column would contain a PEG id (in this case fig|83333.1.peg.4)
21 :     and the second a URL to illustrate conserved contiguity.
22 :    
23 :     ------
24 :    
25 :     The standard input should be a tab-separated table (i.e., each line
26 :     is a tab-separated set of fields). Normally, the last field in each
27 :     line would contain the PEG for which functions are being requested.
28 :     If some other column contains the PEGs, use
29 :    
30 :     -c N
31 :    
32 :     where N is the column (from 1) that contains the PEG in each case.
33 :    
34 :     This is a pipe command. The input is taken from the standard input, and the
35 :     output is to the standard output.
36 :    
37 :     =head2 Command-Line Options
38 :    
39 :     =over 4
40 :    
41 :     =item -c Column
42 :    
43 :     This is used only if the column containing PEGs is not the last.
44 :    
45 :     =item -b Base of URL (to the appropriate SEED) [default=http://pubseed.theseed.org which is the PubSEED]
46 :    
47 :     This allows the user to specify which SEED he wishes the links to point at
48 :    
49 :     =item -m Minimum Functional Coupling Score [default=10]
50 :    
51 :     This parameter lets the user give a minimum FC score (which is the number of
52 :     distinct OTUs in which co-occurrence has been detected).
53 :    
54 :     =item -g GenomesFile [an optional file designating a restricted list of genomes for the links]
55 :    
56 :     This allows the user to restrict the compare region to a limited set of genomes.
57 :    
58 :     =item -show N [default is 10]
59 :    
60 :     This gives the number of genomes to actually show.
61 :    
62 :     =back
63 :    
64 :     =head2 Output Format
65 :    
66 :     The standard output is a tab-delimited file. It consists of the input
67 :     file with an extra column added (the URL to produce the desired compare region
68 :    
69 :     =cut
70 :    
71 :     use SeedUtils;
72 :     use SAPserver;
73 :     my $sapO = SAPserver->new();
74 :     use Getopt::Long;
75 :    
76 :     my $usage = "usage: svr_link_to_compare_regions [-c column] [-m MinFC] [-b BaseURL] [-g GenomesFile] [-show N]";
77 :    
78 :     my $column;
79 :     my $minFC = 10;
80 :     my $show = 10;
81 :     my $base = "http://pubseed.theseed.org";
82 :     my $genomeF;
83 :    
84 :     my $rc = GetOptions('c=i' => \$column,
85 :     'm=i' => \$minFC,
86 :     'b=s' => \$base,
87 :     'g=s' => \$genomeF,
88 :     'show=s' => \$show
89 :     );
90 :     if (! $rc) { print STDERR $usage; exit }
91 :     my @lines = map { chomp; [split(/\t/,$_)] } <STDIN>;
92 :     (@lines > 0) || exit;
93 :     if (! $column) { $column = @{$lines[0]} }
94 :     my @fids = map { $_->[$column-1] } @lines;
95 :     my $fcH = $sapO->conserved_in_neighborhood( -ids => \@fids);
96 :    
97 :     my %genomeH;
98 :     if ($genomeF && (-s $genomeF) && open(G,"<",$genomeF))
99 :     {
100 :     %genomeH = map { ($_ =~ /(\d+\.\d+)/) ? ($1 => 1) : () } <G>;
101 :     close(G);
102 :     }
103 :    
104 :     my @pairs;
105 :     my %urlH;
106 :     foreach my $peg (keys(%$fcH))
107 :     {
108 :     my $hits = $fcH->{$peg};
109 :     my @fc = map { $_->[1] }
110 :     sort { $b->[0] <=> $a->[0] }
111 :     grep { $_->[0] >= $minFC } @{$hits};
112 :     if (@fc)
113 :     {
114 :     push(@pairs,[$peg,[@fc]]);
115 :     }
116 :     }
117 :     my @flattened = map { my($p,$xL) = @$_; map { join(":",($p,$_)) } @$xL } @pairs;
118 :     my $pairH = $sapO->co_occurrence_evidence( -pairs => \@flattened);
119 :     foreach my $pair (@pairs)
120 :     {
121 :     my($peg1,$pegs2) = @$pair;
122 :     my %anchored;
123 :     my $i;
124 :     for ($i=0; ($i < @$pegs2); $i++)
125 :     {
126 :     my $peg2 = $pegs2->[$i];
127 :     my $ev = $pairH->{"$peg1:$peg2"};
128 :     if ($ev)
129 :     {
130 :     foreach my $pair2 (@$ev)
131 :     {
132 :     if ((! $anchored{$pair2->[0]}) &&
133 :     ((! $genomeF) || &genomeH->{&SeedUtils::genome_of($pair2->[0])}))
134 :     {
135 :     $anchored{$pair2->[0]} = $i;
136 :     }
137 :     }
138 :     }
139 :     }
140 :     my @pin = ($peg1,sort { $anchored{$b} <=> $anchored{$b} } keys(%anchored));
141 :     if (@pin > $show) { $#pin = $show-1 }
142 : overbeek 1.2 my $url = "$base/seedviewer.cgi?page=Regions&" . join("&",map { "feature=$_" } @pin);
143 : overbeek 1.1 $urlH{$peg1} = $url;
144 :     }
145 :    
146 :     foreach $_ (@lines)
147 :     {
148 :     my $peg = $_->[$column-1];
149 :     my $url = $urlH{$peg};
150 :     if (! $url) { $url = '' }
151 :     print join("\t",(@$_,$url)),"\n";
152 :     }

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3