← Back to branch summary

~ubuntu-branches/ubuntu/raring/bioperl/raring

« back to all changes in this revision

Viewing changes to examples/tools/extract_genes.pl

  • Committer: Bazaar Package Importer
  • Author(s): Charles Plessy
  • Date: 2008-03-18 14:44:57 UTC
  • mfrom: (4 hardy)
  • mto: This revision was merged to the branch mainline in revision 6.
  • Revision ID: james.westby@ubuntu.com-20080318144457-1jjoztrvqwf0gruk
http://bugs.debian.org/448890
* debian/control:
  - Removed MIA Matt Hope (dopey) from the Uploaders field.
    Thank you for your work, Matt. I hope you are doing well.
  - Downgraded some recommended package to the 'Suggests' priority,
    according to the following discussion on Upstream's mail list.
    http://bioperl.org/pipermail/bioperl-l/2008-March/027379.html
    (Closes: #448890)
* debian/copyright converted to machine-readable format.

Show diffs side-by-side

added added

removed removed

Lines of Context:
examples/tools/extract_genes.pl
 
1
#!/usr/bin/perl -w
 
2
# $Id: extract_genes.pl,v 1.4 2006/02/28 02:52:11 bosborne Exp $
 
3
=pod
 
4
 
 
5
=head1 NAME
 
6
 
 
7
extract_genes.pl - extract genomic sequences from NCBI files
 
8
using BioPerl
 
9
 
 
10
=head1 DESCRIPTION
 
11
 
 
12
This script is a simple solution to the problem of
 
13
extracting genomic regions corresponding to genes. There are other 
 
14
solutions, this particular approach uses genomic sequence 
 
15
files from NCBI and gene coordinates from Entrez Gene.
 
16
 
 
17
The first time this script is run it will be slow as it will
 
18
extract species-specific data from the gene2accession file and create
 
19
a storable hash (retrieving the positional data from this hash is
 
20
significantly faster than reading gene2accession each time the script
 
21
runs). The subsequent runs should be fast.
 
22
 
 
23
=head1 INSTALLATION
 
24
 
 
25
=head2
 
26
 
 
27
Install BioPerl, full instructions at http://bioperl.org.
 
28
 
 
29
=head2 Download gene2accession.gz
 
30
 
 
31
Download this file from ftp://ftp.ncbi.nlm.nih.gov/gene/DATA into 
 
32
your working directory and gunzip it.
 
33
 
 
34
=head2 Download sequence files
 
35
 
 
36
Create one or more species directories in the working directory, the
 
37
directory names do not have to match those at NCBI (e.g. "Sc", "Hs").
 
38
 
 
39
Download the nucleotide fasta files for a given species from its CHR*
 
40
directories at ftp://ftp.ncbi.nlm.nih.gov/genomes and put these files into a 
 
41
species directory. The sequence files will have the suffix ".fna" or 
 
42
"fa.gz", gunzip if necessary.
 
43
 
 
44
=head2 Determine Taxon id
 
45
 
 
46
Determine the taxon id for the given species. This id is the first column
 
47
in the gene2accession file. Modify the %species hash in this script
 
48
such that name of your species directory is a key and the taxon id is the 
 
49
value.
 
50
 
 
51
=head2 Command-line options
 
52
 
 
53
  -i   Gene id
 
54
  -s   Name of species directory
 
55
  -h   Help
 
56
 
 
57
Example:
 
58
 
 
59
  extract_genes.pl -i 850302 -s Sc
 
60
 
 
61
=cut
 
62
 
 
63
use strict;
 
64
use Bio::DB::Fasta;
 
65
use Getopt::Long;
 
66
use Storable;
 
67
 
 
68
my %species = ( "Sc" => 4932,  # Saccharomyces cerevisiae
 
69
                                     "Ec" => 83333, # Escherichia coli K12
 
70
                                          "Hs" => 9606   # H. sapiens
 
71
                                   );
 
72
 
 
73
my ($help,$id,$name);
 
74
 
 
75
GetOptions( "s=s"  =>  \$name,
 
76
            "i=i"  =>  \$id,
 
77
                                "h"    =>  \$help );
 
78
 
 
79
usage() if ($help || !$id || !$name);
 
80
 
 
81
my $storedHash = $name . ".dump";
 
82
 
 
83
# create index for a directory of fasta files
 
84
my $db = Bio::DB::Fasta->new($name, -makeid => \&make_my_id);
 
85
 
 
86
# extract species-specific data from gene2accession
 
87
unless (-e $storedHash) {
 
88
        my $ref;
 
89
        # extract species-specific information from gene2accession
 
90
        open MYIN,"gene2accession" or die "No gene2accession file\n";
 
91
        while (<MYIN>) {
 
92
                my @arr = split "\t",$_;
 
93
                if ($arr[0] == $species{$name} && $arr[9] =~ /\d+/ && $arr[10] =~ /\d+/) {
 
94
                        ($ref->{$arr[1]}->{"start"}, $ref->{$arr[1]}->{"end"}, 
 
95
                         $ref->{$arr[1]}->{"strand"}, $ref->{$arr[1]}->{"id"}) =        
 
96
                                ($arr[9], $arr[10], $arr[11], $arr[7]);
 
97
                }
 
98
        }
 
99
        # save species-specific information using Storable
 
100
        store $ref, $storedHash;
 
101
} 
 
102
 
 
103
# retrieve the species-specific data from a stored hash
 
104
my $ref = retrieve($storedHash);
 
105
 
 
106
# retrieve sequence and sub-sequence
 
107
if (defined $ref->{$id}) {
 
108
        my $chr = $db->get_Seq_by_id($ref->{$id}->{"id"});
 
109
        my $seq = $chr->trunc($ref->{$id}->{"start"},$ref->{$id}->{"end"});
 
110
        $seq = $seq->revcom if ($ref->{$id}->{"strand"} eq "-");
 
111
 
 
112
        # Insert SeqIO options here...
 
113
        print $seq->seq,"\n";
 
114
} else {
 
115
        print "Cannot find id: $id\n";
 
116
}
 
117
 
 
118
sub make_my_id {
 
119
        my $line = shift;
 
120
        $line =~ /ref\|([^|]+)/;
 
121
        $1;
 
122
}
 
123
 
 
124
sub usage {
 
125
        system "perldoc $0";
 
126
        exit;
 
127
}
 
128
 
 
129
__END__