CMsearchWrapper

#!/usr/bin/perl -w -- -*-Perl-*-


##############################################################################
#
#                                 CMsearchWrapper
#
# DESCRIPTION:
# Simple wrapper for CMsearch (Infernal package) in order to create AnnotPairCollection
# for mfannot.
#
##############################################################################

#############################################################################
#                                 HMMsearchWrapper                          #
#                                                                           #
#                          Copyright (C) 2008                               #
#                         Departement de Biochimie,                         #
#                          Universite de Montreal,                          #
#                     C.P. 6128, succursale Centre-ville,                   #
#                      Montreal, Quebec, Canada, H3C 2J7                    #
#                                                                           #
#                Programming:  Natacha Beck, Pierre Rioux.                  #
#                Project management: Franz Lang (OGMP)                      #
#                E-Mail information: Franz.Lang@Umontreal.ca                #
#                                                                           #
#     This software is distributed under the GNU GENERAL PUBLIC LICENSE, as #
# published by the Free Software Foundation. A copy of version 2 of this    #
# license should be included in a file called COPYING. If not, write to the #
# Free Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.   #
#############################################################################

##########################
# Initialization section #
##########################

require 5.00;
use strict;
use vars qw( $VERSION $RCS_VERSION $DEBUG $TMPDIR );
use File::Basename;
use File::Path;
use PirObject;

PirObject->LoadDataModel("CMsearchOutput");
PirObject->LoadDataModel("AnnotPairCollection");

# Default umask
umask 027;
# Default path
my $PATH = $ENV{"PATH"} || "";
my @PATH = split(/:/, $PATH);


