prokka/bin/prokka-uniprot_to_fasta_db at 9e947876121a00a8bf2eb5e54120d0389e22a40c · tseemann/prokka · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
#!/usr/bin/perl -w
use strict;
use SWISS::Entry;
use SWISS::KW;
#use SWISS::OC;
use Data::Dumper;

my(@Options, $verbose, $frag, $evlev, $minlen, $sep, $blank, $term, $hypo);
setOptions();

my $HYPO = 'hypothetical protein';

# Read an entire record at a time
local $/ = "\n//\n";

my $in=0;
my $out=0;

while (<ARGV>)
{
  # Read the entry
  my $entry = SWISS::Entry->fromText($_);
  $in++;

  # Immediately reject partial genes
  next if not $frag and $entry->isFragment;
  next if not $frag and $entry->DEs->hasFragment;

  # Too short?
  next if $minlen > 0 and length($entry->SQs->seq) < $minlen;

  # Reject on evidence level
  # grep ^PE uniprot_sprot.dat | sort | uniq -c
  #  74284 PE   1: Evidence at protein level;
  #  67762 PE   2: Evidence at transcript level;
  # 376894 PE   3: Inferred from homology;
  #  14424 PE   4: Predicted;
  #   1884 PE   5: Uncertain;
  if ($evlev < 5) {
    $entry->PE->text =~ m/^(\d+)/;
    next unless $1 <= $evlev;
  }

  # Only specified organism class
  if ($term) {
    my $tax = $entry->OCs->list or next;
    next unless grep { $_ eq $term } @$tax ;
  }

  # /gene code
  my $gene = $entry->GNs->getFirst || '';
  $gene = '' if $gene =~ m/\d{2}/ or $gene =~ m/\./;

  my $ec = '';
#  my $prod = 'hypothetical protein';
  my $prod = '';

    for my $de ($entry->DEs->elements) {
      if ($de->type eq 'EC') {
    $ec .= $de->text;
#   last;
      }
      elsif ($de->type eq 'Full' and $de->category eq 'RecName') {
    $prod = $de->text;
    if ($prod =~ m/^UPF\d|^Uncharacterized protein|^ORF|^Protein /) {
      $prod = $HYPO;
    }
      }
    }
      $ec =~ s/^\D*//;
      $ec =~ s/EC /;/g;

  $prod ||= $HYPO;

  # skip hypthetical proteins, unless user has overridden this with --hypo
  next if !$hypo and $prod eq $HYPO;

  $ec ||= $blank;
  $gene ||= $blank;
  $prod ||= $blank;

  print STDERR join("\t", $entry->AC, $ec, $gene, $prod), "\n" if $verbose;
  print ">", $entry->AC, " $ec$sep$gene$sep$prod\n", $entry->SQs->seq, "\n";

  $out++;

}

#print STDERR "\n";

#----------------------------------------------------------------------
# Option setting routines

sub setOptions {
  use Getopt::Long;

  @Options = (
    {OPT=>"help",    VAR=>\&usage,             DESC=>"This help"},
    {OPT=>"verbose!", VAR=>\$verbose, DEFAULT=>0, DESC=>"Verbose output"},
    {OPT=>"separator=s",   VAR=>\$sep, DEFAULT=>'~~~', DESC=>"Separator for gene/EC/product"},
    {OPT=>"blank=s",   VAR=>\$blank, DEFAULT=>'', DESC=>"Replace empty gene/EC/product with this"},
    {OPT=>"evidence=i",   VAR=>\$evlev, DEFAULT=>2, DESC=>"1=prot 2=mrna 3=homol 4=pred 5=unsure"},
    {OPT=>"fragments!",   VAR=>\$frag, DEFAULT=>0, DESC=>"Include 'DE Flags: Fragment;' entries"},
    {OPT=>"minlen=i",   VAR=>\$minlen, DEFAULT=>20, DESC=>"Minimum peptide length"},
    {OPT=>"term=s",   VAR=>\$term, DEFAULT=>'', DESC=>"Lineage must contain this term eg. 'Bacteria'"},
    {OPT=>"hypo!",   VAR=>\$hypo, DEFAULT=>0, DESC=>"Don't filter out hypothetical proteins"},
  );

  #(!@ARGV) && (usage());

  &GetOptions(map {$_->{OPT}, $_->{VAR}} @Options) || usage();

  # Now setup default values.
  foreach (@Options) {
    if (defined($_->{DEFAULT}) && !defined(${$_->{VAR}})) {
      ${$_->{VAR}} = $_->{DEFAULT};
    }
  }
}

sub usage {
  print "Usage: $0 [options] <uniprot.dat> > sprot\n";
  foreach (@Options) {
    printf "  --%-13s %s%s.\n",$_->{OPT},$_->{DESC},
           defined($_->{DEFAULT}) ? " (default '$_->{DEFAULT}')" : "";
  }
  exit(1);
}

#----------------------------------------------------------------------