[Bio] / FigMetagenomeTools / exact_duplicates.pl Repository:
ViewVC logotype

Annotation of /FigMetagenomeTools/exact_duplicates.pl

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : olson 1.1 #!/usr/bin/perl -w
2 :    
3 :     # count the exact duplicate reads
4 :    
5 :     use strict;
6 :     use Rob;
7 :    
8 :    
9 :     my $file=shift || die "$0 <file> <full output>\nIf full output is not set (boolean) only a summary is printed\n";
10 :     my $fo=shift;
11 :    
12 :     if ($file =~ /gz$/)
13 :     {
14 :     open(IN, "gunzip -c $file|") || die "Can't open pipe";
15 :     }
16 :     else {open (IN, $file) || die "$file?"}
17 :    
18 :    
19 :    
20 :    
21 :     my $tag; my $seq; my $count;
22 :     while (<IN>) {
23 :     chomp;
24 :     if (/^>/) {
25 :     if ($seq) {
26 :     $seq=uc($seq);
27 :     if ($count->{$seq}) {
28 :     push @{$count->{$seq}}, $tag;
29 :     }
30 :     else {
31 :     # if we don't have it in the forward phase, check the reverse
32 :     my $rc=Rob->rc($seq);
33 :     if ($count->{$rc}) {push @{$count->{$rc}}, $tag}
34 :     else {push @{$count->{$seq}}, $tag;}
35 :     }
36 :     undef $seq;
37 :     }
38 :     $tag = $_;
39 :     $tag =~ s/^>//;
40 :     } else {
41 :     $seq .= $_;
42 :     }
43 :     }
44 :    
45 :     if ($fo)
46 :     {
47 :     foreach my $seq (sort {scalar @{$count->{$b}} <=> scalar @{$count->{$a}}} keys %$count)
48 :     {
49 :     if (scalar @{$count->{$seq}} > 1)
50 :     {
51 :     print join "\t", scalar @{$count->{$seq}}, @{$count->{$seq}}, "\n";
52 :     }
53 :     }
54 :     }
55 :     else
56 :     {
57 :     my %allcounts;
58 :     foreach my $seq (sort {scalar @{$count->{$b}} <=> scalar @{$count->{$a}}} keys %$count)
59 :     {
60 :     if (scalar @{$count->{$seq}} > 1) {$allcounts{scalar @{$count->{$seq}}}++}
61 :     }
62 :     print "No. Seqs.\tNo dups\n";
63 :     my $total;
64 :     map {print "$_\t$allcounts{$_}\n"; $total+=($_-1) * $allcounts{$_}} sort {$allcounts{$b} <=> $allcounts{$a}} keys %allcounts;
65 :     print "TOTAL\t$total\n";
66 :     }
67 :    

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3