[Bio] / Sprout / Sprout.pm Repository:
ViewVC logotype

Diff of /Sprout/Sprout.pm

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

revision 1.21, Fri Sep 9 20:40:41 2005 UTC revision 1.41, Tue Oct 18 02:24:23 2005 UTC
# Line 70  Line 70 
70    
71  * B<maxSequenceLength> maximum number of residues per sequence, (default C<8000>)  * B<maxSequenceLength> maximum number of residues per sequence, (default C<8000>)
72    
73    * B<noDBOpen> suppresses the connection to the database if TRUE, else FALSE
74    
75  =back  =back
76    
77  For example, the following constructor call specifies a database named I<Sprout> and a user name of  For example, the following constructor call specifies a database named I<Sprout> and a user name of
# Line 98  Line 100 
100                                                          # database connection port                                                          # database connection port
101                         maxSegmentLength => 4500,        # maximum feature segment length                         maxSegmentLength => 4500,        # maximum feature segment length
102                         maxSequenceLength => 8000,       # maximum contig sequence length                         maxSequenceLength => 8000,       # maximum contig sequence length
103                           noDBOpen     => 0,               # 1 to suppress the database open
104                        }, $options);                        }, $options);
105      # Get the data directory.      # Get the data directory.
106      my $dataDir = $optionTable->{dataDir};      my $dataDir = $optionTable->{dataDir};
# Line 105  Line 108 
108      $optionTable->{userData} =~ m!([^/]*)/(.*)$!;      $optionTable->{userData} =~ m!([^/]*)/(.*)$!;
109      my ($userName, $password) = ($1, $2);      my ($userName, $password) = ($1, $2);
110      # Connect to the database.      # Connect to the database.
111      my $dbh = DBKernel->new($optionTable->{dbType}, $dbName, $userName, $password, $optionTable->{port});      my $dbh;
112        if (! $optionTable->{noDBOpen}) {
113            $dbh = DBKernel->new($optionTable->{dbType}, $dbName, $userName,
114                                    $password, $optionTable->{port});
115        }
116      # Create the ERDB object.      # Create the ERDB object.
117      my $xmlFileName = "$optionTable->{xmlFileName}";      my $xmlFileName = "$optionTable->{xmlFileName}";
118      my $erdb = ERDB->new($dbh, $xmlFileName);      my $erdb = ERDB->new($dbh, $xmlFileName);
# Line 603  Line 610 
610          if ($prevContig eq $contigID && $dir eq $prevDir) {          if ($prevContig eq $contigID && $dir eq $prevDir) {
611              # Here the new segment is in the same direction on the same contig. Insure the              # Here the new segment is in the same direction on the same contig. Insure the
612              # new segment's beginning is next to the old segment's end.              # new segment's beginning is next to the old segment's end.
613              if (($dir eq "-" && $beg == $prevBeg - $prevLen) ||              if ($dir eq "-" && $beg + $len == $prevBeg) {
614                  ($dir eq "+" && $beg == $prevBeg + $prevLen)) {                  # Here we're merging two backward blocks, so we keep the new begin point
615                  # Here we need to merge two segments. Adjust the beginning and length values                  # and adjust the length.
616                  # to include both segments.                  $len += $prevLen;
617                    # Pop the old segment off. The new one will replace it later.
618                    pop @retVal;
619                } elsif ($dir eq "+" && $beg == $prevBeg + $prevLen) {
620                    # Here we need to merge two forward blocks. Adjust the beginning and
621                    # length values to include both segments.
622                  $beg = $prevBeg;                  $beg = $prevBeg;
623                  $len += $prevLen;                  $len += $prevLen;
624                  # Pop the old segment off. The new one will replace it later.                  # Pop the old segment off. The new one will replace it later.
# Line 615  Line 627 
627          }          }
628          # Remember this specifier for the adjacent-segment test the next time through.          # Remember this specifier for the adjacent-segment test the next time through.
629          ($prevContig, $prevBeg, $prevDir, $prevLen) = ($contigID, $beg, $dir, $len);          ($prevContig, $prevBeg, $prevDir, $prevLen) = ($contigID, $beg, $dir, $len);
630            # Compute the initial base pair.
631            my $start = ($dir eq "+" ? $beg : $beg + $len - 1);
632          # Add the specifier to the list.          # Add the specifier to the list.
633          push @retVal, "${contigID}_$beg$dir$len";          push @retVal, "${contigID}_$start$dir$len";
634      }      }
635      # Return the list in the format indicated by the context.      # Return the list in the format indicated by the context.
636      return (wantarray ? @retVal : join(',', @retVal));      return (wantarray ? @retVal : join(',', @retVal));
# Line 650  Line 664 
664      shift if UNIVERSAL::isa($_[0],__PACKAGE__);      shift if UNIVERSAL::isa($_[0],__PACKAGE__);
665      my ($location) = @_;      my ($location) = @_;
666      # Parse it into segments.      # Parse it into segments.
667      $location =~ /^(.*)_(\d*)([+-_])(\d*)$/;      $location =~ /^(.+)_(\d+)([+\-_])(\d+)$/;
668      my ($contigID, $start, $dir, $len) = ($1, $2, $3, $4);      my ($contigID, $start, $dir, $len) = ($1, $2, $3, $4);
669      # If the direction is an underscore, convert it to a + or -.      # If the direction is an underscore, convert it to a + or -.
670      if ($dir eq "_") {      if ($dir eq "_") {
# Line 758  Line 772 
772          # the start point is the ending. Note that in the latter case we must reverse the DNA string          # the start point is the ending. Note that in the latter case we must reverse the DNA string
773          # before putting it in the return value.          # before putting it in the return value.
774          my ($start, $stop);          my ($start, $stop);
775            Trace("Parse of \"$location\" is $beg$dir$len.") if T(SDNA => 4);
776          if ($dir eq "+") {          if ($dir eq "+") {
777              $start = $beg;              $start = $beg;
778              $stop = $beg + $len - 1;              $stop = $beg + $len - 1;
779          } else {          } else {
780              $start = $beg + $len + 1;              $start = $beg - $len + 1;
781              $stop = $beg;              $stop = $beg;
782          }          }
783            Trace("Looking for sequences containing $start through $stop.") if T(SDNA => 4);
784          my $query = $self->Get(['IsMadeUpOf','Sequence'],          my $query = $self->Get(['IsMadeUpOf','Sequence'],
785              "IsMadeUpOf(from-link) = ? AND IsMadeUpOf(start-position) + IsMadeUpOf(len) > ? AND " .              "IsMadeUpOf(from-link) = ? AND IsMadeUpOf(start-position) + IsMadeUpOf(len) > ? AND " .
786              " IsMadeUpOf(start-position) <= ? ORDER BY IsMadeUpOf(start-position)",              " IsMadeUpOf(start-position) <= ? ORDER BY IsMadeUpOf(start-position)",
# Line 776  Line 792 
792                  $sequence->Values(['IsMadeUpOf(start-position)', 'Sequence(sequence)',                  $sequence->Values(['IsMadeUpOf(start-position)', 'Sequence(sequence)',
793                                     'IsMadeUpOf(len)']);                                     'IsMadeUpOf(len)']);
794              my $stopPosition = $startPosition + $sequenceLength;              my $stopPosition = $startPosition + $sequenceLength;
795                Trace("Sequence is from $startPosition to $stopPosition.") if T(SDNA => 4);
796              # Figure out the start point and length of the relevant section.              # Figure out the start point and length of the relevant section.
797              my $pos1 = ($start < $startPosition ? 0 : $start - $startPosition);              my $pos1 = ($start < $startPosition ? 0 : $start - $startPosition);
798              my $len = ($stopPosition <= $stop ? $stopPosition : $stop) - $startPosition - $pos1;              my $len1 = ($stopPosition < $stop ? $stopPosition : $stop) + 1 - $startPosition - $pos1;
799                Trace("Position is $pos1 for length $len1.") if T(SDNA => 4);
800              # Add the relevant data to the location data.              # Add the relevant data to the location data.
801              $locationDNA .= substr($sequenceData, $pos1, $len);              $locationDNA .= substr($sequenceData, $pos1, $len1);
802          }          }
803          # Add this location's data to the return string. Note that we may need to reverse it.          # Add this location's data to the return string. Note that we may need to reverse it.
804          if ($dir eq '+') {          if ($dir eq '+') {
805              $retVal .= $locationDNA;              $retVal .= $locationDNA;
806          } else {          } else {
807              $locationDNA = join('', reverse split //, $locationDNA);              $retVal .= FIG::reverse_comp($locationDNA);
             $retVal .= $locationDNA;  
808          }          }
809      }      }
810      # Return the result.      # Return the result.
# Line 857  Line 874 
874      # Set it from the sequence data, if any.      # Set it from the sequence data, if any.
875      if ($sequence) {      if ($sequence) {
876          my ($start, $len) = $sequence->Values(['IsMadeUpOf(start-position)', 'IsMadeUpOf(len)']);          my ($start, $len) = $sequence->Values(['IsMadeUpOf(start-position)', 'IsMadeUpOf(len)']);
877          $retVal = $start + $len;          $retVal = $start + $len - 1;
878        }
879        # Return the result.
880        return $retVal;
881    }
882    
883    =head3 ClusterPEGs
884    
885    C<< my $clusteredList = $sprout->ClusterPEGs($sub, \@pegs); >>
886    
887    Cluster the PEGs in a list according to the cluster coding scheme of the specified
888    subsystem. In order for this to work properly, the subsystem object must have
889    been used recently to retrieve the PEGs using the B<get_pegs_from_cell> method.
890    This causes the cluster numbers to be pulled into the subsystem's color hash.
891    If a PEG is not found in the color hash, it will not appear in the output
892    sequence.
893    
894    =over 4
895    
896    =item sub
897    
898    Sprout subsystem object for the relevant subsystem, from the L</get_subsystem>
899    method.
900    
901    =item pegs
902    
903    Reference to the list of PEGs to be clustered.
904    
905    =item RETURN
906    
907    Returns a list of the PEGs, grouped into smaller lists by cluster number.
908    
909    =back
910    
911    =cut
912    #: Return Type $@@;
913    sub ClusterPEGs {
914        # Get the parameters.
915        my ($self, $sub, $pegs) = @_;
916        # Declare the return variable.
917        my $retVal = [];
918        # Loop through the PEGs, creating arrays for each cluster.
919        for my $pegID (@{$pegs}) {
920            my $clusterNumber = $sub->get_cluster_number($pegID);
921            # Only proceed if the PEG is in a cluster.
922            if ($clusterNumber >= 0) {
923                # Push this PEG onto the sub-list for the specified cluster number.
924                push @{$retVal->[$clusterNumber]}, $pegID;
925            }
926      }      }
927      # Return the result.      # Return the result.
928      return $retVal;      return $retVal;
# Line 1007  Line 1072 
1072    
1073  =head3 FeatureAnnotations  =head3 FeatureAnnotations
1074    
1075  C<< my @descriptors = $sprout->FeatureAnnotations($featureID); >>  C<< my @descriptors = $sprout->FeatureAnnotations($featureID, $rawFlag); >>
1076    
1077  Return the annotations of a feature.  Return the annotations of a feature.
1078    
# Line 1017  Line 1082 
1082    
1083  ID of the feature whose annotations are desired.  ID of the feature whose annotations are desired.
1084    
1085    =item rawFlag
1086    
1087    If TRUE, the annotation timestamps will be returned in raw form; otherwise, they
1088    will be returned in human-readable form.
1089    
1090  =item RETURN  =item RETURN
1091    
1092  Returns a list of annotation descriptors. Each descriptor is a hash with the following fields.  Returns a list of annotation descriptors. Each descriptor is a hash with the following fields.
1093    
1094  * B<featureID> ID of the relevant feature.  * B<featureID> ID of the relevant feature.
1095    
1096  * B<timeStamp> time the annotation was made, in user-friendly format.  * B<timeStamp> time the annotation was made.
1097    
1098  * B<user> ID of the user who made the annotation  * B<user> ID of the user who made the annotation
1099    
# Line 1035  Line 1105 
1105  #: Return Type @%;  #: Return Type @%;
1106  sub FeatureAnnotations {  sub FeatureAnnotations {
1107      # Get the parameters.      # Get the parameters.
1108      my ($self, $featureID) = @_;      my ($self, $featureID, $rawFlag) = @_;
1109      # Create a query to get the feature's annotations and the associated users.      # Create a query to get the feature's annotations and the associated users.
1110      my $query = $self->Get(['IsTargetOfAnnotation', 'Annotation', 'MadeAnnotation'],      my $query = $self->Get(['IsTargetOfAnnotation', 'Annotation', 'MadeAnnotation'],
1111                             "IsTargetOfAnnotation(from-link) = ?", [$featureID]);                             "IsTargetOfAnnotation(from-link) = ?", [$featureID]);
# Line 1048  Line 1118 
1118              $annotation->Values(['IsTargetOfAnnotation(from-link)',              $annotation->Values(['IsTargetOfAnnotation(from-link)',
1119                                   'Annotation(time)', 'MadeAnnotation(from-link)',                                   'Annotation(time)', 'MadeAnnotation(from-link)',
1120                                   'Annotation(annotation)']);                                   'Annotation(annotation)']);
1121            # Convert the time, if necessary.
1122            if (! $rawFlag) {
1123                $timeStamp = FriendlyTimestamp($timeStamp);
1124            }
1125          # Assemble them into a hash.          # Assemble them into a hash.
1126          my $annotationHash = { featureID => $featureID,          my $annotationHash = { featureID => $featureID,
1127                                 timeStamp => FriendlyTimestamp($timeStamp),                                 timeStamp => $timeStamp,
1128                                 user => $user, text => $text };                                 user => $user, text => $text };
1129          # Add it to the return list.          # Add it to the return list.
1130          push @retVal, $annotationHash;          push @retVal, $annotationHash;
# Line 1254  Line 1328 
1328          my $query = $self->Get(['IsBidirectionalBestHitOf'],          my $query = $self->Get(['IsBidirectionalBestHitOf'],
1329                                 "IsBidirectionalBestHitOf(from-link) = ? AND IsBidirectionalBestHitOf(genome) = ?",                                 "IsBidirectionalBestHitOf(from-link) = ? AND IsBidirectionalBestHitOf(genome) = ?",
1330                                 [$featureID, $genomeID]);                                 [$featureID, $genomeID]);
1331          # Look for the best hit.          # Peel off the BBHs found.
1332          my $bbh = $query->Fetch;          my @found = ();
1333          if ($bbh) {          while (my $bbh = $query->Fetch) {
1334              my ($targetFeature) = $bbh->Value('IsBidirectionalBestHitOf(to-link)');              push @found, $bbh->Value('IsBidirectionalBestHitOf(to-link)');
             $retVal{$featureID} = $targetFeature;  
1335          }          }
1336            $retVal{$featureID} = \@found;
1337      }      }
1338      # Return the mapping.      # Return the mapping.
1339      return \%retVal;      return \%retVal;
# Line 1538  Line 1612 
1612              Trace("Peg 1 is " . $peg1Data->[1] . " and Peg 2 is " . $peg2Data->[1] . ".") if T(Coupling => 4);              Trace("Peg 1 is " . $peg1Data->[1] . " and Peg 2 is " . $peg2Data->[1] . ".") if T(Coupling => 4);
1613              push @retVal, [$peg1Data->[1], $peg2Data->[1], $peg1Data->[0]];              push @retVal, [$peg1Data->[1], $peg2Data->[1], $peg1Data->[0]];
1614          }          }
1615            Trace("Last index in evidence result is is $#retVal.") if T(Coupling => 4);
1616      }      }
1617      # Return the result.      # Return the result.
1618      return @retVal;      return @retVal;
# Line 1587  Line 1662 
1662                                   [$retVal], ["ParticipatesInCoupling(from-link)", "Coupling(score)"]);                                   [$retVal], ["ParticipatesInCoupling(from-link)", "Coupling(score)"]);
1663      # Check to see if we found anything.      # Check to see if we found anything.
1664      if (!@pegs) {      if (!@pegs) {
1665            Trace("No coupling found.") if T(Coupling => 4);
1666          # No coupling, so undefine the return value.          # No coupling, so undefine the return value.
1667          $retVal = undef;          $retVal = undef;
1668      } else {      } else {
1669          # We have a coupling! Get the score and check for inversion.          # We have a coupling! Get the score and check for inversion.
1670          $score = $pegs[0]->[1];          $score = $pegs[0]->[1];
1671          $inverted = ($pegs[0]->[0] eq $peg1);          my $firstFound = $pegs[0]->[0];
1672            $inverted = ($firstFound ne $peg1);
1673            Trace("Coupling score is $score. First peg is $firstFound, peg 1 is $peg1.") if T(Coupling => 4);
1674      }      }
1675      # Return the result.      # Return the result.
1676      return ($retVal, $inverted, $score);      return ($retVal, $inverted, $score);
# Line 1697  Line 1775 
1775          if ($line =~ m/^>\s*(.+?)(\s|\n)/) {          if ($line =~ m/^>\s*(.+?)(\s|\n)/) {
1776              # Here we have a new header. Store the current sequence if we have one.              # Here we have a new header. Store the current sequence if we have one.
1777              if ($id) {              if ($id) {
1778                  $retVal{$id} = uc $sequence;                  $retVal{$id} = lc $sequence;
1779              }              }
1780              # Clear the sequence accumulator and save the new ID.              # Clear the sequence accumulator and save the new ID.
1781              ($id, $sequence) = ("$prefix$1", "");              ($id, $sequence) = ("$prefix$1", "");
1782          } else {          } else {
1783              # Here we have a data line, so we add it to the sequence accumulator.              # Here we have a data line, so we add it to the sequence accumulator.
1784              # First, we get the actual data out. Note that we normalize to upper              # First, we get the actual data out. Note that we normalize to lower
1785              # case.              # case.
1786              $line =~ /^\s*(.*?)(\s|\n)/;              $line =~ /^\s*(.*?)(\s|\n)/;
1787              $sequence .= $1;              $sequence .= $1;
# Line 1711  Line 1789 
1789      }      }
1790      # Flush out the last sequence (if any).      # Flush out the last sequence (if any).
1791      if ($sequence) {      if ($sequence) {
1792          $retVal{$id} = uc $sequence;          $retVal{$id} = lc $sequence;
1793      }      }
1794      # Close the file.      # Close the file.
1795      close FASTAFILE;      close FASTAFILE;
# Line 2037  Line 2115 
2115      # Get the parameters.      # Get the parameters.
2116      my ($self, $entityName, $entityID) = @_;      my ($self, $entityName, $entityID) = @_;
2117      # Check for the entity instance.      # Check for the entity instance.
2118        Trace("Checking existence of $entityName with ID=$entityID.") if T(4);
2119      my $testInstance = $self->GetEntity($entityName, $entityID);      my $testInstance = $self->GetEntity($entityName, $entityID);
2120      # Return an existence indicator.      # Return an existence indicator.
2121      my $retVal = ($testInstance ? 1 : 0);      my $retVal = ($testInstance ? 1 : 0);
# Line 3137  Line 3216 
3216    
3217  sub FriendlyTimestamp {  sub FriendlyTimestamp {
3218      my ($timeValue) = @_;      my ($timeValue) = @_;
3219      my $retVal = strftime("%a %b %e %H:%M:%S %Y", localtime($timeValue));      my $retVal = localtime($timeValue);
3220      return $retVal;      return $retVal;
3221  }  }
3222    

Legend:
Removed from v.1.21  
changed lines
  Added in v.1.41

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3