← Back to branch summary

~ubuntu-branches/ubuntu/trusty/bioperl/trusty

« back to all changes in this revision

Viewing changes to scripts/utilities/bp_nrdb.PLS

  • Committer: Package Import Robot
  • Author(s): Charles Plessy
  • Date: 2013-09-22 13:39:48 UTC
  • mfrom: (3.1.11 sid)
  • Revision ID: package-import@ubuntu.com-20130922133948-c6z62zegjyp7ztou
Tags: 1.6.922-1
http://bugs.debian.org/722910
* New upstream release.
* Replaces and Breaks grinder (<< 0.5.3-3~) because of overlaping contents.
  Closes: #722910
* Stop Replacing and Breaking bioperl ( << 1.6.9 ): not needed anymore. 

Show diffs side-by-side

added added

removed removed

Lines of Context:
scripts/utilities/bp_nrdb.PLS
1
 
#!perl -w
2
 
 
3
 
# Author Jason Stajich <jason-at-bioperl-dot-org>
4
 
# 
5
 
# Make a non-redundant database based on sequence (not on ID!)
6
 
# This script is still in progress but is intended to mimic what 
7
 
# Warren Gish's nrdb does
8
 
 
9
 
# It requires that Digest::MD5 is installed (for now)
10
 
 
11
 
=head1 NAME
12
 
 
13
 
bp_nrdb.PLS - a script to emulate Warren Gish's nrdb, make a unique sequence database from a set of input databases
14
 
 
15
 
=head1 SYNOPSIS
16
 
 
17
 
 
18
 
Usage: 
19
 
  bp_nrdb.PLS [options] file1 file2 file3
20
 
 
21
 
Alternative usage
22
 
  bp_nrdb.PLS -p [options] file1 id1 file2 id2 file3 id3
23
 
 
24
 
=head1 DESCRIPTION
25
 
 
26
 
This script will create a unique database of sequences
27
 
(quasi-nonredundant).  The options are:
28
 
 
29
 
   -o filename          - the filename the db is written (STDOUT by default)
30
 
   -a filename          - the filename to append the db to
31
 
   -l#                  - minimum required sequence length
32
 
   -i                   - do not check for duplicates
33
 
   -n#                  - max number of descriptions to report per seq
34
 
   -d#                  - delimiter to use between consecutive descriptions
35
 
   -p                   - use database id prefixes from command line
36
 
 
37
 
=head1 AUTHOR
38
 
 
39
 
Jason Stajich, jason-at-bioperl-dot-org
40
 
 
41
 
=cut
42
 
 
43
 
use strict;
44
 
use Bio::SeqIO;
45
 
use Getopt::Long;
46
 
 
47
 
use Digest::MD5 qw(md5_hex);
48
 
my ($output,$append,$min_len, 
49
 
    $no_duplicate_check,$desc_count,
50
 
    $delimiter, $expect_prefixes,$help);
51
 
$delimiter = ';';
52
 
 
53
 
GetOptions(
54
 
           'o|output:s'    => \$output,
55
 
           'a|append:s'    => \$append,
56
 
           'n:s'           => \$desc_count,
57
 
           'l:s'           => \$min_len,
58
 
           'd:s'           => \$delimiter,
59
 
           'p'             => \$expect_prefixes,
60
 
           'i'             => \$no_duplicate_check,
61
 
           'h'             => \$help,
62
 
           );
63
 
 
64
 
die("must supply a positive integer for -d") if ( defined $desc_count &&
65
 
                                                  ( $desc_count !~ /^\d+$/ ||
66
 
                                                    $desc_count < 1) );
67
 
die("must supply a positive integer for -l") if ( defined $min_len &&
68
 
                                                  ( $min_len !~ /^\d+$/ ||
69
 
                                                    $min_len < 1) );
70
 
my @files;
71
 
 
72
 
if( $help || ! @ARGV ) {
73
 
    exec('perldoc',$0);
74
 
    exit(0);
75
 
}
76
 
while( @ARGV ) {
77
 
    
78
 
    my ($file, $id) = (undef,'');
79
 
    if( $expect_prefixes ) {
80
 
        ($file,$id) = (shift @ARGV, shift @ARGV);
81
 
        if( ! $id ) { 
82
 
            die("Must provide 'name id' pairing of dbfile and id");
83
 
        }
84
 
    } else { 
85
 
        $file = shift @ARGV;
86
 
    }
87
 
    push @files, [ $file,$id];
88
 
}
89
 
 
90
 
 
91
 
my $out;
92
 
if( $append ) {
93
 
    $out = new Bio::SeqIO(-file => ">>$append");
94
 
} elsif( $output ) { 
95
 
    $out = new Bio::SeqIO(-file => ">$output");
96
 
} else {
97
 
    $out = new Bio::SeqIO(); # use STDOUT
98
 
}
99
 
 
100
 
my %unique;
101
 
my %seqcount;
102
 
my $counter = 0;
103
 
foreach my $pair ( @files ) {
104
 
    my ($file,$id) = @$pair;
105
 
    my $in = new Bio::SeqIO(-file => $file);
106
 
    while( my $seq = $in->next_seq ) {
107
 
        next if defined $min_len && $seq->length < $min_len;
108
 
        if( $id ) { 
109
 
            $seq->display_id("$id:".$seq->display_id);
110
 
        }
111
 
        my $s = lc($seq->seq());
112
 
        my $md5sum = md5_hex($s);
113
 
        if( $no_duplicate_check ) {
114
 
            $md5sum = $counter++;
115
 
        }
116
 
            
117
 
        if( defined $unique{$md5sum} ) {
118
 
            $seqcount{$md5sum}++;
119
 
            next if defined $desc_count && $seqcount{$md5sum++} > $desc_count;
120
 
            my $desc = $unique{$md5sum}->description;       
121
 
            my $id2 = sprintf("%s %s:%s %s",$delimiter,
122
 
                              $id,$seq->display_id,$seq->description);
123
 
            $unique{$md5sum}->desc($desc . $id2);
124
 
        } else { 
125
 
            $unique{$md5sum} = $seq;    
126
 
        }
127
 
    }
128
 
}
129
 
 
130
 
foreach my $seq ( values %unique ) {
131
 
    $out->write_seq($seq);
132
 
}
133
 
 
134
 
__END__