[Bio] / Kmers2 / mk-kmers.pl Repository:
ViewVC logotype

Annotation of /Kmers2/mk-kmers.pl

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : olson 1.1 use Data::Dumper;
2 :     use File::Basename;
3 :     use lib '../FigKernelPackages';
4 :     use strict;
5 :     use Proc::ParallelLoop;
6 :     use gjoseqlib;
7 :    
8 :     =head1 NAME
9 :    
10 :     mk-kmers
11 :    
12 :     =head1 SYNOPSIS
13 :    
14 :     mk-kmers partition-dir K
15 :    
16 :     =head1 DESCRIPTION
17 :    
18 :     Compute the kmer signatures for the given partition dir.
19 :    
20 :     =cut
21 :    
22 :     @ARGV == 2 or die "Usage: mk-kmers partition-dir k\n";
23 :    
24 :     my $part_dir = shift;
25 :     my $k = shift;
26 :    
27 :     my $sort_tmp = "/space/partitions/tmp";
28 :    
29 :     -d $part_dir or die "Directoy $part_dir does not exist\n";
30 :     $k =~ /^\d+$/ or die "Invalid K '$k'\n";
31 :    
32 :     my $out_kdir = "$part_dir/sigs.$k";
33 :     -d $out_kdir or mkdir($out_kdir) or die "Cannot mkdir $out_kdir: $!\n";
34 :    
35 :     my $nprocs = 6;
36 :    
37 :     opendir(D, $part_dir);
38 :     my @dirs = grep { -d "$part_dir/$_" && /set.\d+/ } readdir(D);
39 :    
40 :     my @work;
41 :     for my $dir (@dirs)
42 :     {
43 :     my($id) = $dir =~ /set\.(\d+)/;
44 :     open(F, "<", "$part_dir/$dir/definition") or die "Cannot open $dir/definition: $!";
45 :    
46 :     my $tmp = "$part_dir/$dir/tmp.$k";
47 :     -d $tmp or mkdir $tmp or die "Cannot mkdir $tmp: $!";
48 :    
49 :     $_ = <F>;
50 :     my @files = <F>;
51 :     chomp @files;
52 :    
53 :     push @work, map { my $base = basename($_); [$id, $_, "$tmp/$base"] } @files;
54 :     }
55 :    
56 :     print scalar(@work) . " work items scheduled for $nprocs workers\n";
57 :    
58 : olson 1.2 #@work = grep { $_->[0] == 3} @work;
59 :     @work = ();
60 :    
61 : olson 1.1 pareach \@work, sub {
62 :     my $ent = shift;
63 :     my($tag, $file, $output) = @$ent;
64 :    
65 :     open(my $fh, "<", $file) or die "Cannot open $file; $!";
66 :     my $base = basename($file);
67 :     open(my $out, "| sort --parallel=1 -S 20G -T $sort_tmp -u > $output") or die "Cannot open $output: $!";
68 :     # open(my $out, ">", $output) or die "Cannot open $output: $!";
69 :    
70 :     while (my($id, $def, $seq) = read_next_fasta_seq($fh))
71 :     {
72 :     my $s = length($seq);
73 :     for (my $i=0; ($i < (length($seq) - $k)); $i++)
74 :     {
75 :     my $kmer = uc substr($seq,$i,$k);
76 :    
77 :     if ($kmer !~ /[^ACDEFGHIKLMNPQRSTVWY]/)
78 :     {
79 :     my $rev = reverse( $kmer );
80 :     $rev =~ tr/ACGTUMRWSYKBDHV/TGCAAKYWSRMVHDB/;
81 :    
82 :     print $out ($kmer lt $rev ? $kmer : $rev), "\t", $tag, "\n";
83 :     }
84 :     }
85 :     }
86 :    
87 :     close($fh);
88 :     }, { Max_Workers => $nprocs };
89 :    
90 :     #
91 :     # We have the raw data. Now sort/unique to prepare for merge.
92 :     #
93 :    
94 :     pareach \@dirs, sub {
95 :     my($dir) = shift;
96 :    
97 :     my @files = <$part_dir/$dir/tmp.$k/*>;
98 :     print "$dir: @files\n";
99 :    
100 :     my $cmd = "sort --parallel=1 -m -T $sort_tmp -S 20G -u @files > $part_dir/$dir/raw.$k";
101 :     print "$cmd\n";
102 :     my $rc = system($cmd);
103 :     if ($rc != 0)
104 :     {
105 :     die "failed with rc=$rc: $cmd\n";
106 :     }
107 :     }, { Max_Workers => $nprocs };
108 :    
109 :     #
110 :     # Now we can create the signatures.
111 :     #
112 :    
113 :     my @raws = map { "$part_dir/$_/raw.$k" } @dirs;
114 :    
115 :     my $cmd = "sort -S 20G -T $sort_tmp -m @raws | uniq -u -w $k > $out_kdir/sigs";
116 :    
117 :     print "$cmd\n";
118 :     my $rc = system($cmd);
119 :     $rc == 0 or die "failed with rc=$rc: $cmd\n";

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3