質問編集履歴
3
systemが指定している最初のファイルのコードを追加
test
CHANGED
|
File without changes
|
test
CHANGED
|
@@ -42,29 +42,375 @@
|
|
|
42
42
|
|
|
43
43
|
system "perl ./Scripts/1_Preprocessing.pl";
|
|
44
44
|
|
|
45
|
-
|
|
46
|
-
|
|
47
|
-
|
|
48
|
-
|
|
49
|
-
|
|
50
|
-
|
|
51
|
-
s
|
|
52
|
-
|
|
53
|
-
s
|
|
54
|
-
|
|
55
|
-
s
|
|
56
|
-
|
|
57
|
-
|
|
58
|
-
|
|
59
|
-
|
|
60
|
-
|
|
61
|
-
|
|
62
|
-
|
|
63
|
-
s
|
|
64
|
-
|
|
65
|
-
|
|
66
|
-
|
|
67
|
-
|
|
45
|
+
以下略
|
|
46
|
+
|
|
47
|
+
|
|
48
|
+
|
|
49
|
+
以下1_Preprocessing.plのコード
|
|
50
|
+
|
|
51
|
+
#!/usr/bin/perl
|
|
52
|
+
|
|
53
|
+
use strict;
|
|
54
|
+
|
|
55
|
+
use warnings;
|
|
56
|
+
|
|
57
|
+
|
|
58
|
+
|
|
59
|
+
#2018/12/04 1_2 Strip_primers change to use the Primer_Cleaner.pl
|
|
60
|
+
|
|
61
|
+
#2018/07/01 1_2_Strip_primers.pl were changed to permit no primer seq
|
|
62
|
+
|
|
63
|
+
#2018/01/29 change Blast to Usearch_global
|
|
64
|
+
|
|
65
|
+
#2018/01/20 output a log file at 1_1 step
|
|
66
|
+
|
|
67
|
+
#2018/01/05 add compress option
|
|
68
|
+
|
|
69
|
+
|
|
70
|
+
|
|
71
|
+
#This script execute below scripts in each sample
|
|
72
|
+
|
|
73
|
+
#1_1 Merge_paird_reads.pl
|
|
74
|
+
|
|
75
|
+
#1_2_Strip_primers.pl
|
|
76
|
+
|
|
77
|
+
#1_3_Quality_filter.pl
|
|
78
|
+
|
|
79
|
+
|
|
80
|
+
|
|
81
|
+
|
|
82
|
+
|
|
83
|
+
#Setting
|
|
84
|
+
|
|
85
|
+
my ($db, $primer, $diff, $separate, $depth, $size, $identity, $identity2, $dictionary, $family, $temporary, $compress, $trim);
|
|
86
|
+
|
|
87
|
+
open (SET, "<", "Setting.txt") or die("error:$!");
|
|
88
|
+
|
|
89
|
+
while(<SET>){
|
|
90
|
+
|
|
91
|
+
if($_ =~ /^DB\s*=\s*(\S+)/){$db = $1;}
|
|
92
|
+
|
|
93
|
+
elsif($_ =~ /^Primers\s*=\s*(\S+)/){my $temp = $1;if($temp =~ /^no$/i){$primer = "No";}else{$primer = $temp;}}
|
|
94
|
+
|
|
95
|
+
elsif($_ =~ /^MaxDiff\s*=\s*(\S+)/){$diff = $1;}
|
|
96
|
+
|
|
97
|
+
elsif($_ =~ /^Divide\s*=\s*(\S+)/){$separate = $1;}
|
|
98
|
+
|
|
99
|
+
elsif($_ =~ /^Depth\s*=\s*(\S+)/){$depth = $1;}
|
|
100
|
+
|
|
101
|
+
elsif($_ =~ /^Length\s*=\s*(\d+)/){$size = $1;}
|
|
102
|
+
|
|
103
|
+
elsif($_ =~ /^UIdentity\s*=\s*(\S+)/){$identity = $1;}
|
|
104
|
+
|
|
105
|
+
elsif($_ =~ /^LIdentity\s*=\s*(\S+)/){$identity2 = $1;}
|
|
106
|
+
|
|
107
|
+
elsif($_ =~ /^Dictionary\s*=\s*(\S+)/){$dictionary = $1;}
|
|
108
|
+
|
|
109
|
+
elsif($_ =~ /^Family\s*=\s*(\S+)/){$family = $1;}
|
|
110
|
+
|
|
111
|
+
elsif($_ =~ /^Temporary\s*=\s*(\S+)/){$temporary = $1;}
|
|
112
|
+
|
|
113
|
+
elsif($_ =~ /^Compress\s*=\s*(\S+)/){$compress = $1;}
|
|
114
|
+
|
|
115
|
+
}
|
|
116
|
+
|
|
117
|
+
close(SET);
|
|
118
|
+
|
|
119
|
+
if($diff =~ /^length$/i){$trim = 1;}
|
|
120
|
+
|
|
121
|
+
unless($diff =~ /\d+/){$diff = 2;}
|
|
122
|
+
|
|
123
|
+
unless($separate){$separate = 0;}
|
|
124
|
+
|
|
125
|
+
if($separate =~ /^yes$/i){$separate = 1;}
|
|
126
|
+
|
|
127
|
+
else{$separate = 0;}
|
|
128
|
+
|
|
129
|
+
unless($db){print "Error!! Check the DB name in Setting.txt.\n"; exit;}
|
|
130
|
+
|
|
131
|
+
unless($primer){print "Error!! Check the Primer file name in Setting.txt.\n"; exit;}
|
|
132
|
+
|
|
133
|
+
unless($depth){print "Error!! Check the number of Depth in Setting.txt.\n"; exit;}
|
|
134
|
+
|
|
135
|
+
unless($size){print "Error!! Check the number of Length in Setting.txt.\n"; exit;}
|
|
136
|
+
|
|
137
|
+
unless($identity){print "Error!! Check the upper cut line of Homology Search in Setting.txt.\n"; exit;}
|
|
138
|
+
|
|
139
|
+
unless($identity2){print "Error!! Check the lower cut line of Homology Search in Setting.txt.\n"; exit;}
|
|
140
|
+
|
|
141
|
+
unless($temporary){$temporary = "YES";}
|
|
142
|
+
|
|
143
|
+
unless($compress){$compress = "YES";}
|
|
144
|
+
|
|
145
|
+
unless(-f "./DataBase/$db"){print "No DataBase file or DataBase name didn't match in Setting.txt\n"; exit;}
|
|
146
|
+
|
|
147
|
+
unless($primer =~ /^no$/i){
|
|
148
|
+
|
|
149
|
+
unless(-f "./DataBase/$primer"){print "No Primers file or file name didn't match in Setting.txt\n"; exit;}
|
|
150
|
+
|
|
151
|
+
}
|
|
152
|
+
|
|
153
|
+
unless(-f "./Dictionary/$dictionary"){
|
|
154
|
+
|
|
155
|
+
if($dictionary =~ /^no$/i){undef($dictionary);}
|
|
156
|
+
|
|
157
|
+
else{print "The file of Dictionary for common name isn't or don't match the name in Setting.txt\n"; exit;}
|
|
158
|
+
|
|
159
|
+
}
|
|
160
|
+
|
|
161
|
+
unless(-f "./Dictionary/$family"){
|
|
162
|
+
|
|
163
|
+
if($family =~ /^no$/i){undef($family);}
|
|
164
|
+
|
|
165
|
+
else{print "The file of Dictionary for Family name isn't or don't match the name in Setting.txt\n"; exit;}
|
|
166
|
+
|
|
167
|
+
}
|
|
168
|
+
|
|
169
|
+
|
|
170
|
+
|
|
171
|
+
#Get primer sequences
|
|
172
|
+
|
|
173
|
+
my $countp = 0;
|
|
174
|
+
|
|
175
|
+
my (%forward, %reverse, $namep);
|
|
176
|
+
|
|
177
|
+
my ($primerF, $primerR);
|
|
178
|
+
|
|
179
|
+
|
|
180
|
+
|
|
181
|
+
unless($primer =~ /^no$/i){
|
|
182
|
+
|
|
183
|
+
open (DATA2, "<", "./DataBase/$primer") or die("error:$!");
|
|
184
|
+
|
|
185
|
+
while(<DATA2>){
|
|
186
|
+
|
|
187
|
+
if($_ =~ /\#/){next;}
|
|
188
|
+
|
|
189
|
+
chomp($_);
|
|
190
|
+
|
|
191
|
+
if($_ =~ /Forward/){$countp = 1;next;}
|
|
192
|
+
|
|
193
|
+
if($_ =~ /Reverse/){$countp = 2;next;}
|
|
194
|
+
|
|
195
|
+
if($countp == 1){
|
|
196
|
+
|
|
197
|
+
if($_ =~ /^>(.+)/){$namep = $1;}
|
|
198
|
+
|
|
199
|
+
elsif($_ =~ /^[a-z]/i){$_ =~ tr/[a-z]/[A-Z]/; $forward{$namep} = $_;}
|
|
200
|
+
|
|
201
|
+
}
|
|
202
|
+
|
|
203
|
+
if($countp == 2){
|
|
204
|
+
|
|
205
|
+
if($_ =~ /^>(.+)/){$namep = $1;}
|
|
206
|
+
|
|
207
|
+
elsif($_ =~ /^[a-z]/i){
|
|
208
|
+
|
|
209
|
+
my $rev = reverse($_);
|
|
210
|
+
|
|
211
|
+
$rev =~ tr/ATGCURYMKDHBV/TACGAYRKMHDVB/;
|
|
212
|
+
|
|
213
|
+
$reverse{$namep} = $rev;}
|
|
214
|
+
|
|
215
|
+
}
|
|
216
|
+
|
|
217
|
+
}
|
|
218
|
+
|
|
219
|
+
close(DATA2);
|
|
220
|
+
|
|
221
|
+
unless(%forward){print "Error!! No primer seqeunce in $primer\n"; exit;}
|
|
222
|
+
|
|
223
|
+
unless(%reverse){print "Error!! No primer seqeunce in $primer\n"; exit;}
|
|
224
|
+
|
|
225
|
+
foreach(keys %forward){
|
|
226
|
+
|
|
227
|
+
unless($reverse{$_}){print "Error!! Paird primer need to be the same name!\nCheck $primer\n"; exit;}
|
|
228
|
+
|
|
229
|
+
}
|
|
230
|
+
|
|
231
|
+
foreach(keys %reverse){
|
|
232
|
+
|
|
233
|
+
unless($forward{$_}){print "Error!! Paird primer need to be the same name!\nCheck $primer\n"; exit;}
|
|
234
|
+
|
|
235
|
+
}
|
|
236
|
+
|
|
237
|
+
|
|
238
|
+
|
|
239
|
+
my $primer_num = keys %forward;
|
|
240
|
+
|
|
241
|
+
if($primer_num == 1){$separate = 0;}
|
|
242
|
+
|
|
243
|
+
if($trim and $primer_num == 1){
|
|
244
|
+
|
|
245
|
+
foreach(keys %forward){$primerF = length($forward{$_}); $primerR = length($reverse{$_});}
|
|
246
|
+
|
|
247
|
+
}elsif($trim and $primer_num > 1){
|
|
248
|
+
|
|
249
|
+
print "Error! \"MaxDiff = length\" in Setting.txt can use when primer pair is one pair\n"; exit;
|
|
250
|
+
|
|
251
|
+
}
|
|
252
|
+
|
|
253
|
+
if($trim){$separate = 0;}
|
|
254
|
+
|
|
255
|
+
}
|
|
256
|
+
|
|
257
|
+
|
|
258
|
+
|
|
259
|
+
#Database and Usearch check
|
|
260
|
+
|
|
261
|
+
opendir (DIR, "./DataBase") or die ("error:$!");
|
|
262
|
+
|
|
263
|
+
my @database = readdir DIR;
|
|
264
|
+
|
|
265
|
+
my $dbcheck = 0;
|
|
266
|
+
|
|
267
|
+
foreach(@database){
|
|
268
|
+
|
|
269
|
+
if ($_ =~ /${db}.udb/){$dbcheck++;}
|
|
270
|
+
|
|
271
|
+
}
|
|
272
|
+
|
|
273
|
+
closedir DIR;
|
|
274
|
+
|
|
275
|
+
|
|
276
|
+
|
|
277
|
+
opendir (DIR, "./Tools") or die ("error:$!");
|
|
278
|
+
|
|
279
|
+
my @tool = readdir DIR;
|
|
280
|
+
|
|
281
|
+
my $usearch;
|
|
282
|
+
|
|
283
|
+
foreach (@tool) {
|
|
284
|
+
|
|
285
|
+
if ($_ =~ /(usearch.+)/){$usearch = $1;}
|
|
286
|
+
|
|
287
|
+
}
|
|
288
|
+
|
|
289
|
+
closedir DIR;
|
|
290
|
+
|
|
291
|
+
unless($usearch){print "Error!! Please put an usearch exectable file in Tools directory!\n"; exit;}
|
|
292
|
+
|
|
293
|
+
|
|
294
|
+
|
|
295
|
+
unless($dbcheck){
|
|
296
|
+
|
|
297
|
+
my $command = "./Tools/$usearch -makeudb_usearch \"./DataBase/$db\" -output \"./DataBase/${db}.udb\"";
|
|
298
|
+
|
|
299
|
+
system $command;
|
|
300
|
+
|
|
301
|
+
}
|
|
302
|
+
|
|
303
|
+
|
|
304
|
+
|
|
305
|
+
|
|
306
|
+
|
|
307
|
+
#Decompressing files
|
|
308
|
+
|
|
309
|
+
opendir (DIR, "./Run") or die ("error:$!");
|
|
310
|
+
|
|
311
|
+
my @run = readdir DIR;
|
|
312
|
+
|
|
313
|
+
my $gzip = 0;
|
|
314
|
+
|
|
315
|
+
foreach (@run) {
|
|
316
|
+
|
|
317
|
+
if ($_ =~ /.gz$/){$gzip++;}
|
|
318
|
+
|
|
319
|
+
}
|
|
320
|
+
|
|
321
|
+
closedir DIR;
|
|
322
|
+
|
|
323
|
+
|
|
324
|
+
|
|
325
|
+
if($gzip){
|
|
326
|
+
|
|
327
|
+
print "Decompressing gz file to fastq file...\n";
|
|
328
|
+
|
|
329
|
+
foreach (@run) {
|
|
330
|
+
|
|
331
|
+
if ($_ =~ /.gz$/){
|
|
332
|
+
|
|
333
|
+
my $comand = "gzip -d \"./Run/$_\"";
|
|
334
|
+
|
|
335
|
+
system $comand;
|
|
336
|
+
|
|
337
|
+
}
|
|
338
|
+
|
|
339
|
+
}
|
|
340
|
+
|
|
341
|
+
}
|
|
342
|
+
|
|
343
|
+
|
|
344
|
+
|
|
345
|
+
#Get data names
|
|
346
|
+
|
|
347
|
+
opendir (DIR, "./Run") or die ("error:$!");
|
|
348
|
+
|
|
349
|
+
my @read = readdir DIR;
|
|
350
|
+
|
|
351
|
+
my $read_count = 0;
|
|
352
|
+
|
|
353
|
+
my %file;
|
|
354
|
+
|
|
355
|
+
foreach (@read) {
|
|
356
|
+
|
|
357
|
+
if ($_ =~ /(.+)_R1/){$file{$1}++; $read_count++;}
|
|
358
|
+
|
|
359
|
+
if ($_ =~ /(.+)_R2/){$file{$1}++; $read_count++;}
|
|
360
|
+
|
|
361
|
+
}
|
|
362
|
+
|
|
363
|
+
closedir DIR;
|
|
364
|
+
|
|
365
|
+
unless(%file){print "None of fastq files in Run file!!\n"; exit;}
|
|
366
|
+
|
|
367
|
+
foreach(keys %file){
|
|
368
|
+
|
|
369
|
+
if($file{$_} < 2){print "Error!! Check the Run directory!The analysis needs both R1 and R2 fastq files.\nR1 or R2 fastq files of $_ is not.\n"; exit;}
|
|
370
|
+
|
|
371
|
+
}
|
|
372
|
+
|
|
373
|
+
unless(-d "./Results"){mkdir "./Results";}
|
|
374
|
+
|
|
375
|
+
|
|
376
|
+
|
|
377
|
+
|
|
378
|
+
|
|
379
|
+
print "============================================================\n";
|
|
380
|
+
|
|
381
|
+
print " 1_Preprocessing \n";
|
|
382
|
+
|
|
383
|
+
print "============================================================\n";
|
|
384
|
+
|
|
385
|
+
#assemble pair seqs
|
|
386
|
+
|
|
387
|
+
mkdir "./Results/1_1_Merge_paird_reads";
|
|
388
|
+
|
|
389
|
+
unless($primer =~ /^no$/i){mkdir "./Results/1_2_Strip_primers";}
|
|
390
|
+
|
|
391
|
+
mkdir "./Results/1_3_Quality_filter";
|
|
392
|
+
|
|
393
|
+
open (LOG, ">", "./Results/log.txt") or die("error:$!");
|
|
394
|
+
|
|
395
|
+
my (%log, %raw);
|
|
396
|
+
|
|
397
|
+
|
|
398
|
+
|
|
399
|
+
my @youbi = ('Sun', 'Mon', 'Tue', 'Wed', 'Thu', 'Fri', 'Sat');
|
|
400
|
+
|
|
401
|
+
my ($sec,$min,$hour,$mday,$mon,$year,$wday,$yday,$isdst) = localtime(time);
|
|
402
|
+
|
|
403
|
+
$year += 1900;
|
|
404
|
+
|
|
405
|
+
$mon += 1;
|
|
406
|
+
|
|
407
|
+
if($min < 10){$min = "0$min";}
|
|
408
|
+
|
|
409
|
+
print LOG "Start: $year/$mon/$mday ($youbi[$wday]) $hour:$min\n";
|
|
410
|
+
|
|
411
|
+
|
|
412
|
+
|
|
413
|
+
以下略
|
|
68
414
|
|
|
69
415
|
|
|
70
416
|
|
2
画像追加
test
CHANGED
|
File without changes
|
test
CHANGED
|
@@ -1,4 +1,4 @@
|
|
|
1
|
-
](848a977b10b17ad124d14d2a8e66bf0e.png)### 前提・実現したいこと
|
|
1
|
+
](848a977b10b17ad124d14d2a8e66bf0e.png)### 前提・実現したいこと
|
|
2
2
|
|
|
3
3
|
|
|
4
4
|
|
1
画像、コード追加
test
CHANGED
|
File without changes
|
test
CHANGED
|
@@ -1,4 +1,4 @@
|
|
|
1
|
-
### 前提・実現したいこと
|
|
1
|
+
](848a977b10b17ad124d14d2a8e66bf0e.png)### 前提・実現したいこと
|
|
2
2
|
|
|
3
3
|
|
|
4
4
|
|
|
@@ -28,6 +28,42 @@
|
|
|
28
28
|
|
|
29
29
|
perl PMiFish.pl
|
|
30
30
|
|
|
31
|
+
|
|
32
|
+
|
|
33
|
+
以下PMiFish.plコード
|
|
34
|
+
|
|
35
|
+
#!/usr/bin/perl
|
|
36
|
+
|
|
37
|
+
use strict;
|
|
38
|
+
|
|
39
|
+
use warnings;
|
|
40
|
+
|
|
41
|
+
|
|
42
|
+
|
|
43
|
+
system "perl ./Scripts/1_Preprocessing.pl";
|
|
44
|
+
|
|
45
|
+
system "perl ./Scripts/1_4_Rarefaction.pl";
|
|
46
|
+
|
|
47
|
+
system "perl ./Scripts/2_1_Find_unique.pl";
|
|
48
|
+
|
|
49
|
+
system "perl ./Scripts/2_2_Denoise.pl";
|
|
50
|
+
|
|
51
|
+
system "perl ./Scripts/2_3_Separate_chimera.pl";
|
|
52
|
+
|
|
53
|
+
system "perl ./Scripts/3_1_Usearch_global.pl";
|
|
54
|
+
|
|
55
|
+
system "perl ./Scripts/4_1_Annotation.pl";
|
|
56
|
+
|
|
57
|
+
system "perl ./Scripts/5_1_Fasta_for_Phylogenetic_Analysis.pl";
|
|
58
|
+
|
|
59
|
+
system "perl ./Scripts/5_2_Summary_table.pl";
|
|
60
|
+
|
|
61
|
+
system "perl ./Scripts/5_3_Fasta_classified_by_family.pl";
|
|
62
|
+
|
|
63
|
+
system "perl ./Scripts/6_1_Portal.pl";
|
|
64
|
+
|
|
65
|
+
system "perl ./Scripts/6_2_Phylogenetic_trees.pl";
|
|
66
|
+
|
|
31
67
|
```
|
|
32
68
|
|
|
33
69
|
|