|
3
|
1 #! /usr/bin/perl
|
|
|
2
|
|
|
3 use strict;
|
|
|
4 use FileHandle;
|
|
|
5 use Bio::SeqIO;
|
|
|
6 #use Statistics::Descriptive;
|
|
|
7
|
|
|
8 #####
|
|
|
9 # Program to count all occurences of a particular REGEX
|
|
|
10 # in a file containing mutiple FASTA sequences.
|
|
|
11 # 11 September 2003. Ian Donaldson.
|
|
|
12 # Revised to convert IUPAC to regex
|
|
|
13 # Revised to read a multiple FASTA file
|
|
|
14 # was CountRegexGFF_IUPAC_1input_nosummary.pl
|
|
|
15 #####
|
|
|
16
|
|
|
17 #### File handles
|
|
|
18 my $input = new FileHandle;
|
|
|
19 my $output = new FileHandle;
|
|
|
20
|
|
|
21 #### Variables
|
|
|
22 my $file_number = 0;
|
|
|
23 my $count_fwd_regex = 0;
|
|
|
24 my $count_rvs_regex = 0;
|
|
|
25 my $count_all_regex = 0;
|
|
|
26 my $seq_tally = 0;
|
|
|
27 my @seq_totals = ();
|
|
|
28
|
|
|
29 #### Command line usage
|
|
|
30 if(@ARGV != 5) {
|
|
|
31 die ("USAGE:
|
|
|
32 $0
|
|
|
33 IUPAC
|
|
|
34 Multiple FASTA input file
|
|
|
35 Output
|
|
|
36 Label
|
|
|
37 Skip palindromic (0=F+R-default|1=F only)\n\n");
|
|
|
38 }
|
|
|
39
|
|
|
40 #### Search forward strand only?
|
|
|
41 my $skip = $ARGV[4];
|
|
|
42 unless($skip =~ /^[01]$/) {
|
|
|
43 die("Only accept 0 or 1 for Skip!\n");
|
|
|
44 }
|
|
|
45
|
|
|
46 #### Process IUPAC string
|
|
|
47 my $iupac = $ARGV[0];
|
|
|
48 chomp $iupac;
|
|
|
49 $iupac = uc($iupac);
|
|
|
50
|
|
|
51 if($iupac !~ /^[ACGTRYMKWSBDHVN]+$/) {
|
|
|
52 die("A non-IUPAC character was detected in your input string!\n");
|
|
|
53 }
|
|
|
54
|
|
|
55 #### Forward strand IUPAC
|
|
|
56 my @fwd_iupac_letters = split(//, $iupac);
|
|
|
57 my @fwd_regex_list = ();
|
|
|
58
|
|
|
59 foreach my $letter (@fwd_iupac_letters) {
|
|
|
60 my $converted_iupac = iupac2regex($letter);
|
|
|
61 push(@fwd_regex_list, $converted_iupac);
|
|
|
62 }
|
|
|
63
|
|
|
64 my $fwd_regex = join('', @fwd_regex_list);
|
|
|
65
|
|
|
66
|
|
|
67 #### Reverse strand IUPAC
|
|
|
68 my $revcomp_iupac = RevCompIUPAC($iupac);
|
|
|
69 my @rev_iupac_letters = split(//, $revcomp_iupac);
|
|
|
70 my @rev_regex_list = ();
|
|
|
71
|
|
|
72 foreach my $letter (@rev_iupac_letters) {
|
|
|
73 my $converted_iupac = iupac2regex($letter);
|
|
|
74 push(@rev_regex_list, $converted_iupac);
|
|
|
75 }
|
|
|
76
|
|
|
77 my $rvs_regex = join('', @rev_regex_list);
|
|
|
78
|
|
|
79 #### Other variables
|
|
|
80 my $label = $ARGV[3];
|
|
|
81
|
|
|
82 if($label !~ /^[\w\d]+$/) {
|
|
|
83 die("A non-letter/number character was detected in your label string!\n");
|
|
|
84 }
|
|
|
85
|
|
|
86 my $length = length($iupac);
|
|
|
87
|
|
|
88 #### Open output file
|
|
|
89 $output->open(">$ARGV[2]") or die "Could not open output file $ARGV[2]!\n";
|
|
|
90 #$output->print("##gff-version 2\n");
|
|
|
91
|
|
|
92 #if($skip == 0) {
|
|
|
93 # $output->print("##Pattern search: $iupac and $revcomp_iupac\n");
|
|
|
94 #}
|
|
|
95
|
|
|
96 #else {
|
|
|
97 # $output->print("##Pattern search: $iupac\n");
|
|
|
98 #}
|
|
|
99
|
|
|
100 #### Work thru FASTA entries in the input file with SeqIO
|
|
|
101 my $seqio = Bio::SeqIO->new(-file => "$ARGV[1]" , '-format' => 'Fasta');
|
|
|
102
|
|
|
103 while( my $seqobj = $seqio->next_seq() ) {
|
|
|
104 $seq_tally++;
|
|
|
105 my $this_seq_tally = 0;
|
|
|
106 my $sequence = $seqobj->seq(); # actual sequence as a string
|
|
|
107 my $seq_id = $seqobj->id(); # header
|
|
|
108 #print(">$seq_id\n$seq\n\n");
|
|
|
109
|
|
|
110 #$output->print(">$seq_id\n");
|
|
|
111
|
|
|
112 #### Clean up $sequence to leave only nucleotides
|
|
|
113 $sequence =~ s/[\s\W\d]//g;
|
|
|
114
|
|
|
115 while ($sequence =~ /($fwd_regex)/ig) {
|
|
|
116 $this_seq_tally++;
|
|
|
117 $count_fwd_regex++;
|
|
|
118 $count_all_regex++;
|
|
|
119
|
|
|
120 my $end_position = pos($sequence);
|
|
|
121 my $start_position = $end_position - ($length - 1);
|
|
|
122 $output->print("$seq_id\tRegexSearch\tCNS\t$start_position\t$end_position\t.\t+\t.\t$label\n");
|
|
|
123 }
|
|
|
124
|
|
|
125 #### Count reverse REGEX
|
|
|
126 unless($skip == 1) {
|
|
|
127 while ($sequence =~ /($rvs_regex)/ig) {
|
|
|
128 $this_seq_tally++;
|
|
|
129 $count_rvs_regex++;
|
|
|
130 $count_all_regex++;
|
|
|
131
|
|
|
132 my $end_position = pos($sequence);
|
|
|
133 my $start_position = $end_position - ($length - 1);
|
|
|
134 $output->print("$seq_id\tRegexSearch\tCNS\t$start_position\t$end_position\t.\t-\t.\t$label\n");
|
|
|
135 }
|
|
|
136
|
|
|
137 push(@seq_totals, $this_seq_tally);
|
|
|
138 #$output->print("$this_seq_tally matches\n");
|
|
|
139 }
|
|
|
140 }
|
|
|
141
|
|
|
142 #### Mean motifs per seq
|
|
|
143 #my $stat = Statistics::Descriptive::Full->new();
|
|
|
144 #$stat->add_data(@seq_totals);
|
|
|
145 #my $mean = $stat->mean();
|
|
|
146
|
|
|
147
|
|
|
148 #### Print a summary file
|
|
|
149 if($skip == 0) {
|
|
|
150 # $output->print("##Forward: $fwd_regex. Reverse: $rvs_regex.\n",
|
|
|
151 # "##$count_fwd_regex on the forward strand.\n",
|
|
|
152 # "##$count_rvs_regex on the reverse strand.\n",
|
|
|
153 # "##$count_all_regex in total.\n",
|
|
|
154 # "##$seq_tally sequences. Mean motifs per seq = $mean\n");
|
|
|
155 #
|
|
|
156 print STDOUT "There were $count_all_regex instances of $fwd_regex and $rvs_regex.\n";
|
|
|
157 }
|
|
|
158
|
|
|
159 if($skip == 1) {
|
|
|
160 # $output->print("##Forward: $fwd_regex.\n",
|
|
|
161 # "##$count_fwd_regex on the forward strand.\n",
|
|
|
162 # "##$seq_tally sequences. Mean motifs per seq = $mean\n");
|
|
|
163 #
|
|
|
164 print STDOUT "There were $count_fwd_regex instances of $fwd_regex on the forward strand.\n";
|
|
|
165 }
|
|
|
166
|
|
|
167 $output->close;
|
|
|
168
|
|
|
169 exit;
|
|
|
170
|
|
|
171 sub iupac2regex {
|
|
|
172 # Convert IUPAC codes to REGEX
|
|
|
173 my $iupac = shift;
|
|
|
174
|
|
|
175 #### Series of regexes to convert
|
|
|
176 if($iupac =~ /A/) { return 'A' }
|
|
|
177 if($iupac =~ /C/) { return 'C' }
|
|
|
178 if($iupac =~ /G/) { return 'G' }
|
|
|
179 if($iupac =~ /T/) { return 'T' }
|
|
|
180 if($iupac =~ /M/) { return '[AC]' }
|
|
|
181 if($iupac =~ /R/) { return '[AG]' }
|
|
|
182 if($iupac =~ /W/) { return '[AT]' }
|
|
|
183 if($iupac =~ /S/) { return '[CG]' }
|
|
|
184 if($iupac =~ /Y/) { return '[CT]' }
|
|
|
185 if($iupac =~ /K/) { return '[GT]' }
|
|
|
186 if($iupac =~ /V/) { return '[ACG]' }
|
|
|
187 if($iupac =~ /H/) { return '[ACT]' }
|
|
|
188 if($iupac =~ /D/) { return '[AGT]' }
|
|
|
189 if($iupac =~ /B/) { return '[CGT]' }
|
|
|
190 if($iupac =~ /N/) { return '[ACGT]' }
|
|
|
191
|
|
|
192 die("IUPAC not recognised by sub iupac2regex!\n");
|
|
|
193 }
|
|
|
194
|
|
|
195 sub RevCompIUPAC {
|
|
|
196 my $iupac_string = shift;
|
|
|
197 my @converted_list = ();
|
|
|
198
|
|
|
199 my @iupac_string_list = split(//, $iupac_string);
|
|
|
200
|
|
|
201 @iupac_string_list = reverse(@iupac_string_list);
|
|
|
202
|
|
|
203 foreach my $letter (@iupac_string_list) {
|
|
|
204 $letter =~ tr/ACGTRYMKWSBDHVN/TGCAYRKMWSVHDBN/;
|
|
|
205 push(@converted_list, $letter);
|
|
|
206 }
|
|
|
207
|
|
|
208 my $joined_list = join('', @converted_list);
|
|
|
209 return $joined_list;
|
|
|
210 }
|
