Skip to content

Instantly share code, notes, and snippets.

/Perldoc Secret

Created January 15, 2012 12:37
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save anonymous/5e57f3e9c9041c7615f1 to your computer and use it in GitHub Desktop.
Save anonymous/5e57f3e9c9041c7615f1 to your computer and use it in GitHub Desktop.
use strict;
use Bio::SearchIO;
use Bio::SeqIO;
use Bio::Seq;
my $input = new Bio::SeqIO (-file => 'lotus_ver25.fa',
-format=>'fasta');
my $in = new Bio::SearchIO(-format => 'blast',
-file => 'Search65.txt');
while( my $result = $in->next_result ) {
## $result is a Bio::Search::Result::ResultI compliant object
while( my $hit = $result->next_hit ) {
## $hit is a Bio::Search::Hit::HitI compliant object
while( my $hsp = $hit->next_hsp ) {
## $hsp is a Bio::Search::HSP::HSPI compliant object
## Let's discard some sequences
if( $hsp->length('total') > 50 ) {
if ( $hsp->percent_identity >= 65 ) {
if ( $hsp->evalue<=1e-25)
{
## Here we +- the start end
my $begin =$hsp->start('hit');
my $ennd =$hsp->end('hit');
my $begin1 = $begin-4000;
my $ennd1 = $ennd+4000;
## Here we read the genome file, and extract the sequence
while( my $seqobj = $input->next_seq ) {
while( my $substring = $seqobj->subseq($begin1,$ennd1)) {
while( my $obj->desc($result)){
open FILE, ">>SEQUENCES.txt";
select FILE;
print ">", $obj->query_name,
"\n",
$substring, "\n";
close FILE
}}}
}
}
}
}
}
}
exit 0;
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment