view release on metacpan or  search on metacpan

lib/Bio/RNA/RNAaliSplit/WrapRscape.pm  view on Meta::CPAN

  my $count = 0;
  open my $stk, "<", $self->ifile or croak "ERROR [$this_function] Cannot open Stockholm file $self->ifile for reading";
  while (<$stk>){
    next if (/^#/);
    next if (/^\/\//);
    $count++;
  }
  close ($stk);
  $self->cseq($count);
}

sub run_rscape {
  my $self = shift;
  my $this_function = (caller(0))[3];
  my ($out_fn,$sout_fn,$out,$sout,$sum);
  my ($rscape_out,$rscape_sout,$rscape_sum);
  my $tag = "";
  if ($self->has_statistic){$tag = ".".$self->statistic};

  if ($self->has_basename){
    $out_fn  = $self->basename.$tag."."."rscape.out";
    $sout_fn = $self->basename.$tag."."."rscape.sorted.out";
  }
  elsif ($self->has_ifilebn){
    $out_fn  = $self->ifilebn.$tag."."."rscape.out";
    $sout_fn = $self->ifilebn.$tag."."."rscape.sorted.out";
  }
  else{
    $out_fn  = $tag."rscape.out";
    $sout_fn = $tag."rscape.sorted.out";
  }
  $out  = file($oodir,$out_fn);  # R-scape stdout
  $sout = file($oodir,$sout_fn); # R-scape sorted stdout

  $rscape_out = "rscape.out";
  $rscape_sout = $rscape_out.".sorted";

  my $rscape_options = " -o $rscape_out --rna --outdir $oodir ";
  if ($self->has_nofigures && $self->nofigures == 1){$rscape_options.=" --nofigures "};
  if ($self->has_statistic){$rscape_options.=" --".$self->statistic." "  }
  my $cmd = $rscape.$rscape_options.$self->ifile;

  my ( $success, $error_message, $full_buf, $stdout_buf, $stderr_buf ) =
    run( command => $cmd, verbose => 0 );
  if( !$success ) {
    print STDERR "ERROR [$this_function] Call to $rscape unsuccessful\n";
    print STDERR "ERROR: $cmd\n";
    print STDERR "ERROR: this is what the command printed:\n";
    print join "", @$full_buf;
    croak $!;
  }

  $self->_parse_rscape($rscape_out);

  rename $rscape_out, $out;
  rename $rscape_sout, $sout;
}

# parse R-scape output
# canonical R-scape output is expected to look like this:
#  MSA myaln nseq 9 (9) alen 45 (45) avgid 79.51 (79.51) nbpairs 11 (11)
#  contacts  11 (11 bpairs 11 wc bpairs)
#  maxD      8.00
#  mind      1
#  Method Target_E-val [cov_min,conv_max] [FP | TP True Found | Sen PPV F]
#  RAFS    0.05           [-1.00,0.83]    [0 | 3 11 3 | 27.27 100.00 42.86]
#        left_pos       right_pos        score   E-value
# ------------------------------------------------------------
#*	        20	        36	0.83	0.000328218
#*	        21	        35	0.83	0.000328218
#*	        22	        34	0.33	0.0497359

sub _parse_rscape {
  my ($self,$out) = @_;
  my $this_function = (caller(0))[3];
  my @buffer = ();
  my $parse1 = 0;
  my $parse2 = 0;
  my $nosbp = 0;
  open my $file, "<", $out or croak "ERROR: [$this_function] Cannot open file $out";
  while(<$file>){
    chomp;
    if (m/^#\s+MSA\s+([a-zA-Z0-9_.|]+)\s+nseq\s+(\d+)\s+\(\d+\)\s+alen\s+(\d+)\s+\(\d+\)\s+avgid\s+(\d+\.\d+)\s+\(\d+\.\d+\)\s+nbpairs\s+(\d+)\s+\(\d+\)/g){
      $self->nseq($2);
      $self->alen($3);
      $self->nbpairs($5);
      $parse1 = 1;
      next;
    }
    if (m/^#\s+([a-zA-Z0-9]+)\s+(\d+\.\d+)\s+\[(\-?\d+\.\d+),(-?\d+\.\d+)\]\s+\[(\d+)\s+\|\s+(\d+)\s+(\d+)\s+(\d+)\s+\|\s+(\d+\.\d+)\s+(\d+\.\d+)\s+(\d+\.\d+)\]/g){
      $self->evalue($2);
      $self->FP($5);
      $self->TP($6);
      $self->T($7);
      $self->F($8);
      $self->Sen($9);
      $self->PPV($10);
      $self->Fmeasure($11);
      $parse2 = 1;
      next;
    }
    if (m/^\*\s+(\d+)\s+(\d+)\s+(\d+\.\d+)\s+(\d+\.\d+)$/){
      my %bp = (i => $1,
		j => $2,
		score => $3,
		evalue => $4);
      push @{$self->sigBP}, \%bp;
    }
    if (m/^no significant pairs$/){
      $nosbp=1;
    }
  }
  close ($file);
  #carp "INFO: [$this_function] parse1:".eval($parse1);
  #carp "INFO: [$this_function] parse2:".eval($parse2);

  if ($nosbp == 1){
    $self->status(1); # no significant basepairs
    $self->_check_attributes();
    return;
  }
  if ($parse1 == 1 && $parse2 == 1){
    $self->_check_attributes();
    $self->status(0); # all OK
    return;
  }
  elsif ($parse1 == 1 && $parse2 == 0){
    $self->status(2); # covariation scores are almost constant, no further analysis
    return;
  }
  else{
    croak "ERROR: [$this_function] ambiguous status when parsing rscape output file. parse1:".eval($parse1)." parse2:".eval($parse2)." This shouldn't happen ...";
    $self->status(3);
    return;
  }
}

sub _check_attributes {
  my ($self,$out) = @_;
  my $this_function = (caller(0))[3];

  $self->has_cseq     ? 1 : croak "ERROR [$this_function] \$self->cseq not set ". $self->ifile;
  $self->has_nseq     ? 1 : croak "ERROR [$this_function] \$self->nseq not set ". $self->ifile;