meme.pl

EM算法的改进

{
  if (!$address) {
    &whine("
      You must include a return e-mail address to receive your results.<BR>
      Please go back and include one.
    ");
  } else {
    $status = valid_address($address);
    if ($status == 0) {
      &whine("
	There is an error in your return email address:<BR>
	&nbsp&nbsp&nbsp&nbsp <TT>$address</TT><BR>
	It is possible that your email address is correct, in which case
	the problem may be that your host is behind a firewall and
	is consequently not found by the nslookup routines.  Consult with
	your systems people to see if you have an nslookup-visible email 
	address.  If none is available, please send email to <BR>
	&nbsp&nbsp&nbsp&nbsp <TT>meme\@nbcr.net</TT> <BR>
	mentioning this problem. 
      ");
    }
  }
} 

#
# check that number of sites is OK
#
sub check_nsites {
  my(
    $minsites, 			# minimum nsites
    $maxsites, 			# maximum nsites
    $dist, 			# type of distribution
    $num			# number of sequences
  ) = @_;

  if ($minsites && $minsites < $MINSITES) {
    whine("You must specify a minimum number of sites >= $MINSITES");
  }
  if ($maxsites && $maxsites > $MAXSITES) {
    whine("You must specify a maximum number of sites <= $MAXSITES");
  }
  if (($minsites && $maxsites) && $minsites > $maxsites) {
    whine(
     "The minimum number of sites is larger than the maximum number of sites");
  }
} # check_nsites

#
# check that width is OK
#
sub check_width
{
  if ($minw < $MINW) {
    whine("The minimum width you specified ($minw) is too small.<BR>")
  }
  if ($maxw > $MAXW) {
    whine("The maximum width you specified ($maxw) is too large.<BR>")
  }
  if ($minw > $maxw) {
    whine("The minimum width you specified ($minw) is larger than the
	maximum width you specified ($maxw).<BR>")
  }
} # check_width

#
# Check to see whether the data are valid.
#
sub check_data
{
  # Convert the data to FASTA format.
  $status = system 
    ("$bin/readseq -a -f=8 $datafile 1>meme.data.$$ 2>meme.error.$$");
  $error = `cat meme.error.$$`; unlink "meme.error.$$";
  if ($status) {
    &whine("An error occurred when the READSEQ program attempted to convert
      your dataset to FASTA format.<BR>
      READSEQ returned: $error
    ");
   return(1);
  }

  # Run the 'getsize' program to get information on the sequence data. 
  $status = system 
    ("$bin/getsize $datafile 1>getsize.$$ 2>getsize.error.$$");
  unlink "getsize.error.$$";
  chop($getsize1 = `cat getsize.$$`); unlink "getsize.$$";
  ($num, $min, $max, $ave, $total, $letters) = split (' ', "$getsize1");

  # Run the 'getsize' program to get information on converted data; will
  # be unchanged if MEME knows the format.
  $status = system 
    ("$bin/getsize -nd meme.data.$$ 1>getsize.$$ 2>getsize.error.$$");
  $error = `cat getsize.error.$$`; unlink "getsize.error.$$";
  if ($error || $status) {
    &whine("After converting to FASTA format using the READSEQ program,
      the following errors in your dataset were detected:<BR>$error
      <BR>Make sure all your sequences are in the same format since READSEQ
      assumes that all sequences are in the same format as the first sequence.
    ");
    unlink "getsize.$$";
    return(1);
  }
  chop($getsize2 = `cat getsize.$$`); unlink "getsize.$$";
  ($num, $min, $max, $ave, $total, $letters) = split (' ', "$getsize2");

  # Use original dataset if MEME recognizes it
  if ($getsize1 eq $getsize2) {
    #print "<P>Using original dataset<P>";
    $error = `cp $datafile meme.data.$$`;
    if ($error ne "") {
      &whine("An error occurred while trying to copy your data.<BR>
        cp returned: $error
      ");
    }
  }

  # check for problem reading the dataset
  if ($num <= 0) {
    &whine("Your dataset appears to be in a format that MEME does not recognize.
      <BR> Please check to be sure that your data is 
      <A HREF=../help_format.html> formatted</A> properly.
    ");
    &copy_stdout ("$datafile");
  }

  # check for bad sequences
  if ($num > 0 && $min == 0) {
    &whine("Your dataset appears to contain one or more zero-length sequences.
      <BR> Please check to be sure that your data is 
      <A HREF=../help_format.html> formatted</A> properly.
    ");
  }

  # get the alphabet used in the sequences
  $alphabet = get_alph($letters);
} # check_data

#
# get the alphabet: DNA or PROTEIN
#
sub get_alph 
{
  local ($letters) = @_;		# get arguments

  $_ = $letters;
  $old = length; 

  # check against allowed dna letters
  $x = $_;
  $x =~ tr/ABCDGHKMNRSTUVWY//cd;
  $new = length $x;
  if ($old == $new) {
    "dna";
  } else {
    # check against allowed protein letters
    $x = $_;
    $x =~ tr/ABCDEFGHIKLMNPQRSTUVWXYZ//cd;
    $new = length $x;
    if ($old == $new) {
      "protein";
    } else {
      # get the unknown letters
      $x = $_;
      $x =~ tr/ABCDEFGHIKLMNPQRSTUVWXYZ//d;
      &whine("
	Your sequences contained the following unrecognized letters: $x.
	<BR>
	Please convert your sequences to one of the sequence
	<A HREF=../help_alphabet.html>alphabets</A> which MEME recognizes.
      ");
    }
  }
} # get_alph

#
# Submit job to meme client 
#
sub submit
{
  # submit it
  if (!$debug) {
    $status = system ("$bin/$MEME_CLIENT meme.msg.$$ 1>meme.client.msg.$$ 2>&1");
    $status /= 256;
  }

  # check for errors
  if ($status) {

    open (ERROR, "meme.client.msg.$$");
    $_ = <ERROR>;
    while (<ERROR>) {$error = $_}               # get last line of message
    if ($status == 12) {$msg = `cat ../meme-downtimes`;}
    &whine("An error occured while submitting your job: status = $status. <BR>
      Error message: $error <BR>
      $msg
      Please try again later.");
    close (ERROR);
    # warn maintainer
    $wd = `pwd`; chop $wd;
    $command = 
      "cat meme.client.msg.$$ | \
        @mail@ -s \'$MEME_CLIENT error: $status  user: $address  dir: $wd\' $maint";
    # warn maintainer unless debugging
    if (!$debug) { system($command); }

  } else {			# print verification form
    print "
      Your MEME search results will be sent to: <b> $address </b><br>
      If you do not receive a confirming email message, there could be an
      error in your email address.
      <HR>
      <UL>
      <LI> E-mail address:<B> $address</B>
      <LI> Sequence file:<B> $datafile_name </B>
      <LI> Description:<B> $subject</B>
      <LI> Distribution of motif occurrences:<B>
    ";
    if ($dist eq "oops") {
      print  "One per sequence</B>\n";
    } elsif ($dist eq "zoops") {
      print  "Zero or one per sequence</B>\n";
    } else {
      print  "Any number of repetitions</B>\n";
    }
    print "
      <LI> Number of different motifs:<B> $nmotifs </B>
      <LI> Minimum number of sites:<B> $minsites</B>
      <LI> Maximum number of sites:<B> $maxsites</B>
      <LI> Minimum motif width:<B> $minw</B>
      <LI> Maximum motif width:<B> $maxw</B>
    ";
    if ($posonly) {
      print "<LI> Searching given strand only";
    }
    if ($pal) {
      print "<LI> Looking for palindromes only";
    }
    if ($text) {
      print "<LI> Output is text (not HTML)";
    }
    if ($shuffle) {
      print "<LI> Shuffling letters in input sequences.";
    }
    if ($bfile) {
      print "<LI> Background model file: $bfile";
    }
    print "
      <LI> Statistics on your dataset:
        <TABLE BORDER>
	  <TR> <TD> type of sequence <TH ALIGN=RIGHT> $alphabet
	  <TR> <TD> number of sequences <TH ALIGN=RIGHT> $num
	  <TR> <TD> shortest sequence (residues) <TH ALIGN=RIGHT> $min 
	  <TR> <TD> longest sequence (residues) <TH ALIGN=RIGHT> $max 
	  <TR> <TD> average sequence length (residues) <TH ALIGN=RIGHT> $ave 
	  <TR> <TD> total dataset size (residues) <TH ALIGN=RIGHT> $total 
        </TABLE>
      </UL>
    ";
  }
  # print the message if debug
  if ($debug) {&copy_stdout ("meme.msg.$$");}
} # submit

#
# Print the output file to be sent to the meme-client.
#
sub print_file
{
  # First put a header on the file.
  open(outfile, "> meme.msg.$$") || whine("Can't open meme.msg.$$"); 
  print outfile <<END;
ADDRESS $address
PROGRAM meme
REMOTENAME $datafile_name
DESCRIPTION $subject $viewer
SWITCHES $switches
LOGINFO	n=$num N=$total shuffle=$shuffle
BEGIN
END
  close (outfile);

  # Then append the sequence data to the message file
  if ($shuffle) {
    $error = `$bin/fasta-shuffle-letters -tod < meme.data.$$ >> meme.msg.$$`;
  } else {
    $error = `cat meme.data.$$ >> meme.msg.$$`;
  }
} # print_file

#
# print an error message, bump the error count and continue
#
sub whine
{
  if ($nerrors == 0) {
    print "
      <H1>Error Report:</H1>
      <HR>
      <OL>
      <B>
    ";
  }
  print "<LI><B>@_<B>\n";
  $nerrors++;
} # whine

#
# copy a file to standard output
#
sub copy_stdout
{
  open (F, "@_");
  print "<PRE>";
  while (<F>) {
    print $_;
  }
  print "</PRE>";
} #copy_stdout