/usr/bin/bp_blast2tree is in bioperl 1.6.924-3.
This file is owned by root:root, with mode 0o755.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 | #!/usr/bin/perl
eval 'exec /usr/bin/perl -S $0 ${1+"$@"}'
if 0; # not running under some shell
# Author: Jason Stajich <jason@bioperl.org>
# Purpose: Blast Report -> MSA -> Tree
# This needs lots more error checking, cmdline input of parameters
# and support for other treebuilding -- only Phylip Neighbor for
# protein data is supported
# Also proper pulling in of the correct sequence data from the
# alignment - multiple hits on different parts of a protein aren't
# going to work properly right now. So this is mostly and example
# starting point which needs a lot more work to be made robust.
use strict;
use warnings;
use Bio::AlignIO;
use Bio::Tools::Run::Alignment::Clustalw;
use Bio::Tools::Run::Phylo::Phylip::ProtDist;
use Bio::Tools::Run::Phylo::Phylip::Neighbor;
use Bio::Tools::Run::Phylo::Molphy::ProtML;
use Bio::Tools::Run::Phylo::Phylip::ProtPars;
use Bio::SearchIO;
use Bio::SimpleAlign;
use Bio::PrimarySeq;
use Bio::TreeIO;
use Getopt::Long;
my $IDLENGTH = 12;
# we could in fact pull the tree out of the guide tree calculated
# by Clustalw in the alignment, but I believe that is UPGMA
# which would *NOT* be a good thing to be giving people.
my $aln_factory = Bio::Tools::Run::Alignment::Clustalw->new
('ktuple' => 2, "quiet" => 1,
'matrix' => 'BLOSUM');
my ($report,$format,$tree_method,$cutoff,$keepall);
$format = 'blast';
$tree_method = 'neighbor';
$cutoff = '0.01';
GetOptions(
'h|help' => sub { exec('perldoc', $0);
exit(0); },
'i|input:s' => \$report,
'f|format:s' => \$format,
'm|method:s' => \$tree_method,
'e|evalue:s' => \$cutoff,
'k|keepall:s' => \$keepall, # keep all HSPs not just best
);
unless( $format =~ /blast|fasta|hmmer/i ) {
die("Must request a valid search report format (fasta,blast,hmmer)");
}
unless ( $tree_method =~ /nj|neighbor/i || $tree_method =~ /protml|ml/i ) {
die("Must request a valid tree building method (neighbor,protml)");
}
my (@alns,@seqs);
my $in = new Bio::SearchIO(-file => $report,
-format => $format);
while( my $r = $in->next_result ) {
# Let's build the simplest system first
die("Cannot process report that does not contain protein sequence")
unless ($r->algorithm =~ /HMMER|BLASTP|FASTP/i );
my @seqs;
while( my $hit = $r->next_hit ) {
next if $hit->significance > $cutoff;
while( my $hsp = $hit->next_hsp ) {
next if $hsp->evalue > $cutoff;
my $seq = $hsp->get_aln->get_seq_by_pos(2)->seq();
push @seqs, new Bio::PrimarySeq(-seq => $seq,
-id => $hsp->hit->seq_id,
-desc => $r->algorithm . " best align to ". $hsp->query->seq_id );
last unless $keepall;
}
}
push @alns, $aln_factory->align(\@seqs);
}
my $out = new Bio::AlignIO(-format => 'phylip',
-interleaved => 1,
-idlength => $IDLENGTH,
-file => ">alignfile.phy");
$out->write_aln(@alns);
$out = undef;
# these need to be parameterized for cmdline arguments
my @params = ('idlength'=>$IDLENGTH,
'model'=>'cat',
'gencode'=>'U',
'category'=>'H',
'probchange'=>'0.4',
'trans'=>'5',
'freq'=>'0.25,0.5,0.125,0.125');
my $dist_factory = Bio::Tools::Run::Phylo::Phylip::ProtDist->new(@params);
$dist_factory->quiet(1);
@params = ('type'=>'NJ',
'outgroup'=>1,
'upptri'=>1,
'jumble'=>17);
my $tree_factory = Bio::Tools::Run::Phylo::Phylip::Neighbor->new(@params);
$tree_factory->quiet(1);
my $count = 1;
my $outtrees = new Bio::TreeIO(-file => ">trees.tre",
-format => 'newick');
foreach my $aln ( @alns ) {
# NOTE NOTE NOTE
# This interface is probably going to change from create_tree to
# next_tree per some discussions I'm having with Shawn - we may need
# to tweak any scripts before you publish
# also may move the create_distance_matrix method around some
# and need to write in the switched support for Molphy's ProtML
my $matrix = $dist_factory->create_distance_matrix($aln);
my @seqnames = keys %$matrix;
open my $MATRIX, '>', "Group$count.dist" or die "Could not write file 'Group$count.dist': $!\n";
printf $MATRIX "%4d\n",scalar @seqnames;
for(my $i =0; $i< (scalar @seqnames - 1); $i++ ) {
printf $MATRIX "%-12s ", $seqnames[$i];
for( my $j = $i+1; $j < scalar @seqnames; $j++ ) {
print $MATRIX $matrix->{$seqnames[$i]}->{$seqnames[$j]}," ";
}
print $MATRIX "\n";
}
close $MATRIX;
my $tree = $tree_factory->create_tree("Group$count.dist");
$outtrees->write_tree($tree);
$count++;
}
=head1 NAME
tree_from_seqsearch - builds a phylogenetic tree based on a sequence
search (FastA,BLAST,HMMER)
=head1 DESCRIPTION
This script requires that the bioperl-run pkg be also installed.
=cut
|