[Bio] / FigWebServices / index.cgi Repository:
ViewVC logotype

Diff of /FigWebServices/index.cgi

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1.64, Wed Mar 16 17:52:34 2005 UTC revision 1.65, Sun Mar 27 22:20:39 2005 UTC
# Line 2  Line 2 
2    
3  use FIG;  use FIG;
4    
5    use strict;
6    
7    use FIGjs        qw( toolTipScript );
8    use GenoGraphics qw( render );
9    use IPC::Open2   qw( open2 );
10    
11  use POSIX;  use POSIX;
12  use HTML;  use HTML;
13  use strict;  
14  use CGI;  use CGI;
15  my $cgi = new CGI;  my $cgi = new CGI;
16    
# Line 979  Line 985 
985          my $db = "$FIG_Config::organisms/$org/contigs";          my $db = "$FIG_Config::organisms/$org/contigs";
986          &verify_db( $db, "n" );                               ### fix to get all contigs          &verify_db( $db, "n" );                               ### fix to get all contigs
987          @out = `$blastall -i $tmp_seq -d $db -p blastn -r 1 -q -1 $blast_opt`;          @out = `$blastall -i $tmp_seq -d $db -p blastn -r 1 -q -1 $blast_opt`;
988            push @$html, blast_graphics( $fig, $cgi, $org, \@out, $tool );
989      }      }
990    
991      elsif ( $tool eq "tblastn" )      elsif ( $tool eq "tblastn" )
# Line 986  Line 993 
993          my $db = "$FIG_Config::organisms/$org/contigs";          my $db = "$FIG_Config::organisms/$org/contigs";
994          &verify_db( $db, "n" );                               ### fix to get all contigs          &verify_db( $db, "n" );                               ### fix to get all contigs
995          @out = `$blastall -i $tmp_seq -d $db -p tblastn $blast_opt`;          @out = `$blastall -i $tmp_seq -d $db -p tblastn $blast_opt`;
996            push @$html, blast_graphics( $fig, $cgi, $org, \@out, $tool );
997      }      }
998    
999      elsif ( $tool eq 'blastp against complete genomes' )     ### this tool gets nonstandard treatment: RAO      elsif ( $tool eq 'blastp against complete genomes' )     ### this tool gets nonstandard treatment: RAO
# Line 1031  Line 1039 
1039      &format_sims($fig,$cgi,$html,\@sims);      &format_sims($fig,$cgi,$html,\@sims);
1040  }  }
1041    
1042    
1043    #------------------------------------------------------------------------------
1044    #  Graphically display searches against contigs
1045    #
1046    #  use FIGjs        qw( toolTipScript );
1047    #  use GenoGraphics qw( render );
1048    #  use IPC::Open2   qw( open2 );
1049    #------------------------------------------------------------------------------
1050    #  Fields produced by rationalize_blast:
1051    #
1052    #  0    1      2    3    4        5       6       7      8      9   10  11   12   13  14  15
1053    # HSP  score  exp  p_n  p_val  n_match  n_ident  n_sim  n_gap  dir  q1  q2  q_sq  s1  s2  s_sq
1054    #------------------------------------------------------------------------------
1055    
1056    sub blast_graphics {
1057        my ( $fig_or_sprout, $cgi, $genome, $out, $tool ) = @_;
1058    
1059        my $e_min = 0.1;
1060        my $gg = [];
1061        my @html = ();;
1062    
1063        # Run rationalize_blast:
1064    
1065        my( $pid, $rd, $wr );
1066        if ( $pid = open2( $rd, $wr, "rationalize_blast" ) )
1067        {
1068            my $outlen = 0;
1069            foreach ( @$out ) { $outlen += length( $_ ) }
1070    
1071            $wr->write( join( "", @$out ), $outlen );
1072            close( $wr );
1073    
1074            my ( $qid, $qdef, $qlen, $contig, $sdef, $slen );
1075            my @rational = <$rd>;
1076            foreach ( map { chomp; $_ } @rational )
1077            {
1078                if    ( /^Query=/ ) { ( undef, $qid,    undef, $qlen ) = split /\t/ }
1079                elsif ( /^>/ )      { ( undef, $contig, undef, $slen ) = split /\t/ }
1080                elsif ( /^HSP/ && $qid && $qlen && $contig && $slen )
1081                {
1082                    my @hsp = split /\t/;
1083                    next if $hsp[2] > $e_min;
1084                    my ( $e_val, $q1, $q2, $s1, $s2 ) = @hsp[ 2, 10, 11, 13, 14 ];
1085                    my ( $genes, $min, $max ) = hsp_context( $fig_or_sprout, $cgi, $genome,
1086                                                             $e_val, 100 * $hsp[6] / $hsp[5],
1087                                                             $qid,    $q1, $q2, $qlen,
1088                                                             $contig, $s1, $s2, $slen
1089                                                           );
1090                    push @$gg, [ substr( $contig, 0, 18 ), $min, $max, $genes ];
1091                }
1092            }
1093            close( $rd );
1094    
1095            # $gene  = [ $beg, $end, $shape, $color, $text, $url, $pop-up, $alt_action, $pop-up_title ];
1096            # $genes = [ $gene, $gene, ... ];
1097            # $map   = [ $label, $min_coord, $max_coord, $genes ];
1098            # $gg    = [ $map, $map, ... ];
1099            # render( $gg, $width, $obj_half_heigth, $save, $img_index_number )
1100    
1101            if ( @$gg )
1102            {
1103                # print STDERR Dumper( $gg );
1104                my $gs = $fig_or_sprout->genus_species( $genome );
1105                my $space = "&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;";
1106                my $legend = "<TABLE>\n"
1107                           . "    <TR>\n"
1108                           . "        <TD>Q = Query sequence$space</TD>\n"
1109                           . "        <TD Bgcolor='#FF0000'>$space</TD><TD>Frame 1 translation$space</TD>\n"
1110                           . "        <TD Bgcolor='#00FF00'>$space</TD><TD>Frame 2 translation$space</TD>\n"
1111                           . "        <TD Bgcolor='#0000FF'>$space</TD><TD>Frame 3 translation$space</TD>\n"
1112                           . "        <TD Bgcolor='#808080'>$space</TD><TD>Untranslated sequence</TD>\n"
1113                           . "    </TR>\n"
1114                           . "</TABLE><P />";
1115    
1116                push @html, "\n", FIGjs::toolTipScript(), "\n",
1117                            $cgi->h2( "Results of $tool search of contigs from $gs\n"),
1118                            $legend,
1119                            @{ GenoGraphics::render( $gg, 600, 4, 0, 1 ) },
1120                            $cgi->hr, "\n";
1121            }
1122    
1123            waitpid $pid, 0;
1124        }
1125    
1126        return @html;
1127    }
1128    
1129    
1130    sub hsp_context {
1131        my( $fig_or_sprout, $cgi, $genome, $e_val, $pct_id,
1132            $qid,    $q1, $q2, $qlen,
1133            $contig, $s1, $s2, $slen ) = @_;
1134        my $half_sz = 5000;
1135    
1136        my( $from, $to, $features, $fid, $beg, $end );
1137        my( $link, $lbl, $isprot, $function, $uniprot, $info, $prot_query );
1138    
1139        my $user   = $cgi->param( 'user' ) || "";
1140        my $sprout = $cgi->param( 'SPROUT' ) ? '&SPROUT=1' : '';
1141    
1142        my @genes  = ();
1143    
1144        #  Based on the match position of the query, select the context region:
1145    
1146        ( $from, $to ) = ( $s1 <= $s2 ) ? ( $s1 - $half_sz, $s2 + $half_sz )
1147                                        : ( $s2 - $half_sz, $s1 + $half_sz );
1148        $from = 1      if ( $from < 1 );
1149        $to   = $slen  if ( $to > $slen );
1150    
1151        #  Get the genes in the region, and adjust the ends to include whole genes:
1152    
1153        ( $features, $from, $to ) = genes_in_region( $fig_or_sprout, $cgi, $genome, $contig, $from, $to );
1154    
1155    
1156        #  Add the other features:
1157    
1158        foreach $fid ( @$features )
1159        {
1160            my $contig1;
1161            ( $contig1, $beg, $end ) = boundaries_of( $fig_or_sprout, feature_locationS( $fig_or_sprout, $fid ) );
1162            next if $contig1 ne $contig;
1163    
1164            $link = "";
1165            if ( ( $lbl ) = $fid =~ /peg\.(\d+)$/ ) {
1166                ( $link = $cgi->url() ) =~ s/index\.cgi/protein.cgi/;
1167                $link .= "?prot=$fid&user=$user$sprout";
1168                $isprot = 1;
1169            } elsif ( ( $lbl ) = $fid =~ /\.([a-z]+)\.\d+$/ ) {
1170                $lbl = uc $lbl;
1171                $isprot = 0;
1172            } else {
1173                $lbl = "";
1174                $isprot = 0;
1175            }
1176    
1177            $function = function_ofS( $fig_or_sprout, $fid );
1178    
1179            $uniprot = join ", ", grep { /^uni\|/ } feature_aliasesL( $fig_or_sprout, $fid);
1180    
1181            $info = join( '<br />', "<b>Feature:</b> $fid",
1182                                    "<b>Contig:</b> $contig",
1183                                    "<b>Begin:</b> $beg",
1184                                    "<b>End:</b> $end",
1185                                    $function ? "<b>Function:</b> $function" : '',
1186                                    $uniprot ? "<b>Uniprot ID:</b> $uniprot" : ''
1187                        );
1188    
1189            push @genes, [ feature_graphic( $beg, $end, $isprot ),
1190                           $lbl, $link, $info,
1191                           $isprot ? () : ( undef, "Feature information" )
1192                         ];
1193        }
1194    
1195        #  Draw the query.  The subject coordinates are always DNA.  If the query
1196        #  is protein, it is about 3 times shorter than the matching contig DNA.
1197        #  Splitting the difference, if 1.7 times the query length is still less
1198        #  than the subject length, we will call it a protein query (and reading
1199        #  frame in the contig coordinates has meaning).  If it is nucleotides,
1200        #  there is no defined frame.
1201    
1202        $info = join( '<br />', $qid ne 'query ' ? "<b>Query:</b> $qid" : (),
1203                                "<b>Length:</b> $qlen",
1204                                "<b>E-value:</b> $e_val",
1205                                "<b>% identity:</b> " . sprintf( "%.1f", $pct_id ),
1206                                "<b>Region of similarity:</b> $q1 &#150; $q2"
1207                    );
1208        $prot_query = ( 1.7 * abs( $q2 - $q1 ) < abs( $s2 - $s1 ) ) ? 1 : 0;
1209    
1210        push @genes, [ feature_graphic( $s1, $s2, $prot_query ),
1211                       'Q', undef, $info, undef, 'Query and match information'
1212                     ];
1213    
1214        return \@genes, $from, $to;
1215    }
1216    
1217    
1218    sub feature_graphic {
1219        my ( $beg, $end, $isprot ) = @_;
1220        my ( $min, $max, $symb, $color );
1221    
1222        ( $min, $max, $symb ) = ( $beg <= $end ) ? ( $beg, $end, "rightArrow" )
1223                                                 : ( $end, $beg, "leftArrow" );
1224    
1225        #  Color proteins by translation frame
1226    
1227        $color = $isprot ? qw( blue red green )[ $beg % 3 ] : 'grey';
1228    
1229        ( $min, $max, $symb, $color );
1230    }
1231    
1232    
1233    sub genes_in_region {
1234        my( $fig_or_sprout, $cgi, $genome, $contig, $min, $max ) = @_;
1235    
1236        if ( $cgi->param( 'SPROUT' ) )
1237        {
1238            my( $x, $feature_id );
1239            my( $feat, $min, $max ) = $fig_or_sprout->genes_in_region( $genome, $contig, $min, $max );
1240            my @tmp = sort { ($a->[1] cmp $b->[1]) or
1241                             (($a->[2]+$a->[3]) <=> ($b->[2]+$b->[3]))
1242                           }
1243                      map  { $feature_id = $_;
1244                             $x = feature_locationS( $fig_or_sprout, $feature_id );
1245                             $x ? [ $feature_id, boundaries_of( $fig_or_sprout, $x )]  : ()
1246                           }
1247                      @$feat;
1248            return ( [map { $_->[0] } @tmp ], $min, $max );
1249        }
1250        else
1251        {
1252            return $fig_or_sprout->genes_in_region( $genome, $contig, $min, $max );
1253        }
1254    }
1255    
1256    
1257    sub feature_locationS {
1258        my ( $fig_or_sprout, $peg ) = @_;
1259        scalar $fig_or_sprout->feature_location( $peg );
1260    }
1261    
1262    
1263    sub boundaries_of {
1264        my( $fig_or_sprout, $loc ) = @_;
1265        $fig_or_sprout->boundaries_of( $loc );
1266    }
1267    
1268    
1269    sub function_ofS {
1270        my( $fig_or_sprout, $peg, $user ) = @_;
1271        scalar $fig_or_sprout->function_of( $peg, $user );
1272    }
1273    
1274    
1275    sub feature_aliasesL {
1276        my( $fig_or_sprout, $fid ) = @_;
1277        my @tmp = $fig_or_sprout->feature_aliases( $fid );
1278        @tmp
1279    }
1280    
1281    
1282  sub format_sims {  sub format_sims {
1283      my($fig,$cgi,$html,$sims) = @_;      my($fig,$cgi,$html,$sims) = @_;
1284      my($col_hdrs,$table,@ids,$ids,$sim,%seen);      my($col_hdrs,$table,@ids,$ids,$sim,%seen);

Legend:
Removed from v.1.64  
changed lines
  Added in v.1.65

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3