scripts/remove_stop_codons.pl at master · reubwn/scripts · 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
#!/usr/bin/env perl

## author: reubwn May 2019

use strict;
use warnings;

use Bio::SeqIO;
use Getopt::Long;

my $usage = "
SYNOPSIS:
  Strips trailing asterisks (stop codons; '*') from protein sequences
  Also removes sequences with internal stop codons by default
  Can read from gzip

OPTIONS:
  -f|--fasta [FILE]   : fasta file [required]
  -s|--stop  [STRING] : set stop string (default '*')
  -i|--inplace        : do in-place replacement of file (with '*.bak' backup)
  -h|--help           : prints this help message
\n";

my ($fasta, $stop, $inplace, $help);

GetOptions (
  'f|fasta=s' => \$fasta,
  's|stop:s'  => \$stop,
  'i|inplace' => \$inplace,
  'h|help'    => \$help
);

die $usage if $help;
die $usage unless ($fasta);

my $search_string = quotemeta ($stop);

my (%stripped, %removed);

my $in;
if ($fasta =~ m/gz$/) { ## read from gzip
  open (my $zcat, "zcat $fasta |") or die $!;
  $in = Bio::SeqIO -> new ( -fh => $zcat, -format => "fasta" );
} else {
  $in = Bio::SeqIO -> new ( -file => $fasta, -format => "fasta" );
}

if ($inplace) { ## do in-place replacement of file
  print STDERR "[INFO] Doing in-place file replacement (backup saved to $fasta.bak)\n";
  ## save a backup
  `cp $fasta $fasta.bak`;

  open (my $TMP, ">$fasta.tmp") or die $!;
  while (my $seq_obj = $in->next_seq()) {
    my $seq_string = $seq_obj->seq();
    if ($seq_string =~ m/$search_string$/) { ## has terminal stop codon
      $seq_string =~ s/$search_string$//; ## remove it
      if ($seq_string =~ m/$search_string/) { ## also has internal stop codon
        $removed{$seq_obj->display_id()}++; ## dont print it
      } else {
        print $TMP ">" . $seq_obj->display_id() . "\n" . $seq_string . "\n";
        $stripped{$seq_obj->display_id()}++;
      }

    } elsif ($seq_string =~ m/$search_string/) { ## stop codon is internal; don't print
      $removed{$seq_obj->display_id()}++; ## dont print it
    } else {
      print $TMP ">" . $seq_obj->display_id() . "\n" . $seq_string . "\n";
    }
  }
  ## do the replacement
  `mv $fasta.tmp $fasta`;

} else { ## print to STDOUT
  while (my $seq_obj = $in->next_seq()) {
    my $seq_string = $seq_obj->seq();
    if ($seq_string =~ m/$search_string$/) { ## has terminal stop codon
      $seq_string =~ s/$search_string$//; ## remove it
      if ($seq_string =~ m/$search_string/) { ## also has internal stop codon
        $removed{$seq_obj->display_id()}++; ## dont print it
      } else {
        print STDOUT ">" . $seq_obj->display_id() . "\n" . $seq_string . "\n";
        $stripped{$seq_obj->display_id()}++;
      }

    } elsif ($seq_string =~ m/$search_string/) { ## stop codon is internal; don't print
      $removed{$seq_obj->display_id()}++; ## dont print it
    } else {
      print STDOUT ">" . $seq_obj->display_id() . "\n" . $seq_string . "\n";
    }
  }

}

print STDERR "[INFO] Number of seqs with trailing '$stop' stripped: " . (keys %stripped) . "\n";
print STDERR "[INFO] Number of seqs with internal '$stop' removed: " . (keys %removed) . "\n";
print STDERR "[INFO] Done " . `date`;