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

Annotation of /Kmers2/mk-kmers.pl

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.1 - (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 :     pareach \@work, sub {
59 :     my $ent = shift;
60 :     my($tag, $file, $output) = @$ent;
61 :    
62 :     open(my $fh, "<", $file) or die "Cannot open $file; $!";
63 :     my $base = basename($file);
64 :     open(my $out, "| sort --parallel=1 -S 20G -T $sort_tmp -u > $output") or die "Cannot open $output: $!";
65 :     # open(my $out, ">", $output) or die "Cannot open $output: $!";
66 :    
67 :     while (my($id, $def, $seq) = read_next_fasta_seq($fh))
68 :     {
69 :     my $s = length($seq);
70 :     for (my $i=0; ($i < (length($seq) - $k)); $i++)
71 :     {
72 :     my $kmer = uc substr($seq,$i,$k);
73 :    
74 :     if ($kmer !~ /[^ACDEFGHIKLMNPQRSTVWY]/)
75 :     {
76 :     my $rev = reverse( $kmer );
77 :     $rev =~ tr/ACGTUMRWSYKBDHV/TGCAAKYWSRMVHDB/;
78 :    
79 :     print $out ($kmer lt $rev ? $kmer : $rev), "\t", $tag, "\n";
80 :     }
81 :     }
82 :     }
83 :    
84 :     close($fh);
85 :     }, { Max_Workers => $nprocs };
86 :    
87 :     #
88 :     # We have the raw data. Now sort/unique to prepare for merge.
89 :     #
90 :    
91 :     pareach \@dirs, sub {
92 :     my($dir) = shift;
93 :    
94 :     my @files = <$part_dir/$dir/tmp.$k/*>;
95 :     print "$dir: @files\n";
96 :    
97 :     my $cmd = "sort --parallel=1 -m -T $sort_tmp -S 20G -u @files > $part_dir/$dir/raw.$k";
98 :     print "$cmd\n";
99 :     my $rc = system($cmd);
100 :     if ($rc != 0)
101 :     {
102 :     die "failed with rc=$rc: $cmd\n";
103 :     }
104 :     }, { Max_Workers => $nprocs };
105 :    
106 :     #
107 :     # Now we can create the signatures.
108 :     #
109 :    
110 :     my @raws = map { "$part_dir/$_/raw.$k" } @dirs;
111 :    
112 :     my $cmd = "sort -S 20G -T $sort_tmp -m @raws | uniq -u -w $k > $out_kdir/sigs";
113 :    
114 :     print "$cmd\n";
115 :     my $rc = system($cmd);
116 :     $rc == 0 or die "failed with rc=$rc: $cmd\n";

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3