ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/gclib/scripts/geanno_seqLoad
Revision: 24
Committed: Tue Jul 26 21:46:39 2011 UTC (8 years, 1 month ago) by gpertea
File size: 2763 byte(s)
Log Message:
Line File contents
1 #!/usr/bin/perl
2 use strict;
3 use Getopt::Std;
4 use FindBin;use lib $FindBin::Bin;
5 use dbSession;
6
7 my $usage = q/Usage:
8 geanno_seqLoad -t <taxon_id> -g <genome_release_name> <fasta_file(s)...>
9
10 Load the given fasta file(s) with the genomic sequences into the
11 geanno database, creating new entries in the ginfo and gseq tables
12 as needed.
13 For each chromosome\/contig loaded, ginfo.name and ginfo.accession
14 will be set as the sequence name found in the input fasta file(s).
15 /;
16 umask 0002;
17 getopts('t:g:c:d:') || die($usage."\n");
18 my ($taxonID, $genome_name, $descr, $srcdb)=
19 ($Getopt::Std::opt_t, $Getopt::Std::opt_g,
20 $Getopt::Std::opt_c, $Getopt::Std::opt_d);
21 die("$usage Error: -t and -g options are required!\n")
22 unless $taxonID && $genome_name;
23
24 my $ds=dbSession->new('geanno@NEOSYBASE');
25 my ($genus, $species);
26 if ($taxonID !~ m/^\d+$/) {
27 (my $t, $genus, $species)=$ds->getval("select taxon_id, genus, species from taxon where common_name='$taxonID'");
28 if ($t>0) {
29 $taxonID=$t;
30 }
31 else {
32 die("Error: cannot find taxon ID for organism '$taxonID'\n");
33 }
34 } #common name given, not numeric taxon ID
35 else {
36 ($genus, $species)=$ds->getval(qq/select genus, species
37 from taxon where taxon_id=$taxonID/);
38 die("Error: cannot locate taxon# $taxonID in the database!\n".
39 "Please create the entry first.\n") unless $genus;
40 }
41 print STDERR "Loading genomic sequence for taxon $taxonID : '$genus $species'\n";
42
43 my $chunksize=1024; # gseq.seq field length
44
45 my $seqname;
46 my $seqid=$ds->getval('select max(gseq_id) from ginfo');
47 $seqid=0 unless $seqid>0; #the very first time
48 my $seqchunk; #max 1024 seqchunk to load;
49 my $chunkid; #starting at 0 for each sequence <=> gseq.seg_id
50
51 while(<>) {
52 if (m/^>(\S+)/) { #defline
53 my $newseqname=$1;
54 if ($seqname) {
55 loadSeqChunk($seqchunk) if $seqchunk;
56 print STDERR " .. loaded $chunkid sequence chunks.\n";
57 }
58 $seqname=$newseqname;
59 $seqid++;
60 $chunkid=0;
61 $seqchunk='';
62 $ds->do(qq/insert ginfo (gseq_id, taxon_id, release, name)
63 values ($seqid, $taxonID, '$genome_name', '$seqname')/);
64 print STDERR "Loading sequence data for $seqname..\n";
65 next;
66 }
67 # sequence line
68 next unless $seqname;
69 tr/\r\n\t //d;
70 $seqchunk.=$_;
71 if ((my $lendif=length($seqchunk)-$chunksize)>=0) {
72 loadSeqChunk(substr($seqchunk,0,$chunksize)); #loads only the first $chunksize characters always!
73 $seqchunk = ($lendif>0) ? substr($seqchunk, $chunksize) : '';
74 }
75 }
76 if ($seqname) {
77 loadSeqChunk($seqchunk) if $seqchunk;
78 print STDERR " .. loaded $chunkid sequence chunks.\n";
79 }
80
81 $ds->logout();
82
83 sub loadSeqChunk {
84 my $seqdata=shift(@_);
85 $ds->do("insert gseq values ($seqid, $chunkid, '$seqdata')");
86 $chunkid++;
87 }

Properties

Name Value
svn:executable *