|
0
|
1 #!/usr/bin/perl -w
|
|
|
2
|
|
|
3 # a program to compute the frequency of each motif at each window in both upstream and downstream sequences flanking indels
|
|
|
4 # in a chromosome/genome.
|
|
|
5 # the first input is a TABULAR format file containing the motif names and sequences, such that the file consists of two
|
|
|
6 # columns: the left column represents the motif names and the right column represents the motif sequence, one line per motif.
|
|
|
7 # the second input is a TABULAR format file containing the upstream and downstream sequences flanking indels, one line per indel.
|
|
|
8 # the fourth input is an integer number representing the window size according to which the upstream and downstream sequences
|
|
|
9 # flanking each indel will be divided.
|
|
|
10 # the first output is a TABULAR format file containing the windows into which both upstream and downstream sequences flanking
|
|
|
11 # indels are divided.
|
|
|
12 # the second output is a TABULAR format file containing the motifs and their corresponding frequencies at each window in both
|
|
|
13 # upstream and downstream sequences flanking indels, one line per motif.
|
|
|
14
|
|
|
15 use strict;
|
|
|
16 use warnings;
|
|
|
17
|
|
|
18 #variable to handle the falnking sequences information
|
|
|
19 my $sequence = "";
|
|
|
20 my $upstreamFlankingSequence = "";
|
|
|
21 my $downstreamFlankingSequence = "";
|
|
|
22 my $discardedSequenceLength = 0;
|
|
|
23 my $lengthOfDownstreamFlankingSequenceAfterTrimming = 0;
|
|
|
24
|
|
|
25 #variable to handle the window information
|
|
|
26 my $window = "";
|
|
|
27 my $windowStartIndex = 0;
|
|
|
28 my $windowNumber = 0;
|
|
|
29 my $totalWindowsNumber = 0;
|
|
|
30 my $totalNumberOfWindowsInUpstreamSequence = 0;
|
|
|
31 my $totalNumberOfWindowsInDownstreamSequence = 0;
|
|
|
32 my $totalWindowsNumberInBothFlankingSequences = 0;
|
|
|
33 my $totalWindowsNumberInMotifCountersTwoDimArray = 0;
|
|
|
34 my $upstreamAndDownstreamFlankingSequencesWindows = "";
|
|
|
35
|
|
|
36 #variable to handle the motif information
|
|
|
37 my $motif = "";
|
|
|
38 my $motifSequence = "";
|
|
|
39 my $motifNumber = 0;
|
|
|
40 my $totalMotifsNumber = 0;
|
|
|
41
|
|
|
42 #arrays to sotre window and motif data
|
|
|
43 my @windowsArray = ();
|
|
|
44 my @motifNamesArray = ();
|
|
|
45 my @motifSequencesArray = ();
|
|
|
46 my @motifCountersTwoDimArray = ();
|
|
|
47
|
|
|
48 #variables to store line counter values
|
|
|
49 my $lineCounter1 = 0;
|
|
|
50 my $lineCounter2 = 0;
|
|
|
51
|
|
|
52 # check to make sure having correct files
|
|
|
53 my $usage = "usage: compute_motifs_frequency.pl [TABULAR.in] [TABULAR.in] [windowSize] [TABULAR.out] [TABULAR.out]\n";
|
|
|
54 die $usage unless @ARGV == 5;
|
|
|
55
|
|
|
56 #get the input and output arguments
|
|
|
57 my $motifsInputFile = $ARGV[0];
|
|
|
58 my $indelFlankingSequencesInputFile = $ARGV[1];
|
|
|
59 my $windowSize = $ARGV[2];
|
|
|
60 my $indelFlankingSequencesWindowsOutputFile = $ARGV[3];
|
|
|
61 my $motifFrequenciesOutputFile = $ARGV[4];
|
|
|
62
|
|
|
63 #open the input and output files
|
|
|
64 open (INPUT1, "<", $motifsInputFile) || die("Could not open file $motifsInputFile \n");
|
|
|
65 open (INPUT2, "<", $indelFlankingSequencesInputFile) || die("Could not open file $indelFlankingSequencesInputFile \n");
|
|
|
66 open (OUTPUT1, ">", $indelFlankingSequencesWindowsOutputFile) || die("Could not open file $indelFlankingSequencesWindowsOutputFile \n");
|
|
|
67 open (OUTPUT2, ">", $motifFrequenciesOutputFile) || die("Could not open file $motifFrequenciesOutputFile \n");
|
|
|
68
|
|
|
69 #store the motifs input file in the array @motifsData
|
|
|
70 my @motifsData = <INPUT1>;
|
|
|
71
|
|
|
72 #iterated through the motifs (lines) of the motifs input file
|
|
|
73 foreach $motif (@motifsData){
|
|
|
74 chomp ($motif);
|
|
|
75 #print ($motif . "\n");
|
|
|
76
|
|
|
77 #split the motif data into its name and its sequence
|
|
|
78 my @motifNameAndSequenceArray = split(/\t/, $motif);
|
|
|
79
|
|
|
80 #store the name of the motif into the array @motifNamesArray
|
|
|
81 push @motifNamesArray, $motifNameAndSequenceArray[0];
|
|
|
82
|
|
|
83 #store the sequence of the motif into the array @motifSequencesArray
|
|
|
84 push @motifSequencesArray, $motifNameAndSequenceArray[1];
|
|
|
85 }
|
|
|
86
|
|
|
87 #compute the size of the motif names array
|
|
|
88 $totalMotifsNumber = @motifNamesArray;
|
|
|
89
|
|
|
90 #store the input file in the array @sequencesData
|
|
|
91 my @sequencesData = <INPUT2>;
|
|
|
92
|
|
|
93 #iterated through the sequences of the second input file in order to create windwos file
|
|
|
94 foreach $sequence (@sequencesData){
|
|
|
95 chomp ($sequence);
|
|
|
96 $lineCounter1++;
|
|
|
97
|
|
|
98 my @indelAndSequenceArray = split(/\t/, $sequence);
|
|
|
99
|
|
|
100 #get the upstream falnking sequence
|
|
|
101 $upstreamFlankingSequence = $indelAndSequenceArray[3];
|
|
|
102
|
|
|
103 #if the window size is 0, then the whole upstream will be one window only
|
|
|
104 if ($windowSize == 0){
|
|
|
105 $totalNumberOfWindowsInUpstreamSequence = 1;
|
|
|
106 $windowSize = length ($upstreamFlankingSequence);
|
|
|
107 }
|
|
|
108 else{
|
|
|
109 #compute the total number of windows into which the upstream flanking sequence will be divided
|
|
|
110 $totalNumberOfWindowsInUpstreamSequence = length ($upstreamFlankingSequence) / $windowSize;
|
|
|
111
|
|
|
112 #compute the length of the subsequence to be discared from the upstream flanking sequence if any
|
|
|
113 $discardedSequenceLength = length ($upstreamFlankingSequence) % $windowSize;
|
|
|
114
|
|
|
115 #check if the sequence could be split into windows of equal sizes
|
|
|
116 if ($discardedSequenceLength != 0) {
|
|
|
117 #trim the upstream flanking sequence
|
|
|
118 $upstreamFlankingSequence = substr($upstreamFlankingSequence, $discardedSequenceLength);
|
|
|
119 }
|
|
|
120 }
|
|
|
121
|
|
|
122 #split the upstream flanking sequence into windows
|
|
|
123 for ($windowNumber = 0; $windowNumber < $totalNumberOfWindowsInUpstreamSequence; $windowNumber++){
|
|
|
124 $windowStartIndex = $windowNumber * $windowSize;
|
|
|
125 print OUTPUT1 (substr($upstreamFlankingSequence, $windowStartIndex, $windowSize) . "\t");
|
|
|
126 }
|
|
|
127
|
|
|
128 #add a column representing the indel
|
|
|
129 print OUTPUT1 ("indel" . "\t");
|
|
|
130
|
|
|
131 #get the downstream falnking sequence
|
|
|
132 $downstreamFlankingSequence = $indelAndSequenceArray[4];
|
|
|
133
|
|
|
134 #if the window size is 0, then the whole upstream will be one window only
|
|
|
135 if ($windowSize == 0){
|
|
|
136 $totalNumberOfWindowsInDownstreamSequence = 1;
|
|
|
137 $windowSize = length ($downstreamFlankingSequence);
|
|
|
138 }
|
|
|
139 else{
|
|
|
140 #compute the total number of windows into which the downstream flanking sequence will be divided
|
|
|
141 $totalNumberOfWindowsInDownstreamSequence = length ($downstreamFlankingSequence) / $windowSize;
|
|
|
142
|
|
|
143 #compute the length of the subsequence to be discared from the upstream flanking sequence if any
|
|
|
144 $discardedSequenceLength = length ($downstreamFlankingSequence) % $windowSize;
|
|
|
145
|
|
|
146 #check if the sequence could be split into windows of equal sizes
|
|
|
147 if ($discardedSequenceLength != 0) {
|
|
|
148 #compute the length of the sequence to be discarded
|
|
|
149 $lengthOfDownstreamFlankingSequenceAfterTrimming = length ($downstreamFlankingSequence) - $discardedSequenceLength;
|
|
|
150
|
|
|
151 #trim the downstream flanking sequence
|
|
|
152 $downstreamFlankingSequence = substr($downstreamFlankingSequence, 0, $lengthOfDownstreamFlankingSequenceAfterTrimming);
|
|
|
153 }
|
|
|
154 }
|
|
|
155
|
|
|
156 #split the downstream flanking sequence into windows
|
|
|
157 for ($windowNumber = 0; $windowNumber < $totalNumberOfWindowsInDownstreamSequence; $windowNumber++){
|
|
|
158 $windowStartIndex = $windowNumber * $windowSize;
|
|
|
159 print OUTPUT1 (substr($downstreamFlankingSequence, $windowStartIndex, $windowSize) . "\t");
|
|
|
160 }
|
|
|
161
|
|
|
162 print OUTPUT1 ("\n");
|
|
|
163 }
|
|
|
164
|
|
|
165 #compute the total number of windows on both upstream and downstream sequences flanking the indel
|
|
|
166 $totalWindowsNumberInBothFlankingSequences = $totalNumberOfWindowsInUpstreamSequence + $totalNumberOfWindowsInDownstreamSequence;
|
|
|
167
|
|
|
168 #add an additional cell to store the name of the motif and another one for the indel itself
|
|
|
169 $totalWindowsNumberInMotifCountersTwoDimArray = $totalWindowsNumberInBothFlankingSequences + 1 + 1;
|
|
|
170
|
|
|
171 #initialize the two dimensional array $motifCountersTwoDimArray. the first column will be initialized with motif names
|
|
|
172 for ($motifNumber = 0; $motifNumber < $totalMotifsNumber; $motifNumber++){
|
|
|
173
|
|
|
174 for ($windowNumber = 0; $windowNumber < $totalWindowsNumberInMotifCountersTwoDimArray; $windowNumber++){
|
|
|
175
|
|
|
176 if ($windowNumber == 0){
|
|
|
177 $motifCountersTwoDimArray [$motifNumber] [0] = $motifNamesArray[$motifNumber];
|
|
|
178 }
|
|
|
179 elsif ($windowNumber == $totalNumberOfWindowsInUpstreamSequence + 1){
|
|
|
180 $motifCountersTwoDimArray [$motifNumber] [$windowNumber] = "indel";
|
|
|
181 }
|
|
|
182 else{
|
|
|
183 $motifCountersTwoDimArray [$motifNumber] [$windowNumber] = 0;
|
|
|
184 }
|
|
|
185 }
|
|
|
186 }
|
|
|
187
|
|
|
188 close(OUTPUT1);
|
|
|
189
|
|
|
190 #open the file the contains the windows of the upstream and downstream flanking sequences, which is actually the first output file
|
|
|
191 open (INPUT3, "<", $indelFlankingSequencesWindowsOutputFile) || die("Could not open file $indelFlankingSequencesWindowsOutputFile \n");
|
|
|
192
|
|
|
193 #store the first output file containing the windows of both upstream and downstream flanking sequences in the array @windowsData
|
|
|
194 my @windowsData = <INPUT3>;
|
|
|
195
|
|
|
196 #iterated through the lines of the first output file. Each line represents
|
|
|
197 #the windows of the upstream and downstream flanking sequences of an indel
|
|
|
198 foreach $upstreamAndDownstreamFlankingSequencesWindows (@windowsData){
|
|
|
199
|
|
|
200 chomp ($upstreamAndDownstreamFlankingSequencesWindows);
|
|
|
201 $lineCounter2++;
|
|
|
202
|
|
|
203 #split both upstream and downstream flanking sequences into their windows
|
|
|
204 my @windowsArray = split(/\t/, $upstreamAndDownstreamFlankingSequencesWindows);
|
|
|
205
|
|
|
206 $totalWindowsNumber = @windowsArray;
|
|
|
207
|
|
|
208 #iterate through the windows to search for matched motifs and increment their corresponding counters accordingly
|
|
|
209 WINDOWS:
|
|
|
210 for ($windowNumber = 0; $windowNumber < $totalWindowsNumber; $windowNumber++){
|
|
|
211
|
|
|
212 #get the window
|
|
|
213 $window = $windowsArray[$windowNumber];
|
|
|
214
|
|
|
215 #if the window is the one that contains the indel, then skip the indel window
|
|
|
216 if ($window eq "indel") {
|
|
|
217 next WINDOWS;
|
|
|
218 }
|
|
|
219 else{ #iterated through the motif sequences to check their occurrences in the winodw
|
|
|
220 #and increment their corresponding counters accordingly
|
|
|
221
|
|
|
222 for ($motifNumber = 0; $motifNumber < $totalMotifsNumber; $motifNumber++){
|
|
|
223 #get the motif sequence
|
|
|
224 $motifSequence = $motifSequencesArray[$motifNumber];
|
|
|
225
|
|
|
226 #if the motif is found in the window, then increment its corresponding counter
|
|
|
227 if ($window =~ m/$motifSequence/i){
|
|
|
228 $motifCountersTwoDimArray [$motifNumber] [$windowNumber + 1]++;
|
|
|
229 }
|
|
|
230 }
|
|
|
231 }
|
|
|
232 }
|
|
|
233 }
|
|
|
234
|
|
|
235 #store the motif counters values in the second output file
|
|
|
236 for ($motifNumber = 0; $motifNumber < $totalMotifsNumber; $motifNumber++){
|
|
|
237
|
|
|
238 for ($windowNumber = 0; $windowNumber <= $totalWindowsNumber; $windowNumber++){
|
|
|
239
|
|
|
240 print OUTPUT2 $motifCountersTwoDimArray [$motifNumber] [$windowNumber] . "\t";
|
|
|
241 #print ($motifCountersTwoDimArray [$motifNumber] [$windowNumber] . " ");
|
|
|
242 }
|
|
|
243 print OUTPUT2 "\n";
|
|
|
244 #print ("\n");
|
|
|
245 }
|
|
|
246
|
|
|
247 #close the input and output files
|
|
|
248 close(OUTPUT2);
|
|
|
249 close(OUTPUT1);
|
|
|
250 close(INPUT3);
|
|
|
251 close(INPUT2);
|
|
|
252 close(INPUT1); |
