← Back to branch summary

~ubuntu-branches/ubuntu/hoary/bioperl/hoary

« back to all changes in this revision

Viewing changes to examples/run_genscan.pl

  • Committer: Bazaar Package Importer
  • Author(s): Matt Hope
  • Date: 2004-04-18 14:24:11 UTC
  • mfrom: (1.2.1 upstream) (2.1.1 warty)
  • Revision ID: james.westby@ubuntu.com-20040418142411-gr92uexquw4w8liq
Tags: 1.4-1
* New upstream release
* Examples and working code are installed by default to usr/bin,
  this has been moved to usr/share/doc/bioperl/bin

Show diffs side-by-side

added added

removed removed

Lines of Context:
examples/run_genscan.pl
1
 
#!/usr/bin/perl -w
2
 
# Brian Osborne
3
 
# script to run genscan on all nucleotide sequences in a fasta file
4
 
# and save results as fasta, creates <file>.gs.pept and <file>.gs.cds
5
 
 
6
 
use Bio::SeqIO;
7
 
use Getopt::Long;
8
 
use Bio::Tools::Genscan;
9
 
use Carp;
10
 
use strict;
11
 
 
12
 
# directory with GENSCAN matrices
13
 
my $genscanDir = "/dbs/genscan";
14
 
# GENSCAN location
15
 
my $binDir = "/usr/local/bin";
16
 
# GENSCAN matrix
17
 
my $matrix = "HumanIso.smat";
18
 
 
19
 
my ($file,$in);
20
 
 
21
 
GetOptions( "f|file=s" => \$file );
22
 
usage() if ( !$file );
23
 
 
24
 
# create output files for predicted protein and CDS
25
 
open PEPT, ">>$file.gs.pept" or croak "Error opening $file.gs.pept: $!\n";
26
 
open DNA, ">>$file.gs.cds" or croak "Error opening $file.gs.cds: $!\n";
27
 
 
28
 
$in = Bio::SeqIO->new(-file => $file , -format => 'Fasta');
29
 
 
30
 
while ( my $seq = $in->next_seq() ) {
31
 
    croak "Input sequence is protein\n" if ( $seq->moltype eq 'protein' );
32
 
    my $str = $seq->seq;
33
 
    my $id = $seq->id;
34
 
    $id = $1 if ( $id =~ /gi\|(\d+)/ );
35
 
    # create temp file, input to GENSCAN
36
 
    open OUT,">$id.temp.fa" or croak "Error opening $id.temp.fa: $!\n";
37
 
    print OUT ">$id.fa\n$str\n\n";
38
 
 
39
 
    # the contents of the *raw file will be parsed
40
 
    system "genscan $genscanDir/$matrix $id.temp.fa -cds > $id.gs.raw";
41
 
    unlink "$id.temp.fa";
42
 
    my $genscan = Bio::Tools::Genscan->new( -file => "$id.gs.raw");
43
 
    while ( my $gene = $genscan->next_prediction() ) {
44
 
        my $prt = $gene->predicted_protein;
45
 
        my $cds = $gene->predicted_cds;
46
 
 
47
 
        if ( defined $cds  ) {
48
 
            $cds->display_id =~ /predicted_(CDS_\d+)/;
49
 
            print DNA ">" . $id . "_" . $1 . " " . $cds->display_id . "\n"
50
 
              . $cds->seq . "\n";
51
 
        }
52
 
        if ( defined $prt ) {
53
 
            $prt->display_id =~ /predicted_(peptide_\d+)/;
54
 
            print PEPT ">" . $id . "_" . $1 . " " . $prt->display_id . "\n"
55
 
              . $prt->seq . "\n";
56
 
        }
57
 
    }
58
 
    $genscan->close();
59
 
    unlink "$id.gs.raw";
60
 
}
61
 
 
62
 
 
63
 
sub usage {
64
 
    print "
65
 
Usage    : $0 -f <file>
66
 
Function : run genscan on all nucleotide sequences in a multiple fasta file
67
 
Output   : <file>.gs.pept and <file>.gs.cds
68
 
 
69
 
";
70
 
    exit;
71
 
}