[Bio] / FigKernelScripts / index_feature_dir_fasta.pl Repository:
ViewVC logotype

Annotation of /FigKernelScripts/index_feature_dir_fasta.pl

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.1 - (view) (download) (as text)

1 : olson 1.1 #
2 :     # Normalize the given feature directory per blastall rules and index with formatdb.
3 :     #
4 :     # Create btree index of the tbl data as well.
5 :     #
6 :    
7 :     use FIG;
8 :     use FIG_Config;
9 :     use strict;
10 :     use DB_File;
11 :     use Data::Dumper;
12 :    
13 :     @ARGV == 1 or die "Usage: $0 feature-directory\n";
14 :    
15 :     my $dir = shift;
16 :    
17 :     -d $dir or die "Given directory $dir is not a directory\n";
18 :    
19 :     my $fasta = "$dir/fasta";
20 :     my $tbl = "$dir/tbl";
21 :    
22 :     open(IN, "<", $fasta) or die "Cannot open $fasta: $!";
23 :     open(TBL, "<", $tbl) or die "Cannot open $tbl: $!";
24 :    
25 :     #
26 :     # Normalize & create blast db.
27 :     #
28 :    
29 :     my $norm = "$fasta.norm";
30 :     open(OUT, ">", $norm) or die "Cannot open $norm for writing: $!";
31 :    
32 :     while (<IN>)
33 :     {
34 :     s/^>fig\|/>gnl|fig|/;
35 :     print OUT;
36 :     }
37 :     close(IN);
38 :     close(OUT);
39 :    
40 :     my @cmd = ("$FIG_Config::ext_bin/formatdb", "-o", "T", "-i", $norm);
41 :     my $rc = system(@cmd);
42 :     if ($rc != 0)
43 :     {
44 :     die "Error $rc running formatdb: @cmd\n";
45 :     }
46 :    
47 :     #
48 :     # Create btree index of tbl file.
49 :     #
50 :     # $idx{$fid} = join($;, contig, loc, index, aliases)
51 :     # $idx{$contig} = join($;, index-start, index-stop)
52 :     #
53 :     # The second entry records the beginning and ending indices
54 :     # of the given contig in the recno list of pegs in contig order.
55 :     #
56 :     # We also create a recno keyed on index as above, in order to
57 :     # accelerate genes_in_region.
58 :     #
59 :     # $list[order] = join($;, $fid, $contig, $beg, $end, $strand)
60 :     #
61 :    
62 :     my $btree = "$tbl.btree";
63 :     my $recno = "$tbl.recno";
64 :     my %idx;
65 :     my @list;
66 :    
67 :     if (-f $btree)
68 :     {
69 :     unlink($btree);
70 :     }
71 :    
72 :     if (-f $recno)
73 :     {
74 :     unlink($recno);
75 :     }
76 :    
77 :     my $t = tie %idx, 'DB_File', $btree, O_RDWR | O_CREAT, 0666, $DB_BTREE;
78 :     $t or die "Cannot tie $btree: $!";
79 :    
80 :     my $r = tie @list, 'DB_File', $recno, O_RDWR | O_CREAT, 0666, $DB_RECNO;
81 :     $r or die "Cannot tie $recno: $!";
82 :    
83 :     my @tmp;
84 :     while (<TBL>)
85 :     {
86 :     chomp;
87 :     my($id, $loc, @aliases) = split(/\t/);
88 :    
89 :     my($contig, $beg, $end, $strand) = &FIG::boundaries_of($loc);
90 :     push(@tmp, [$id, $contig, $beg, $end, $strand, $loc, \@aliases]);
91 :     }
92 :    
93 :     my $idx = 0;
94 :     my $last_contig;
95 :     my $last_contig_start;
96 :     for my $ent (sort { $a->[1] cmp $b->[1] or $a->[2] <=> $b->[2] } @tmp)
97 :     {
98 :     my($id, $contig, $beg, $end, $strand, $loc, $alist) = @$ent;
99 :    
100 :     if ($contig ne $last_contig)
101 :     {
102 :     if (defined($last_contig))
103 :     {
104 :     $idx{$last_contig} = join($;, $last_contig_start, $idx - 1);
105 :     }
106 :     $last_contig = $contig;
107 :     $last_contig_start = $idx;
108 :     }
109 :     $list[$idx] = join($;, $id, $contig, $beg, $end, $strand);
110 :     $idx{$id} = join($;, $loc, $idx, @$alist);
111 :     $idx++;
112 :     }
113 :     if (defined($last_contig))
114 :     {
115 :     $idx{$last_contig} = join($;, $last_contig_start, $idx - 1);
116 :     }
117 :    
118 :     untie $r;
119 :     untie $t;

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3