# Program's name and version number.
$RCS_VERSION='$Id: CMsearchWrapper,v 1.4 2011/12/28 16:14:54 nbeck Exp $';
($VERSION) = ($RCS_VERSION =~ m#,v ([\w\.]+)#);
my ($BASENAME) = ($0 =~ /([^\/]+)$/);

##################################
# Global variables and constants #
##################################

$|=1;
$DEBUG=0;
my $TMPDIR="/tmp/$BASENAME.$$";
if (! -d $TMPDIR) {
    mkdir($TMPDIR,0700) or die "Error: can't create work directory '$TMPDIR': $!\n";
}

# Basic command-line options
my $GENOME="";     # Genome file
my $CM_MODEL="";   # CM model file
my $NAME="";       # Gene name
my $TBL_OUT="";    # TBL output file of cmsearch
my $ALI_OUT="";    # ALI output file of cmsearch
my $APC_OUT="";    # AnnotPairCollection output file
my $LOG_OUT="";    # Log output file
my $COMMENT="";    # Comment for mfannot
my $EVALUE=0.0001; # E-value threshold for CMsearch
my $CMsearch_path= &GetPath("cmsearch"); # Path to cmsearch

#####################################
# Command-Line Arguments Processing #
#####################################

for (;@ARGV;) {
    my ($opt,$arg) = ($ARGV[0] =~ /^-([dgnmtaolce])(.*)$/o);
    last if ! defined $opt;
    # See for option... with argumentcm
    if ($opt =~ /[gnmtaolce]/ && $arg eq "") {
        if (@ARGV < 2) {
            print STDERR "Argument required for switch \"-$opt\".\n";
            exit 1;
        }
        shift @ARGV;
        $arg=$ARGV[0];
    }

    #     CMsearchWrapper -c "test" -o %TMPDIR%/tmRNA/tmRNA.xml -g %PLAINFASTAFILE% -m %MODPATH%/CM_models/tmRNA.cm -n %GENENAME%
    $DEBUG     = 1                         if $opt eq 'd';
    $GENOME    = $arg                      if $opt eq 'g';
    $NAME      = $arg                      if $opt eq 'n';
    $CM_MODEL  = $arg                      if $opt eq 'm';
    $TBL_OUT   = $arg                      if $opt eq 't';
    $ALI_OUT   = $arg                      if $opt eq 'a';
    $APC_OUT   = $arg                      if $opt eq 'o';
    $LOG_OUT   = $arg                      if $opt eq 'l';
    $EVALUE    = $arg                      if $opt eq 'e';
    $COMMENT   = $arg                      if $opt eq 'c';

    shift @ARGV;
}

sub Usage {
    my $message = shift || "";
    print STDERR <<USAGE;
Basic Usage: $BASENAME

where:
  -d          : debug mode
  -g <file>   : Genome file
  -n <string> : Gene name
  -c <string> : Comment for mfannot
  -m <file>   : CM model file
  -a <file>   : ALI output file of cmsearch
  -t <file>   : TBL output file of cmsearch
  -l <file>   : Log file a human readable log file
  -o <file>   : AnnotPairCollection output file (XML format)
  -e <float>  : E-value threshold for CMsearch (default 0.0001)

USAGE
    print STDERR $message,"\n" if $message;
    exit 20;
}


###########################################
# Search for program                      #
###########################################

sub GetPath {
  my $name_prog = shift;

  foreach my $dir (@PATH) {
    if (-f "$dir/$name_prog") {
      if (-r _ && -x _) {
        return "$dir/$name_prog";
      }
      else {
        die "   -> ERROR: $name_prog is not readable and executable! Please run:\n",
            "             chmod 755 \"$dir/ $name_prog\"\n";
      }
    last;
    }
  }
  die "-> ERROR: Could not find '$name_prog' in your search path. Please install\n",
      "   $name_prog from the source (see INSTALL.txt).\n";
}
###########################################
# Validate remaining command-line options #
###########################################

# Check genome file
&Usage("Error: Missing genome file (-g option)") if $GENOME eq "";
&Usage("Error: the genome file '$GENOME' does not exist") if ! -f $GENOME;

# Check model file
&Usage("Error: Missing CM model file (-m option)") if $CM_MODEL eq "";
&Usage("Error: the CM model file '$CM_MODEL' does not exist") if ! -f $CM_MODEL;

# Check output files
&Usage("Error: File '$TBL_OUT' give for option -t is not writable by you\n") if (-f ($TBL_OUT) && (!(-w ($TBL_OUT))));
&Usage("Error: File '$APC_OUT' give for option -o is not writable by you\n") if (-f ($APC_OUT) && (!(-w ($APC_OUT))));
&Usage("Error: File '$ALI_OUT' give for option -a is not writable by you\n") if (-f ($ALI_OUT) && (!(-w ($ALI_OUT))));
&Usage("Error: File '$LOG_OUT' give for option -l is not writable by you\n") if (-f ($LOG_OUT) && (!(-w ($LOG_OUT))));

# Check other options
&Usage("Error: Missing gene name (-n option)") if $NAME eq "";

############################################
# Main program                             #
############################################

my ($SEQ,$ID2SEQ)       = &CleanFastaAndChangeHeader($GENOME);
my ($RES_FILE,$RES_ALI_FILE) = &RunCMsearch($SEQ,$CM_MODEL);
my ($APC)         = &CreateAPC($RES_FILE,$ID2SEQ);
print "Temporary work directory $TMPDIR cleaned up ...\n" if $DEBUG;
if ($APC_OUT) {
  print "APC (.xml) file: $APC_OUT\n";
  system("cp $APC $APC_OUT");
}

if ($TBL_OUT) {
    print "TBL file: $TBL_OUT";
    system("cp $RES_FILE $TBL_OUT");
    print "\n";
}

if ($ALI_OUT) {
    print "ALI file: $ALI_OUT";
    system("cp $RES_ALI_FILE $ALI_OUT");
    print "\n";
}

if ($LOG_OUT) {
    &CreateLogFile($LOG_OUT,$RES_FILE,$RES_ALI_FILE);
}

END {
    # With exit, programme will go here
    # Cleanup temp directory when program exits.
    return unless defined($TMPDIR) and $TMPDIR =~ m#^/tmp/#;
    print "Temporary work directory $TMPDIR NOT cleaned up ...\n" if $DEBUG;
    rmtree($TMPDIR) unless $DEBUG;
}

#####################################
# ReverseFasta                      #
#####################################

sub CleanFastaAndChangeHeader {
    my $file = shift;

    # First step read file
    my $fh_infile  = new IO::File "<$file"
        or die "Can't read from file '$file': $!\n";

    my ($cnt,$id2seq) = ({},{});
    my $id      = "";
    my $tab_seq = [];
    while (my $line = <$fh_infile>) {
        if ($line =~ m#^>(.*)#) {
            if ($id ne "") {
                push(@$tab_seq,$id2seq);
                $id2seq = {};
            }
            $id = "$1";
            $cnt->{$id}++;
            if ($cnt->{$id} > 1) {
                print "2 seq have same name $id\n";
            }
            $id2seq->{Header} = "$id";
            next;
        }
        $id2seq->{Seq} .= "$line" if $id;
    }
    push(@$tab_seq,$id2seq); # last seq
    $fh_infile->close();

    # Second Create new File
    my $Seq     = "$TMPDIR/Seq";
    my $fh_Seq = new IO::File ">$Seq"
        or die "Can't write in file '$Seq': $!\n";

    $id2seq = {};
    my $num_header = 0;
    foreach my $entry (@$tab_seq)  {
        my $new_header = "cg_$num_header";
        my $name       = $entry->{Header};
        my $seq        = $entry->{Seq};
        my $seq_wo_bl = $seq;
           $seq_wo_bl =~ s/\s+//g;

        $id2seq->{$new_header}->{Seq} = $seq;
        $id2seq->{$new_header}->{Length} = length($seq_wo_bl);
        $id2seq->{$new_header}->{OriName} = $name;
        print $fh_Seq    ">$new_header\n$seq\n\n\n";
        $num_header++;
    }

    $fh_Seq->close();

    return ($Seq,$id2seq);
}

#####################################
# RunCMsearch                       #
#####################################

sub RunCMsearch {
    my ($seq,$cm_model) = @_;

    # Run on + strand
    my $tbl_output = "$TMPDIR/CMres.tbl";
    my $ali_output = "$TMPDIR/CMres.ali";
    my $cmd    = "nice -19 $CMsearch_path -E $EVALUE -A $ali_output --tblout $tbl_output $cm_model $seq >/dev/null 2>/dev/null";
    print "$cmd\n" if $DEBUG;
    my $resultat  = system($cmd);
    &PrintErrorSystem($resultat, "Error when proccessing CMsearch\n");

    return ($tbl_output,$ali_output);
}

#####################################
# Create AnnotPairCollection    #
#####################################


# Based on the output of CMsearch,
# create an AnnotPairCollection
sub CreateAPC {
    my ($infile,$id2seq) = @_;
    my $comment = $COMMENT;

    my $CMsearchItems = &ReadCMsearch($infile);

    my $outfile = "$TMPDIR/CMres.apc";
    my $fh_outfile = new IO::File ">$outfile"
      or die "Can't write in file '$outfile': $!\n";

    my $APs = [];
    my $APC_by_contig = {};
    foreach my $CMsearchItem (@$CMsearchItems) {
        my $contigname = $CMsearchItem->get_TargetName();
        my $isMinus    = $CMsearchItem->get_Strand() eq "-";
        my $dir        = !$isMinus ? "==>" : "<==";
        my $evalue     = $CMsearchItem->get_Evalue();

        my $startline = ";     G-$NAME $dir start;; mfannot: $comment, CM evalue: $evalue";
        my $endline   = ";     G-$NAME $dir end";

        my $AP = new PirObject::AnnotPair(
            type      => "G",
            genename  => $NAME,
            score     => $CMsearchItem->get_Evalue(),
            startpos  => $CMsearchItem->get_SequenceFrom(),
            endpos    => $CMsearchItem->get_SequenceTo(),
            direction => $CMsearchItem->get_Strand() eq "+" ? "==>" : "<==",
            startline => $startline,
            endline   => $endline,
            idbyHMM   => 1,
            );
        $APC_by_contig->{$contigname} = [] if !$APC_by_contig->{$contigname};
        push(@{$APC_by_contig->{$contigname}},$AP);
    }

    foreach my $contigname (keys %$APC_by_contig) {
        my $APC = new PirObject::AnnotPairCollection();
        $APC->set_genename($NAME);
        $APC->set_contigname($ID2SEQ->{$contigname}->{"OriName"});
        $APC->set_annotpairlist($APC_by_contig->{$contigname});
        my $xml_APC = $APC->ObjectToXML();
        print $fh_outfile $xml_APC;
    }

    $fh_outfile->close();
    print "APC created in $outfile\n";
    return ($outfile);
}

sub CreateLogFile {
    my ($outfile,$tbl_file,$ali_file) = @_;

    my $CMsearchItems = &ReadCMsearch($RES_FILE);
    my $nb_sol = scalar(@$CMsearchItems);

    my $comment = $COMMENT || "";
    # Remove the 'note=' in comment if present
    $comment =~ s/note=//g;
    my $identified = "$NAME";
    $identified .= " ($comment)" if $comment;

    my $short_summary = "";

    $short_summary .= "\nIdentified RNA structure with CM model for $identified\n";
    $short_summary .= sprintf("\n%-30s %15s \n%-30s %15s", "","Retained #","","of solution");
    $short_summary .= sprintf("\n%-30.30s %15s",$identified,$nb_sol);

    $short_summary  .= "\n\n******************************************************************************\n";
    $short_summary  .= "List of identified $identified prediction :\n\n";

    $short_summary .=  "                                                RNA    RNA\n";
    $short_summary .=  "Species name                    RNA # Strand  Begin    End             E-value\n";
    $short_summary .=  "------------------------------ ------ ------ ------ ------ ---- ----- --------\n";

    # Iteshort_summaryover all CMsearchItems
    my $nb = 1;
    foreach my $CMsearchItem (@$CMsearchItems) {
        # Add one line to the short summary for each CMsearchItem
        my $contigname = $CMsearchItem->get_TargetName();
        my $strand     = $CMsearchItem->get_Strand() eq "+" ? "FW" : "RC";
        my $start      = $CMsearchItem->get_SequenceFrom();
        my $end        = $CMsearchItem->get_SequenceTo();
        my $evalue     = $CMsearchItem->get_Evalue();

        $short_summary .= sprintf("%-30.30s %6s %6s %6s %6s %4s %5s %8s\n",
                                  $ID2SEQ->{$contigname}->{"OriName"},
                                  $nb,
                                  $strand,
                                  $start,
                                  $end,
                                  "",
                                  "",
                                  $evalue);

        $nb = $nb + 1;
    }

    $short_summary .= "\n\n*************************************************\n";
    $short_summary .= "Condensed alignment for $NAME\n";
    $short_summary .= "*************************************************\n";
    $short_summary .= "$nb_sol match(es)\n";
    $short_summary .= "\n";

    # Cat the content of the alignment file
    my $io_infile  = new IO::File "<$ali_file"
        or die "Can't read from file '$ali_file': $!\n";
    my @tab = <$io_infile>;
    $io_infile->close();
    $short_summary .= join("",@tab);

    # Save the short summary in the log file
    my $fh_outfile = new IO::File ">$outfile"
        or die "Can't write in file '$outfile': $!\n";

    print $fh_outfile $short_summary;
    $fh_outfile->close();

    print "Log file created in $outfile\n";
}

sub ReadCMsearch {
    my $file = shift;

    my $io_infile  = new IO::File "<$file"
        or die "Can't read from file '$file': $!\n";
    my @tab = <$io_infile>;
    $io_infile->close();

    my $CMsearchOutput = new PirObject::CMsearchOutput();
    my $CMRes          = $CMsearchOutput->FillCMsearchItems(\@tab);
}

#####################################
# PrintErrorSystem                  #
#####################################

sub PrintErrorSystem {
    my ($result,$message) = @_;

    my $hascoredump = ($result & 128) >> 7;  # 0 if no core dump, 1 if core dump
    my $signal      = $result & 127;  # SIGNAL received by subprocess, from 0 to 127;
    my $returncode  = $result >> 8;   # exit status of subprogram
    if ($returncode > 1 || $signal > 0 || $hascoredump == 1) {
        print "$message";
    }
}