Mercurial > repos > pjbriggs > motif_tools

comparison Scan_IUPAC_output_each_match.pl @ 3:b4e22d2d7fa7 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add missing Perl scripts from previous uploads.
author pjbriggs
date Mon, 09 Apr 2018 04:45:06 -0400
parents
children
comparison
equal deleted inserted replaced
2:96a4ee27b08a 3:b4e22d2d7fa7
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 }