use Bio::Tools::Run::Phylo::Phast::PhastCons;

  # Make a PhastCons factory
  $factory = Bio::Tools::Run::Phylo::Phast::PhastCons->new();

  # Pass the factory an alignment and the corresponding species tree
  $align_filename = 't/data/apes.multi_fasta';
  $species_tree_filename = 't/data/apes.newick';
  @features = $factory->run($align_filename, $species_tree_filename);

  # or get a Bio::Align::AlignI (SimpleAlign) object from somewhere, and
  # generate the species tree automatically using a Bio::DB::Taxonomy database
  $tdb = Bio::DB::Taxonomy->new(-source => 'entrez');
  @features = $factory->run($aln_obj, $tdb);

  # @features is an array of Bio::SeqFeature::Annotated, one feature per
  # alignment sequence and prediction

 1. Use phyloFit to generate a tree model for initialization of the nonconserved
    model from the supplied alignment (all data) and species tree
 2. Run phastCons in 'training' mode for parameter estimation using all the
    alignment data and the model from step 1
 3. Run phastCons with the trees from step 2 to discover the most conserved
    regions

 1. Make sure the phastCons and phyloFit executables are in your path.
 2. Define an environmental variable PHASTDIR which is a 
    directory which contains the phastCons and phyloFit applications:
    In bash:

    export PHASTDIR=/home/username/phast/bin

    In csh/tcsh:

    setenv PHASTDIR /home/username/phast/bin

 3. Include a definition of an environmental variable PHASTDIR in
    every script that will use this PhastCons wrapper module, e.g.:

    BEGIN { $ENV{PHASTDIR} = '/home/username/phast/bin' }
    use Bio::Tools::Run::Phylo::Phast::PhastCons;

    return $PROGRAM_NAME;

    return $PROGRAM_DIR;

    my ($class, @args) = @_;
    my $self = $class->SUPER::new(@args);
    
    $self->_set_from_args(\@args, -methods => {(map { $_ => $PARAMS{$_} } keys %PARAMS),
                                               (map { $_ => $SWITCHES{$_} } keys %SWITCHES)},
                                  -create => 1);
    
    return $self;

    my ($self, $num) = @_;
    if (defined ($num)) {
        ($num > 0 && $num < 1) || $self->throw("target_coverage value must be between 0 and 1, exclusive");
        $self->{coverage} = $num;
    }
    return $self->{coverage} || return;

    my ($self, $num) = @_;
    if (defined ($num)) {
        ($num > 0 && $num < 1) || $self->throw("rho value must be between 0 and 1, exclusive");
        $self->{rho} = $num;
    }
    return $self->{rho} || return;

    my ($self, $aln, $tree) = @_;
    
    ($aln && $tree) || $self->throw("alignment and tree must be supplied");
    my $aln_obj = $self->_alignment($aln);
    $tree = $self->_tree($tree);
    
    # if aln was a file, set the alignment id to match file name
    if (-e $aln) {
        my $aln_id = basename($aln);
        ($aln_id) = $aln_id =~ /^([^\.]+)/;
        $aln_obj->id($aln_id);
    }
    
    return $self->_run;

    my $self = shift;
    
    my $exe = $self->executable || return;
    
    # use phyloFit to generate tree model initialization (?) using species tree
    # and alignment
    my $pf = Bio::Tools::Run::Phylo::Phast::PhyloFit->new(-verbose => $self->verbose, -quiet => $self->quiet);
    my $init_mod = $pf->run($self->_alignment, $self->_tree) || $self->throw("phyloFit failed to work as expected, is it installed?");
    
    # cd to a temp dir
    my $temp_dir = $self->tempdir;
    my $cwd = Cwd->cwd();
    chdir($temp_dir) || $self->throw("Couldn't change to temp dir '$temp_dir'");
    
    my $aln_file = $self->_write_alignment;
    
    # do training for parameter estimation
    my $command = $exe.$self->_setparams($aln_file, $init_mod);
    $self->debug("phastCons training command = $command\n");
    system($command) && $self->throw("phastCons training call ($command) crashed: $?");
    
    # do the final analysis
    $command = $exe.$self->_setparams($aln_file);
    $self->debug("phastCons command = $command\n");
    system($command) && $self->throw("phastCons call ($command) crashed: $?");
    
    # read in most_cons.bed as the result
    my $bedin = Bio::FeatureIO->new(-format => 'bed', -file => 'most_cons.bed');
    
    # cd back to orig dir
    chdir($cwd) || $self->throw("Couldn't change back to working directory '$cwd'");
    
    my @feats = ();
    my $aln = $self->_alignment;
    while (my $feat = $bedin->next_feature) {
        $feat->source_tag('phastCons');
        my $sv = Bio::Annotation::SimpleValue->new(-tagname => 'predicted', -value => 1);
        $feat->annotation->add_Annotation($sv);
        # $feat->type('TF_binding_site'); causes seg fault in subsequent clone()
        
        # features are in zero-based alignment coords; make a feature for each
        # alignment sequence
        foreach my $seq ($aln->each_seq) {
            my $clone = clone($feat);
            # $clone->type('TF_binding_site'); causes massive slowdown if you later store/retrieve these features from Bio::DB::SeqFeature database
            
            # give it the correct id
            $clone->seq_id($seq->id);
            
            # check and correct the coords (sequence may not have the feature)
            my $sloc = $seq->location_from_column($feat->start + 1) || next;
            my $eloc = $seq->location_from_column($feat->end + 1) || next;
            $clone->start($sloc->start - 1);
            $clone->end($eloc->end - 1);
            
            push(@feats, $clone);
        }
    }
    return @feats;

    my ($self, $aln_file, $init_mod) = @_;
    
    my $param_string = $self->SUPER::_setparams(-params => [keys %PARAMS],
                                                -switches => [keys %SWITCHES],
                                                -double_dash => 1,
                                                -underscore_to_dash => 1);
    
    $param_string .= ' --no-post-probs';
    my $aln_id = $self->_alignment->id;
    $param_string .= " --seqname $aln_id --idpref $aln_id" if $aln_id;
    $param_string .= ' --refidx 0';
    
    my $input = ' --msa-format FASTA '.$aln_file;
    if ($init_mod) {
        $param_string .= ' --estimate-trees mytrees '.$input.' '.$init_mod;
    }
    else {
        $param_string .= $input.' --most-conserved most_cons.bed --score mytrees.cons.mod,mytrees.noncons.mod';
    }
    
    return $param_string;
}

1;

  bioperl-l@bioperl.org                  - General discussion
  http://bioperl.org/wiki/Mailing_lists  - About the mailing lists

  http://redmine.open-bio.org/projects/bioperl/

 Title   : expected_length
 Usage   : $factory->expected_length(5);
 Function: Set transition probabilities such that the expected length of a
           conserved element is the supplied value. target_coverage() must also
           be set.
 Returns : int (default undef)
 Args    : None to get, int to set