-
Notifications
You must be signed in to change notification settings - Fork 0
/
prot_motif_count.pl
executable file
·52 lines (40 loc) · 1.04 KB
/
prot_motif_count.pl
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
#!/usr/bin/perl
use strict;
use warnings;
# a script to find a defined motif in a protein sequence.
# wildcards in the motif should be represented by a .
# requires:
my $protein_fasta = $ARGV[0];
my $motif = uc $ARGV[1];
my $out_tab = $ARGV[2];
if (@ARGV != 3){
die "A script to get motif count in sequences from a protein fasta file\n\n requires 3 arguments:\n\n1) protein fasta file\n2) user-defined motif\n3) user-defined outfile\n\n";
}
#read fasta into hash
my %prot;
my $prot_id;
open PROT, $protein_fasta or die "no protein fasta file provided";
while (my $line = <PROT>){
chomp $line;
my $fc = substr($line, 0, 1);
if ($fc eq ">"){
$prot_id = substr($line,1);
}
else{
$prot{$prot_id} .= uc $line;
}
}
close PROT;
# print out seqs that match input motif
my $count;
close OUT;
open OUT, ">> $out_tab";
print OUT "Sequence header\tMotif count\n";
foreach $prot_id (keys %prot){
if ($prot{$prot_id} =~ /$motif/){
$count = () = $prot{$prot_id} =~ /$motif/gi;
close OUT;
open OUT, ">> $out_tab";
print OUT ">$prot_id\t$count\n";
}
}