[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.3 - (view) (download) (as text)

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

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3