Mercurial > repos > dcouvin > catchsequenceinfo

comparison catchsequence/catchsequence.pl @ 0:37d48392bf22 draft default tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author dcouvin
date Tue, 21 Sep 2021 16:44:26 +0000
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:37d48392bf22
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 my $res = 90.00;
17 my $plas = 90.00;
18 my $vf = 80.00;
19
20 my $percentage = "ident"; # other possibility is "cov"
21 my $columPerc = 10; # other possibility is 9
22
23
24
25 #Other parameters
26 for (my $i = 0; $i <= $#ARGV; $i++) {
27 if ($ARGV[$i]=~/-percent/i or $ARGV[$i]=~/-perc/i) {
28 $percentage = $ARGV[$i+1];
29 }
30 elsif ($ARGV[$i]=~/-res/i) {
31 $res = $ARGV[$i+1];
32 }
33 elsif ($ARGV[$i]=~/-plas/i) {
34 $plas = $ARGV[$i+1];
35 }
36 elsif ($ARGV[$i]=~/-vf/i) {
37 $vf = $ARGV[$i+1];
38 }
39 }
40
41 ##########################################################################################
42
43 if ($percentage eq "ident"){
44 $columPerc = 10;
45 }
46 elsif ($percentage eq "cov"){
47 $columPerc = 9;
48 }
49
50
51 #open (OUT, ">$output");
52 print "Sequence\tResistance genes\tPlasmids\tVirulence genes\tST (MLST)\tAlleles (MLST)\n";
53
54 foreach my $sequence (@list_seq) {
55 my $Result_RES = `abricate --db resfinder $sequence > $sequence.RES.txt`; #appel système de la commande abricate avec la BDD ResFinder
56 my $Result_PLA = `abricate --db plasmidfinder $sequence > $sequence.PLA.txt`; #appel système de la commande abricate avec la BDD PlasmidFinder
57 my $Result_VIR = `abricate --db vfdb $sequence > $sequence.VIR.txt`;
58 my $Result_MLST = `mlst $sequence > $sequence.MLST.txt`;
59
60 open (RES, "$sequence.RES.txt");
61 print "$sequence\t";
62
63 while (<RES>) {
64
65 chomp();
66 if ($_ !~ m/^#/) {
67 my @infos = split(/\t/,$_);
68 my $geneRes = $infos[5]; # resistance gene name (ancienne valeur $infos[4])
69 my $identity = $infos[$columPerc]; # identity % (ancienne valeur $infos[9])
70
71 if ($identity > $res) {
72 print "$geneRes;";
73 }
74 }
75
76 }
77
78 close (RES);
79 print "\t";
80
81
82 open (PLA, "$sequence.PLA.txt") or die "could not open $!";
83
84 while (<PLA>) {
85 chomp();
86 if ($_ !~ m/^#/) {
87 my @infos = split(/\t/,$_);
88 my $plasmid = $infos[5]; # plasmid name
89 my $identity = $infos[$columPerc]; # identity %
90
91 if ($identity > $plas) {
92 print"$plasmid;";
93 }
94 }
95
96 }
97 close (PLA);
98 print "\t";
99
100 open (VIR, "$sequence.VIR.txt") or die "could not open $!";
101
102 while (<VIR>) {
103 chomp();
104
105 if ($_ !~ m/^#/) {
106 my @infos = split(/\t/,$_);
107 my $geneVir = $infos[5]; # virulence gene name
108 my $identity = $infos[$columPerc]; # identity %
109
110 if ($identity > $vf) {
111 print "$geneVir;";
112 }
113 }
114
115 }
116 close (VIR);
117 print "\t";
118
119
120 open (MLST, "$sequence.MLST.txt") or die "could not open $!";
121
122 while (<MLST>) {
123 chomp();
124
125 my @infos = split(/\t/,$_);
126 my $numMLST = $infos[2];
127 print "$numMLST\t";
128
129 for (my $i=3; $i <= $#infos; $i++){
130 print "$infos[$i];";
131 }
132
133 }
134 close (MLST);
135 print "\n";
136 }
137
138 #close (OUT);
139
140 unlink glob ('*.VIR.txt');
141 unlink glob ('*.PLA.txt');
142 unlink glob ('*.RES.txt');
143 unlink glob ('*.MLST.txt');
144