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

Annotation of /FigKernelScripts/CSA_make_unique_kmers.pl

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : overbeek 1.1 use strict;
2 :     use Data::Dumper;
3 :     use gjoseqlib;
4 :     use SeedAware;
5 :     use SeedEnv;
6 :    
7 :     my $N;
8 :     my $rpts;
9 :     my $usage = "usage: CSA_make_unique_kmers N RepeatKmers < Contigs > Kmers";
10 :     (
11 :     ($N = shift @ARGV) &&
12 :     ($rpts = shift @ARGV)
13 :     )
14 :     || die $usage;
15 :    
16 :     my @contigs = &gjoseqlib::read_fasta;
17 :     my $sz = 0;
18 :     foreach $_ (@contigs) { $sz += length($_->[2]) }
19 :     my $tmp_dir = '.';
20 :     my $tmp1 = "$tmp_dir/tmp1_$$.fasta";
21 : overbeek 1.2 open(TMP,"| sort | CSA_filter_unique > $tmp1 2> $rpts") || die "could not set up pipeline";
22 : overbeek 1.1 foreach $_ (@contigs)
23 :     {
24 :     my($id,undef,$seq) = @$_;
25 :     my $seqR = &SeedUtils::rev_comp($seq);
26 :     my $ln = length($seq);
27 :     my $i;
28 :     for ($i=0; ($i < (length($seq) - $N)); $i++)
29 :     {
30 :     &print_kmer(uc substr($seq,$i,$N),$id,'+',$i+1,\*TMP);
31 :     &print_kmer(uc substr($seqR,$i,$N),$id,'-',$ln-$i,\*TMP);
32 :     }
33 :     }
34 :     close(TMP);
35 :     open(KMERS,"<$tmp1") || die "could not open the kmers";
36 :     while (defined($_ = <KMERS>)) { print $_ }
37 :     close(KMERS);
38 :     unlink($tmp1);
39 :    
40 :     sub print_kmer {
41 :     my($seq,$id,$strand,$pos,$fh) = @_;
42 :    
43 :     my $i;
44 :     my $ln = length($seq);
45 :    
46 :     for ($i=2; ($i < $ln); $i += 3)
47 :     {
48 :     substr($seq,$i,1) = 'N';
49 :     }
50 :     print $fh join("\t",($seq,$id,$strand,$pos)),"\n";
51 :     }

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3