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

Annotation of /FigKernelScripts/get_close_pegs_kmer_memcache.pl

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : olson 1.1 use strict;
2 :     use Data::Dumper;
3 :     use gjoseqlib;
4 :     use IO::Compress::Gzip;
5 :     use Getopt::Long;
6 :     use Cache::Memcached::Fast;
7 :    
8 :     my $dataD = "/vol/ross/GeneralKmerCompareRegions/Data";
9 :    
10 : olson 1.2 open(SIZE, "<", "$dataD/hash.size") || die "Cannot open $dataD/hash.size: $!";
11 : olson 1.1 my $hash_size = <SIZE>;
12 :     chomp $hash_size;
13 :     close(SIZE);
14 :    
15 :     our @doctypes = qw(peg
16 :     rna
17 :     atn
18 :     att
19 :     bs
20 :     opr
21 :     pbs
22 :     pi
23 :     pp
24 :     prm
25 :     pseudo
26 :     rsw
27 :     sRNA
28 :     trm
29 :     box
30 :     );
31 :    
32 :     $| = 1;
33 :    
34 :     our %tmap;
35 :     for my $i (0..$#doctypes)
36 :     {
37 :     $tmap{$doctypes[$i]} = $i;
38 :     $tmap{$i} = $doctypes[$i];
39 :     }
40 :    
41 :     my $cache = new Cache::Memcached::Fast({
42 :     servers => [ { address => 'elm.mcs.anl.gov:11211' } ],
43 :     compress_threshold => 10_000,
44 :     compress_ratio => 0.9,
45 :     compress_methods => [ \&IO::Compress::Gzip::gzip,
46 :     \&IO::Uncompress::Gunzip::gunzip ],
47 :     max_failures => 3,
48 :     failure_timeout => 2,
49 :     });
50 :    
51 :     while (my($fid, $def, $seq) = read_next_fasta_seq())
52 :     {
53 :     $seq = uc $seq;
54 :     my $kmers = &get_kmers($seq);
55 :     my %hits;
56 :     my $res = $cache->get_multi(@$kmers);
57 :     while (my($kmer, $pegsN) = each(%$res))
58 :     {
59 :     foreach my $pegN (split(/,/,$pegsN))
60 :     {
61 :     $hits{$pegN}++;
62 :     }
63 :     }
64 :     my @sorted = map { [$_,$hits{$_}] } sort { $hits{$b} <=> $hits{$a} } keys(%hits);
65 :    
66 :     my @output;
67 :     foreach my $tuple (@sorted)
68 :     {
69 :     my($pegE,$n) = @$tuple;
70 :     my $peg = docid_to_fid($pegE);
71 :    
72 :     push(@output,[$pegE,$peg,$n]);
73 :     }
74 :     foreach $_ (@output)
75 :     {
76 :     print join("\t",@$_),"\n";
77 :     }
78 :     print "//\n";
79 :     }
80 :    
81 :     sub get_kmers {
82 :     my($seq) = @_;
83 :    
84 :     open(TMP,">tmp.$$") || die "could not open tmp.$$";
85 :     print TMP ">peg\n$seq\n";
86 :     close(TMP);
87 :    
88 :     my @hits = map { ($_ =~ /^HIT\s+\d+\s+(\d+)/) ? $1 : () }
89 :     `kmer_guts -d 1 -D $dataD/Data.Kmers -s $hash_size -a < tmp.$$`;
90 :     unlink("tmp.$$");
91 :     my %hits = map { ($_ => 1) } @hits;
92 :     return [sort { $a <=> $b } keys(%hits)];
93 :     }
94 :    
95 :    
96 :     sub fid_to_docid
97 :     {
98 :     my($fid) = @_;
99 :    
100 :     if ($fid =~ /^fig\|(\d+)\.(\d+)\.([^.]+)\.(\d+)$/)
101 :     {
102 :     my ($g, $ext, $type, $num) = ($1, $2, $3, $4);
103 :     my $tnum = $tmap{$type};
104 :    
105 :     #
106 :     # right to left: (cumulative)
107 :     # 17 bits for feature number (0)
108 :     # 4 bits for type (17)
109 :     # 15 bits for ext (21)
110 :     # Rest for genome (36)
111 :    
112 :     my $enc;
113 :    
114 :     $enc = $g << 36| $ext << 21 | $tnum << 17 | $num;
115 :    
116 :     return $enc;
117 :     }
118 :    
119 :     return undef;
120 :     }
121 :    
122 :     sub docid_to_fid
123 :     {
124 :     my($doc) = @_;
125 :    
126 :     my($g, $e, $t, $n);
127 :    
128 :     $g = $doc >> 36;
129 :     $e = ($doc >> 21) & 0x7fff;
130 :     $t = ($doc >> 17) & 0xf;
131 :     $n = $doc & 0x1ffff;
132 :    
133 :     my $type = $tmap{$t};
134 :     my $genome = "$g.$e";
135 :     my $fid = "fig|$genome.$type.$n";
136 :    
137 :     return $fid;
138 :     }
139 :    

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3