Mercurial > repos > mgarnier > catchsequence

comparison catchsequence/catchsequence.pl @ 0:c30eb2050ef5 draft default tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author mgarnier
date Thu, 26 Aug 2021 13:41:33 +0000
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:c30eb2050ef5
1 #!/usr/bin/perl
2
3 use strict;
4 use warnings;
5
6 #INPUTS_
7 # my $Result_RES = $ARGV[0];
8 my $sequences = $ARGV[0];
9
10 #OUTPUT_
11 #my $output = $ARGV[1];
12
13 #my @list_seq = split(/,/,$sequences);
14 my @list_seq = @ARGV;
15
16 ##########################################################################################
17
18 #open (OUT, ">$output");
19 print "SEQUENCE\tRESISTANCE GENES\tPLASMIDS\tVIRULENCE GENES\tMLST NUMBER\tGENE(ALLELE)\n";
20
21 foreach my $sequence (@list_seq) {
22 my $Result_RES = `abricate --db resfinder $sequence > $sequence.RES.txt`; #appel système de la commande abricate avec la BDD ResFinder
23 my $Result_PLA = `abricate --db plasmidfinder $sequence > $sequence.PLA.txt`; #appel système de la commande abricate avec la BDD PlasmidFinder
24 my $Result_VIR = `abricate --db vfdb $sequence > $sequence.VIR.txt`;
25 my $Result_MLST = `mlst $sequence > $sequence.MLST.txt`;
26
27 open (RES, "$sequence.RES.txt");
28 print "$sequence\t";
29
30 while (<RES>) {
31
32 chomp();
33 if ($_ !~ m/^#/) {
34 my @infos = split(/\t/,$_);
35 my $geneRes = $infos[5]; # resistance gene name (ancienne valeur $infos[4])
36 my $identity = $infos[10]; # identity % (ancienne valeur $infos[9])
37
38 if ($identity > 90.000) {
39 print "$geneRes;";
40 }
41 }
42
43 }
44
45 close (RES);
46 print "\t";
47
48
49 open (PLA, "$sequence.PLA.txt") or die "could not open $!";
50
51 while (<PLA>) {
52 chomp();
53 if ($_ !~ m/^#/) {
54 my @infos = split(/\t/,$_);
55 my $plasmid = $infos[5]; # plasmid name
56 my $identity = $infos[10]; # identity %
57
58 if ($identity > 90.000) {
59 print"$plasmid;";
60 }
61 }
62
63 }
64 close (PLA);
65 print "\t";
66
67 open (VIR, "$sequence.VIR.txt") or die "could not open $!";
68
69 while (<VIR>) {
70 chomp();
71
72 if ($_ !~ m/^#/) {
73 my @infos = split(/\t/,$_);
74 my $geneVir = $infos[5]; # virulence gene name
75 my $identity = $infos[10]; # identity %
76
77 if ($identity > 80.000) {
78 print "$geneVir;";
79 }
80 }
81
82 }
83 close (VIR);
84 print "\t";
85
86
87 open (MLST, "$sequence.MLST.txt") or die "could not open $!";
88
89 while (<MLST>) {
90 chomp();
91
92 my @infos = split(/\t/,$_);
93 my $numMLST = $infos[2];
94 print "$numMLST\t";
95
96 for (my $i=3; $i <= $#infos; $i++){
97 print "$infos[$i];";
98 }
99
100 }
101 close (MLST);
102 print "\n";
103 }
104
105 #close (OUT);
106
107 unlink glob ('*.VIR.txt');
108 unlink glob ('*.PLA.txt');
109 unlink glob ('*.RES.txt');
110 unlink glob ('*.MLST.txt');
111