Mercurial > repos > kkonganti > cfsan_bettercallsal

annotate 0.6.1/assets/abricate-get_db @ 14:b0a37e88ecb5

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
"planemo upload"
author kkonganti
date Thu, 07 Sep 2023 10:13:31 -0400
parents 749faef1caa9
children
rev   line source
kkonganti@11 1 #!/usr/bin/env perl
kkonganti@11 2
kkonganti@11 3 use strict;
kkonganti@11 4 use FindBin;
kkonganti@11 5 use Bio::SeqIO;
kkonganti@11 6 use Bio::Seq;
kkonganti@11 7 use Path::Tiny;
kkonganti@11 8 use File::Basename;
kkonganti@11 9 use File::Spec;
kkonganti@11 10 use File::Path qw(make_path remove_tree);
kkonganti@11 11 use List::Util qw(first);
kkonganti@11 12 use Cwd qw(abs_path);
kkonganti@11 13 use Data::Dumper;
kkonganti@11 14 use LWP::Simple;
kkonganti@11 15 use JSON;
kkonganti@11 16
kkonganti@11 17 #..............................................................................
kkonganti@11 18 # Globals
kkonganti@11 19
kkonganti@11 20 my $EXE = basename($0);
kkonganti@11 21 my $ABX_SEP = ';';
kkonganti@11 22
kkonganti@11 23 my %DATABASE = (
kkonganti@11 24 'resfinder' => \&get_resfinder,
kkonganti@11 25 'plasmidfinder' => \&get_plasmidfinder,
kkonganti@11 26 'megares' => \&get_megares,
kkonganti@11 27 'argannot' => \&get_argannot,
kkonganti@11 28 'card' => \&get_card,
kkonganti@11 29
kkonganti@11 30 # 'ncbibetalactamase' => \&get_ncbibetalactamase,
kkonganti@11 31 'ncbi' => \&get_ncbi,
kkonganti@11 32 'vfdb' => \&get_vfdb,
kkonganti@11 33 'ecoli_vf' => \&get_ecoli_vf, # https://github.com/phac-nml/ecoli_vf
kkonganti@11 34 'ecoh' => \&get_ecoh,
kkonganti@11 35 'bacmet2' => \&get_bacmet2,
kkonganti@11 36 'victors' => \&get_victors,
kkonganti@11 37
kkonganti@11 38 # 'serotypefinder' => \&get_serotypefinder,
kkonganti@11 39 );
kkonganti@11 40 my $DATABASES = join( ' ', sort keys %DATABASE );
kkonganti@11 41
kkonganti@11 42 #..............................................................................
kkonganti@11 43 # Command line options
kkonganti@11 44
kkonganti@11 45 my ( @Options, $debug, $outdir, $db, $force );
kkonganti@11 46 setOptions();
kkonganti@11 47
kkonganti@11 48 $db or err("Please choose a --db from: $DATABASES");
kkonganti@11 49 exists $DATABASE{$db} or err("Unknown --db '$db', choose from: $DATABASES ");
kkonganti@11 50 -d $outdir or err("--outdir '$outdir' does not exist");
kkonganti@11 51
kkonganti@11 52 my $dir = abs_path( File::Spec->catdir( $outdir, $db ) );
kkonganti@11 53 make_path($dir);
kkonganti@11 54 msg("Setting up '$db' in '$dir'");
kkonganti@11 55
kkonganti@11 56 #my $tmpdir = tempdir("$db-XXXXXXXX", DIR=>$dir, CLEANUP=>0);
kkonganti@11 57 #my $tmpdir = "/home/tseemann/git/abricate/db/resfinder/resfinder-6Kuphtvv";
kkonganti@11 58 my $tmpdir = "$dir/src";
kkonganti@11 59 make_path($tmpdir);
kkonganti@11 60
kkonganti@11 61 # run the specific function from --db
kkonganti@11 62 chdir $tmpdir;
kkonganti@11 63 my $seq = $DATABASE{$db}->();
kkonganti@11 64 map { is_full_gene($_) } @$seq; # doesn't do anything?
kkonganti@11 65 $seq = dedupe_seq($seq);
kkonganti@11 66
kkonganti@11 67 #print Dumper($seq);
kkonganti@11 68 msg("Sorting sequences by ID");
kkonganti@11 69 $seq = [ sort { $a->{ID} cmp $b->{ID} } @$seq ];
kkonganti@11 70 save_fasta( "$dir/sequences", $seq );
kkonganti@11 71
kkonganti@11 72 msg("Formatting BLASTN database: $dir/sequences");
kkonganti@11 73 my $logfile = "$tmpdir/makeblastdb.log";
kkonganti@11 74 my $ncbi_title = $db;
kkonganti@11 75 if ( "$db" eq "ncbi" ) {
kkonganti@11 76 $ncbi_title = "ncbiamrplus";
kkonganti@11 77 }
kkonganti@11 78 my $ec = system(
kkonganti@11 79 "makeblastdb -in '$dir/sequences' -title '$ncbi_title' -dbtype nucl -hash_index -logfile $logfile"
kkonganti@11 80 );
kkonganti@11 81 if ( $ec != 0 ) {
kkonganti@11 82 system("tail '$logfile'");
kkonganti@11 83 err("Error with makign BLAST database. See $logfile");
kkonganti@11 84 }
kkonganti@11 85
kkonganti@11 86 #msg("Run 'abricate --setupdb' to format the database");
kkonganti@11 87
kkonganti@11 88 msg("Done.");
kkonganti@11 89
kkonganti@11 90 #..............................................................................
kkonganti@11 91
kkonganti@11 92 sub download {
kkonganti@11 93 my ( $url, $dest ) = @_;
kkonganti@11 94 if ( -r $dest and not $force ) {
kkonganti@11 95 msg("Won't re-download existing $dest (use --force)");
kkonganti@11 96
kkonganti@11 97 #exit(1);
kkonganti@11 98 }
kkonganti@11 99 else {
kkonganti@11 100 msg("Downloading: $url");
kkonganti@11 101 my $ec = mirror( $url, $dest );
kkonganti@11 102 msg("HTTP Result: $ec");
kkonganti@11 103 ( $ec == 200 or $ec = 304 )
kkonganti@11 104 or err("HTTP $ec | failed to download $url"); # is HTTP OK ?
kkonganti@11 105 }
kkonganti@11 106 msg("Destination: $dest");
kkonganti@11 107 msg( "Filesize:", ( -s $dest ), "bytes" );
kkonganti@11 108 }
kkonganti@11 109
kkonganti@11 110 #..............................................................................
kkonganti@11 111 sub trim_spaces {
kkonganti@11 112 my ($s) = @_;
kkonganti@11 113 $s =~ s/^\s+//;
kkonganti@11 114 $s =~ s/\s+$//;
kkonganti@11 115 return $s;
kkonganti@11 116 }
kkonganti@11 117
kkonganti@11 118 #..............................................................................
kkonganti@11 119 sub get_resfinder {
kkonganti@11 120 my $name = "resfinder_db";
kkonganti@11 121
kkonganti@11 122 # FIXME - can we just get HEAD.zip like in plasmidfinder?
kkonganti@11 123 my $url = "https://bitbucket.org/genomicepidemiology/$name.git";
kkonganti@11 124
kkonganti@11 125 # if (-r $name and not $force) {
kkonganti@11 126 # msg("Won't overwrite existing $name (use --force)");
kkonganti@11 127 # # exit(1);
kkonganti@11 128 # }
kkonganti@11 129 # else {
kkonganti@11 130 # msg("Nuking existing folder: $name");
kkonganti@11 131 # remove_tree("./$name");
kkonganti@11 132 # msg("Cloning $url to $name");
kkonganti@11 133 # system("git clone --quiet $url $name");
kkonganti@11 134 # }
kkonganti@11 135
kkonganti@11 136 #<*.fsa>
kkonganti@11 137 #>aac(6')-Ib_2_M23634
kkonganti@11 138 #>blaNDM-19_1_MF370080
kkonganti@11 139 #>mcr-1.1_1_KP347127
kkonganti@11 140 #>fosB1_1_CP001903
kkonganti@11 141 #>fusB_1_AY373761
kkonganti@11 142 #>VanHAX_1_FJ866609
kkonganti@11 143 #>ere(A)_6_DQ157752
kkonganti@11 144 #>nimA_1_X71444
kkonganti@11 145 #>cfr_1_AM408573
kkonganti@11 146 #>catB3_2_U13880
kkonganti@11 147 #>qnrA1_1_AY070235
kkonganti@11 148 #>ARR-2_1_HQ141279
kkonganti@11 149 #>sul1_2_U12338
kkonganti@11 150 #>tet_1_M74049
kkonganti@11 151 #>dfrA19_1_EU855687
kkonganti@11 152
kkonganti@11 153 #<notes.txt>
kkonganti@11 154 #aac(6')-Iv:Aminoglycoside resistance:
kkonganti@11 155 #aac(6')-Iw:Aminoglycoside resistance:Alternate name; aac(6')-Ix
kkonganti@11 156 #sul3:Sulphonamide resistance:
kkonganti@11 157 ##Tetracycline:
kkonganti@11 158 #ort(B):Tetracycline resistance:
kkonganti@11 159 #blaCMY-59:Beta-lactam resistance:
kkonganti@11 160
kkonganti@11 161 #<phenotypes.txt>
kkonganti@11 162 #Gene_accession no. Class Phenotype PMID Mechanism of resistance Notes Required_gene
kkonganti@11 163 #ant(2'')-Ia_1_X04555 Aminoglycoside Gentamicin, Tobramycin 3024112 Enzymatic modification Alternative name aadB
kkonganti@11 164 #ant(2'')-Ia_2_JF826500 Aminoglycoside Gentamicin, Tobramycin 22271862 Enzymat
kkonganti@11 165
kkonganti@11 166 $name = "~/apps/bettercallsal/assets/abricate_dbs/$name";
kkonganti@11 167 my $metafn = "$name/phenotypes.txt";
kkonganti@11 168 my @meta = path($metafn)->lines( { chomp => 1 } );
kkonganti@11 169 my %anno;
kkonganti@11 170 foreach (@meta) {
kkonganti@11 171 next if m/^#/;
kkonganti@11 172 my @x = split m/\t/;
kkonganti@11 173
kkonganti@11 174 #msg("$metafn: @x");
kkonganti@11 175 my ($gene) = ( $x[0] =~ m/^(.*?)_\w+$/ );
kkonganti@11 176 $anno{$gene}{ABX} = [
kkonganti@11 177 map { trim_spaces($_) }
kkonganti@11 178 grep { !m/(unknown|notes|^none)/i }
kkonganti@11 179 split m/,\s*/,
kkonganti@11 180 $x[2]
kkonganti@11 181 ];
kkonganti@11 182
kkonganti@11 183 #msg("$metafn: $gene |", $anno{$gene}{ABX}->@*);
kkonganti@11 184 }
kkonganti@11 185 msg( "get_resfinder: $metafn", scalar( keys %anno ), "genes" );
kkonganti@11 186
kkonganti@11 187 #print Dumper(\%anno);
kkonganti@11 188
kkonganti@11 189 my @seq;
kkonganti@11 190 for my $fasta (<$name/*.fsa>) {
kkonganti@11 191
kkonganti@11 192 # Issue #62 - repair broken fasta files like this:
kkonganti@11 193 # GCTTTAAATTGGAAAAAAGATAGTCAAACTCTTTAA>cmr_1_U43535
kkonganti@11 194 # inline replacement
kkonganti@11 195 system( 'sed', '-i.bak', 's/\([A-Z]\)>/\1\n>/gi', $fasta );
kkonganti@11 196 my $args = load_fasta($fasta);
kkonganti@11 197
kkonganti@11 198 # use name of fasta file as antibiotic name
kkonganti@11 199 #my $abx = basename($fasta, '.fsa');
kkonganti@11 200 #msg("$fasta: Assigning '$abx' to all genes");
kkonganti@11 201 #push @{$_->{ABX}}, $abx for (@$args);
kkonganti@11 202 push @seq, @$args;
kkonganti@11 203 }
kkonganti@11 204
kkonganti@11 205 # https://github.com/tseemann/abricate/issues/92
kkonganti@11 206 # mcr-9_1_NZ_NAAN01000063.1
kkonganti@11 207 #>mcr-9_1_NZ_NAAN01000063.1
kkonganti@11 208 # mcr-9.1:Colistin resistance:
kkonganti@11 209
kkonganti@11 210 for my $seq (@seq) {
kkonganti@11 211 my ( $id, $copy, $acc ) = $seq->{ID} =~ m/^(.*?)_(\d+)_(\S+)$/;
kkonganti@11 212
kkonganti@11 213 #msg("resfinder: $1 $2 $3", $anno{$1});
kkonganti@11 214 $seq->{ID} = "${id}_${copy}";
kkonganti@11 215 $seq->{ACC} = $acc;
kkonganti@11 216 $seq->{DESC} = $anno{$id}{DESC} || $id;
kkonganti@11 217 push @{ $seq->{ABX} }, @{ $anno{$id}{ABX} } if $anno{$id}{ABX};
kkonganti@11 218 }
kkonganti@11 219
kkonganti@11 220 return \@seq;
kkonganti@11 221 }
kkonganti@11 222
kkonganti@11 223 #..............................................................................
kkonganti@11 224 sub get_serotypefinder {
kkonganti@11 225 my $name = "serotypefinder_db";
kkonganti@11 226 my $url = "https://bitbucket.org/genomicepidemiology/$name.git";
kkonganti@11 227
kkonganti@11 228 if ( -r $name and not $force ) {
kkonganti@11 229 msg("Won't overwrite existing $name (use --force)");
kkonganti@11 230
kkonganti@11 231 # exit(1);
kkonganti@11 232 }
kkonganti@11 233 else {
kkonganti@11 234 msg("Nuking existing folder: $name");
kkonganti@11 235 remove_tree("./$name");
kkonganti@11 236 msg("Cloning $url to $name");
kkonganti@11 237 system("git clone --quiet $url $name");
kkonganti@11 238 }
kkonganti@11 239
kkonganti@11 240 my @seq;
kkonganti@11 241 for my $fasta (<$name/*.fsa>) {
kkonganti@11 242 push @seq, @{ load_fasta($fasta) };
kkonganti@11 243 }
kkonganti@11 244
kkonganti@11 245 # >fliC_44444_AY250028_H52
kkonganti@11 246 # FIXME - this is already in EcOH database!
kkonganti@11 247 for my $seq (@seq) {
kkonganti@11 248 my ( $id, $copy, $acc ) = $seq->{ID} =~ m/^(.*)_(\d+)_(\w+)$/;
kkonganti@11 249
kkonganti@11 250 #msg("serotypefinder: $1 $2 $3", $anno{$1});
kkonganti@11 251 $seq->{ID} = "${id}_${copy}";
kkonganti@11 252 $seq->{ACC} = $acc;
kkonganti@11 253
kkonganti@11 254 #$seq->{DESC} = $anno{$id} || '';
kkonganti@11 255 }
kkonganti@11 256
kkonganti@11 257 return \@seq;
kkonganti@11 258 }
kkonganti@11 259
kkonganti@11 260 #..............................................................................
kkonganti@11 261 sub get_tag {
kkonganti@11 262 my ( $f, $tag ) = @_;
kkonganti@11 263 if ( $f->has_tag($tag) ) {
kkonganti@11 264 my ($val) = $f->get_tag_values($tag);
kkonganti@11 265 return $val;
kkonganti@11 266 }
kkonganti@11 267 return '';
kkonganti@11 268 }
kkonganti@11 269
kkonganti@11 270 #..............................................................................
kkonganti@11 271 sub get_ncbi {
kkonganti@11 272 my $AFP =
kkonganti@11 273 "https://ftp.ncbi.nlm.nih.gov/pathogen/Antimicrobial_resistance/AMRFinderPlus/database/latest";
kkonganti@11 274
kkonganti@11 275 #my $src = "https://ftp.ncbi.nlm.nih.gov/pathogen/Antimicrobial_resistance/AMRFinderPlus/data/latest/AMR_CDS";
kkonganti@11 276 my $src = "$AFP/AMR_CDS";
kkonganti@11 277 my $name = "amr_cds.ffn";
kkonganti@11 278
kkonganti@11 279 #my $src2 = "https://ftp.ncbi.nlm.nih.gov/pathogen/Antimicrobial_resistance/AMRFinderPlus/data/latest/ReferenceGeneCatalog.txt";
kkonganti@11 280 my $src2 = "$AFP/ReferenceGeneCatalog.txt";
kkonganti@11 281 my $name2 = "amr_cds.tsv";
kkonganti@11 282
kkonganti@11 283 if ( -r $name and -r $name2 and not $force ) {
kkonganti@11 284 msg("Won't overwrite existing $name/$name2 (use --force)");
kkonganti@11 285
kkonganti@11 286 # exit(1);
kkonganti@11 287 }
kkonganti@11 288 else {
kkonganti@11 289 download( $src, $name );
kkonganti@11 290 download( $src2, $name2 );
kkonganti@11 291 }
kkonganti@11 292
kkonganti@11 293 #1 allele
kkonganti@11 294 #2 gene_family ble
kkonganti@11 295 #3 whitelisted_taxa
kkonganti@11 296 #4 product_name BLMA family bleomycin binding protein
kkonganti@11 297 #5 scope core
kkonganti@11 298 #6 type AMR
kkonganti@11 299 #7 subtype AMR
kkonganti@11 300 #8 class BLEOMYCIN
kkonganti@11 301 #9 subclass BLEOMYCIN
kkonganti@11 302 #10 refseq_protein_accession WP_063842967.1
kkonganti@11 303 #11 refseq_nucleotide_accession NG_047554.1
kkonganti@11 304 #12 curated_refseq_start No
kkonganti@11 305 #13 genbank_protein_accession CAA02068.1
kkonganti@11 306 #14 genbank_nucleotide_accession A31900.1
kkonganti@11 307 #15 genbank_strand_orientation +
kkonganti@11 308 #16 genbank_cds_start 6
kkonganti@11 309 #17 genbank_cds_stop 374
kkonganti@11 310 #18 pubmed_reference
kkonganti@11 311 #19 blacklisted_taxa
kkonganti@11 312 #20 db_version 2019-08-27.1
kkonganti@11 313
kkonganti@11 314 my $tsv = load_tabular( $name2, 10 ); # refseq_nucleotide_accession
kkonganti@11 315 msg( "[$name2] Loaded", scalar keys %$tsv, "records" );
kkonganti@11 316
kkonganti@11 317 # https://github.com/ncbi/amr/wiki/AMRFinderPlus-database#amrprot
kkonganti@11 318 # 0 1 2 3 4 5 6 7
kkonganti@11 319 # >1000909371|WP_061158039.1|NG_050200|1|1|blaTEM-156|blaTEM|class_A_beta-lactamase_TEM-156 NG_050200:101-961
kkonganti@11 320 my @seq;
kkonganti@11 321 my $in = Bio::SeqIO->new( -file => $name, -format => "fasta" );
kkonganti@11 322 while ( my $rec = $in->next_seq ) {
kkonganti@11 323
kkonganti@11 324 # parse ID
kkonganti@11 325 my ( $gi, $pi, $acc, $fp, $fn, $gene, $fam, $prod ) = split m/\|/,
kkonganti@11 326 $rec->id;
kkonganti@11 327
kkonganti@11 328 # skip fusion genes
kkonganti@11 329 next unless $fp == 1 and $fn == 1;
kkonganti@11 330
kkonganti@11 331 # only keep true ARGs
kkonganti@11 332 $acc .= ".1" unless $acc =~ m/\.\d+$/;
kkonganti@11 333 my $t = $tsv->{$acc} or next;
kkonganti@11 334 next unless ( $t->{scope} eq 'plus' || $t->{scope} eq 'core' );
kkonganti@11 335
kkonganti@11 336 # next unless $t->{type} eq 'VIRULENCE';
kkonganti@11 337 # next unless $t->{subtype} eq 'VIRULENCE';
kkonganti@11 338 # construct sequence record
kkonganti@11 339 $prod =~ s/_/ /g;
kkonganti@11 340 err("$pi: gene is empty") unless $gene;
kkonganti@11 341 err("$pi: product is empty") unless $prod;
kkonganti@11 342 my $s = {
kkonganti@11 343 ID => $gene,
kkonganti@11 344 ACC => $acc,
kkonganti@11 345 DESC => $prod,
kkonganti@11 346 SEQ => $rec->seq,
kkonganti@11 347 ABX => [ split m'/', $t->{subclass} ]
kkonganti@11 348 };
kkonganti@11 349 push @seq, $s;
kkonganti@11 350 msg( "[$name]", 0 + @seq, "|", $s->{ID}, "|", $s->{ACC}, "|",
kkonganti@11 351 $s->{DESC} );
kkonganti@11 352
kkonganti@11 353 #msg(Dumper($s));
kkonganti@11 354 msg( $s->{ID}, " is fusion $fp/$fn" ) if "$fp$fn" ne '11';
kkonganti@11 355 }
kkonganti@11 356 return \@seq;
kkonganti@11 357 }
kkonganti@11 358
kkonganti@11 359 #..............................................................................
kkonganti@11 360 sub get_plasmidfinder {
kkonganti@11 361 my $name = "plasmidfinder";
kkonganti@11 362 my $zip = "$name.zip";
kkonganti@11 363
kkonganti@11 364 # download("https://cge.cbs.dtu.dk/cge/download_data.php?folder=$name&filename=$zip&submit=$zip", $zip);
kkonganti@11 365 download(
kkonganti@11 366 "https://bitbucket.org/genomicepidemiology/plasmidfinder_db/get/HEAD.zip",
kkonganti@11 367 $zip
kkonganti@11 368 );
kkonganti@11 369 system("unzip -j -u $zip");
kkonganti@11 370
kkonganti@11 371 my @seq;
kkonganti@11 372 for my $fasta (<*.fsa>) {
kkonganti@11 373 push @seq, @{ load_fasta($fasta) };
kkonganti@11 374 }
kkonganti@11 375
kkonganti@11 376 for my $seq (@seq) {
kkonganti@11 377 $seq->{DESC} = $seq->{ID}; # no desc, so use ORIGINAL ID as desc
kkonganti@11 378 my ( $id, $acc ) =
kkonganti@11 379 ( $seq->{ID} =~ m/^(.*)_(([A-Z]+|NC_)\d+(\.\d+)?)$/ );
kkonganti@11 380 $id =~ s/_+$//g;
kkonganti@11 381 $seq->{ID} = $id || $seq->{ID};
kkonganti@11 382 $seq->{ACC} = $acc || '';
kkonganti@11 383 wrn( "Parsed empty ID:",
kkonganti@11 384 $seq->{DESC},
kkonganti@11 385 "=> id='$id' acc='$acc' seq=" . substr( $seq->{SEQ}, 0, 10 ) )
kkonganti@11 386 if not $id;
kkonganti@11 387 }
kkonganti@11 388
kkonganti@11 389 return \@seq;
kkonganti@11 390 }
kkonganti@11 391
kkonganti@11 392 #..............................................................................
kkonganti@11 393 sub get_megares {
kkonganti@11 394 my $zip = "megares.zip";
kkonganti@11 395 download( 'https://www.meglab.org/downloads/megares_v3.00.zip', $zip );
kkonganti@11 396 system("unzip -j -u $zip");
kkonganti@11 397 my $seqs = load_fasta( glob("megares_drugs_*.fasta") );
kkonganti@11 398 my @okseq;
kkonganti@11 399
kkonganti@11 400 # >MEG_372|Multi-compound|Drug_and_biocide_resistance|Drug_and_biocide_MATE_efflux_pumps|ABEM
kkonganti@11 401 # >MEG_411|Multi-compound|Drug_and_biocide_resistance|Drug_and_biocide_RND_efflux_regulator|ACRR|RequiresSNPConfirmation
kkonganti@11 402 # >MEG_7860|Drugs|betalactams|Class_B_betalactamases|ZOG
kkonganti@11 403 # >MEG_7439|Drugs|Glycopeptides|VanI-type_resistance_protein|VANI
kkonganti@11 404 # >MEG_7245|Drugs|Tetracyclines|Tetracycline_resistance_MFS_efflux_pumps|TETY
kkonganti@11 405 # >MEG_9|Drugs|Aminoglycosides|Aminoglycoside-resistant_16S_ribosomal_subunit_protein|A16S|RequiresSNPConfirmation
kkonganti@11 406
kkonganti@11 407 for my $s (@$seqs) {
kkonganti@11 408 my ( $id, $type, $class, $mech, $group, $note ) = split m/\|/, $s->{ID};
kkonganti@11 409 if ($note) {
kkonganti@11 410
kkonganti@11 411 # "RequiresSNPConfirmation" is the common one; we can't do that
kkonganti@11 412 msg("Skipping $id due to: $note");
kkonganti@11 413 next;
kkonganti@11 414 }
kkonganti@11 415 $s->{ID} = $group;
kkonganti@11 416 $s->{ACC} = $id;
kkonganti@11 417 $s->{DESC} = join( ':', $type, $class, $mech, $group );
kkonganti@11 418 push @okseq, $s;
kkonganti@11 419 }
kkonganti@11 420 return [@okseq];
kkonganti@11 421
kkonganti@11 422 #return $seqs;
kkonganti@11 423 }
kkonganti@11 424
kkonganti@11 425 #..............................................................................
kkonganti@11 426 sub get_argannot {
kkonganti@11 427 my $fasta = 'arg-annot.fa';
kkonganti@11 428 download(
kkonganti@11 429 # 'http://www.mediterranee-infection.com/arkotheque/client/ihumed/_depot_arko/articles/1425/argannot-aa-v3-march2017_doc.fasta',
kkonganti@11 430 # 'http://www.mediterranee-infection.com/arkotheque/client/ihumed/_depot_arko/articles/691/argannot-nt_doc.fasta',
kkonganti@11 431 'https://www.mediterranee-infection.com/wp-content/uploads/2019/06/ARG_ANNOT_V5_Nt_JUNE2019.txt',
kkonganti@11 432 $fasta
kkonganti@11 433 );
kkonganti@11 434
kkonganti@11 435 # fix syntax errors in the FASTA file...
kkonganti@11 436 path($fasta)->edit( sub { s/\\//g; $_ } );
kkonganti@11 437
kkonganti@11 438 my $seqs = load_fasta($fasta);
kkonganti@11 439
kkonganti@11 440 # 0 1 2 3
kkonganti@11 441 # >(AGly)Aac2-Ie:NC_011896:3039059-3039607:549
kkonganti@11 442 for my $s (@$seqs) {
kkonganti@11 443 my @x = split m/:/, $s->{ID};
kkonganti@11 444 $s->{ID} = $x[0];
kkonganti@11 445 $s->{ACC} = $x[1] . ':' . $x[2];
kkonganti@11 446 $s->{DESC} = '';
kkonganti@11 447 }
kkonganti@11 448
kkonganti@11 449 return $seqs;
kkonganti@11 450 }
kkonganti@11 451
kkonganti@11 452 #..............................................................................
kkonganti@11 453 sub get_bacmet2 {
kkonganti@11 454 my $fasta = 'bacmet2.fa';
kkonganti@11 455 download(
kkonganti@11 456 'http://bacmet.biomedicine.gu.se/download/BacMet2_EXP_database.fasta',
kkonganti@11 457 $fasta );
kkonganti@11 458
kkonganti@11 459 # This is a PROTEIN file
kkonganti@11 460 my $seqs = load_fasta($fasta);
kkonganti@11 461
kkonganti@11 462 # 0 1 2 3 4 ^
kkonganti@11 463 # >BAC0098|ctpC|sp|P0A502|CTPC_MYCTU Probable manganese/zinc-exporting
kkonganti@11 464 for my $s (@$seqs) {
kkonganti@11 465 my @x = split m/\|/, $s->{ID};
kkonganti@11 466 $s->{ID} = $x[1] . '-' . $x[0];
kkonganti@11 467 $s->{ACC} = $x[2] . ':' . $x[3];
kkonganti@11 468 }
kkonganti@11 469
kkonganti@11 470 return $seqs;
kkonganti@11 471 }
kkonganti@11 472
kkonganti@11 473 #..............................................................................
kkonganti@11 474 sub get_card {
kkonganti@11 475
kkonganti@11 476 # https://github.com/tseemann/abricate/issues/25
kkonganti@11 477 my $tarball = 'card.tar.bz2';
kkonganti@11 478 download(
kkonganti@11 479 #'https://card.mcmaster.ca/download/0/broadstreet-v2.0.2.tar.gz',
kkonganti@11 480 'https://card.mcmaster.ca/latest/data',
kkonganti@11 481
kkonganti@11 482 #'https://card.mcmaster.ca/latest/data/card-data.tar.bz2',
kkonganti@11 483 $tarball # yes, it's really BZ2 not GZ ...
kkonganti@11 484 );
kkonganti@11 485
kkonganti@11 486 # my $fasta = "./nucleotide_fasta_protein_homolog_model.fasta";
kkonganti@11 487 my $jsonfile = "./card.json";
kkonganti@11 488 system( "tar", "xf", $tarball, $jsonfile ) == 0
kkonganti@11 489 or err("Problem with tar xf $tarball $jsonfile");
kkonganti@11 490 -r $jsonfile or err("Could not extract $jsonfile from $tarball");
kkonganti@11 491
kkonganti@11 492 # JSON
kkonganti@11 493 my $json = path($jsonfile)->slurp_utf8;
kkonganti@11 494 my $card = from_json( $json, { latin1 => 1 } );
kkonganti@11 495 my @seq;
kkonganti@11 496 for my $g ( values %$card ) {
kkonganti@11 497 next unless ref($g) eq 'HASH';
kkonganti@11 498
kkonganti@11 499 # msg(Dumper($g));
kkonganti@11 500 next
kkonganti@11 501 unless $g->{model_type} eq
kkonganti@11 502 "protein homolog model"; # only 'acquired' genes
kkonganti@11 503 my $id = $g->{model_name};
kkonganti@11 504 err("$id has {model_param}{snp}") if exists $g->{model_param}{snp};
kkonganti@11 505
kkonganti@11 506 # msg("CARD: $id");
kkonganti@11 507 # print STDERR Dumper($g);
kkonganti@11 508 my $dna = $g->{model_sequences}{sequence}
kkonganti@11 509 or err("$id: no {model_sequences}{sequence} found");
kkonganti@11 510 my ($key) = sort keys %$dna; # first key
kkonganti@11 511 $dna = $dna->{$key} or err("$id: invalid key '$key'");
kkonganti@11 512 $dna = $dna->{dna_sequence} or err("$id: no dna_sequence");
kkonganti@11 513
kkonganti@11 514 # msg(Dumper($dna)) if $id eq 'OXA-25';
kkonganti@11 515
kkonganti@11 516 # ARO_category => {
kkonganti@11 517 # 'category_aro_name' => 'cephalosporin',
kkonganti@11 518 # 'category_aro_class_name' => 'Drug Class',
kkonganti@11 519 my $is_amr_gene = 0;
kkonganti@11 520 my @abx;
kkonganti@11 521 for my $key ( keys $g->{ARO_category}->%* ) {
kkonganti@11 522 my $c = $g->{ARO_category}{$key};
kkonganti@11 523 if ( $c->{category_aro_class_name} eq 'Drug Class' ) {
kkonganti@11 524 my $abx = $c->{category_aro_name};
kkonganti@11 525 $abx =~ s/ antibiotic//;
kkonganti@11 526 $abx =~ s/\s/_/g;
kkonganti@11 527 push @abx, $abx;
kkonganti@11 528 }
kkonganti@11 529 if ( $c->{category_aro_class_name} eq 'AMR Gene Family' ) {
kkonganti@11 530 $is_amr_gene++;
kkonganti@11 531 }
kkonganti@11 532 }
kkonganti@11 533
kkonganti@11 534 #err("CARD | $id | ", Dumper($g->{ARO_category}) ) unless $is_amr_gene;
kkonganti@11 535 #msg("ABX=$_") for @abx;
kkonganti@11 536
kkonganti@11 537 # put coordinates into normal form
kkonganti@11 538 my ( $start, $stop ) =
kkonganti@11 539 $dna->{strand} eq '-'
kkonganti@11 540 ? ( $dna->{fmax}, $dna->{fmin} )
kkonganti@11 541 : ( $dna->{fmin}, $dna->{fmax} );
kkonganti@11 542
kkonganti@11 543 $id =~ s/\s+/_/g;
kkonganti@11 544 push @seq,
kkonganti@11 545 {
kkonganti@11 546 ID => $id,
kkonganti@11 547 ACC => $dna->{accession} . ":$start-$stop",
kkonganti@11 548 DESC => ( $g->{ARO_description} || $g->{ARO_accession} ),
kkonganti@11 549 SEQ => $dna->{sequence},
kkonganti@11 550 ABX => [@abx],
kkonganti@11 551 };
kkonganti@11 552
kkonganti@11 553 # msg(Dumper($seq[-1]));
kkonganti@11 554 }
kkonganti@11 555
kkonganti@11 556 return \@seq;
kkonganti@11 557 }
kkonganti@11 558
kkonganti@11 559 #..............................................................................
kkonganti@11 560 sub get_victors {
kkonganti@11 561
kkonganti@11 562 # the CDS data is in .ffn and has source GI and coords
kkonganti@11 563 # the PROT data is in .faa and has the protein ref and /product
kkonganti@11 564
kkonganti@11 565 #>gi|115534241:2616-3152 Campylobacter jejuni plasmid pCJ01, complete sequence
kkonganti@11 566 download( 'http://www.phidias.us/victors/downloads/gen_downloads.php',
kkonganti@11 567 'victors.ffn' );
kkonganti@11 568
kkonganti@11 569 #>gi|115534244|ref|YP_783826.1| hypothetical protein pCJ01p4 [Campylobacter jejuni]
kkonganti@11 570 download(
kkonganti@11 571 'http://www.phidias.us/victors/downloads/gen_downloads_protein.php',
kkonganti@11 572 'victors.faa' );
kkonganti@11 573
kkonganti@11 574 my %gi;
kkonganti@11 575 open my $FAA, '<', 'victors.faa';
kkonganti@11 576 while (<$FAA>) {
kkonganti@11 577
kkonganti@11 578 #>gi|115534244|ref|YP_783826.1| hypothetical protein pCJ01p4 [Campylobacter jejuni]
kkonganti@11 579 next unless m"^>gi.(\d+).ref.([^|]+). ([^[]+)";
kkonganti@11 580 $gi{$1}{ACC} = $2;
kkonganti@11 581 $gi{$1}{DESC} = $3;
kkonganti@11 582 }
kkonganti@11 583
kkonganti@11 584 my $seqs = load_fasta("victors.ffn");
kkonganti@11 585
kkonganti@11 586 for my $s (@$seqs) {
kkonganti@11 587
kkonganti@11 588 #>gi|115534241:2616-3152 Campylobacter jejuni plasmid pCJ01, complete sequence
kkonganti@11 589 $s->{ID} =~ m/gi.(\d+):(\d+)-(\d+)/;
kkonganti@11 590 $s->{ACC} = $gi{$1}{ACC} || "gi|$1:$2-$3";
kkonganti@11 591 $s->{DESC} =~ $gi{$1}{DESC} || 'hypothetical protein';
kkonganti@11 592 }
kkonganti@11 593
kkonganti@11 594 # print Dumper($seqs); exit;
kkonganti@11 595
kkonganti@11 596 return $seqs;
kkonganti@11 597 }
kkonganti@11 598
kkonganti@11 599 #..............................................................................
kkonganti@11 600 sub get_vfdb {
kkonganti@11 601 download( 'http://www.mgc.ac.cn/VFs/Down/VFDB_setA_nt.fas.gz',
kkonganti@11 602 'vfdb.fa.gz' );
kkonganti@11 603 system("gzip -f -d -c vfdb.fa.gz > vfdb.fa");
kkonganti@11 604 my $seqs = load_fasta("vfdb.fa");
kkonganti@11 605
kkonganti@11 606 # >VFG000676(gb|AAD32411) (lef) anthrax toxin lethal factor precursor [Anthrax toxin (VF0142)] [Bacillus anthracis str. Sterne]
kkonganti@11 607 for my $s (@$seqs) {
kkonganti@11 608
kkonganti@11 609 # https://github.com/tseemann/abricate/issues/64#issuecomment-421895159 by @VGalata
kkonganti@11 610 $s->{ID} =~ m/^(\w+)\(\w+\|(\w+)(\.\d+)?\)$/; #
kkonganti@11 611 #$s->{ID} =~ m/^(\w+)\(\w+\|(\w+)\)$/;
kkonganti@11 612 $s->{ACC} = $2 if $2;
kkonganti@11 613 $s->{DESC} =~ m/^\((.*?)\)/;
kkonganti@11 614 $s->{ID} = $1 if $1;
kkonganti@11 615
kkonganti@11 616 # print STDERR Dumper($s); exit;
kkonganti@11 617 }
kkonganti@11 618
kkonganti@11 619 return $seqs;
kkonganti@11 620 }
kkonganti@11 621
kkonganti@11 622 #..............................................................................
kkonganti@11 623 sub get_ncbibetalactamase {
kkonganti@11 624 my $fasta = "ncbi.fa";
kkonganti@11 625 download(
kkonganti@11 626 'ftp://ftp.ncbi.nlm.nih.gov/pathogen/betalactamases/Allele-dna.fa',
kkonganti@11 627 $fasta );
kkonganti@11 628 my $tab = "ncbi.tab";
kkonganti@11 629 download( 'ftp://ftp.ncbi.nlm.nih.gov/pathogen/betalactamases/Allele.tab',
kkonganti@11 630 $tab );
kkonganti@11 631
kkonganti@11 632 # >ACD12694.1 EU650653.1:1-1173
kkonganti@11 633 my $seqs = load_fasta($fasta);
kkonganti@11 634
kkonganti@11 635 # ACC-1 ACD12694.1 EU650653.1 blaACC-1 1 1173 + cephalosporin-hydrolyzing class C beta-lactamase ACC-1
kkonganti@11 636 my %anno;
kkonganti@11 637 my @anno = grep { !m/^#/ } path($tab)->lines( { chomp => 1 } );
kkonganti@11 638 msg( "Read", 0 + @anno, "annotations" );
kkonganti@11 639 foreach (@anno) {
kkonganti@11 640 my ( $name, $id, $acc, $gene, $begin, $end, undef, $product ) =
kkonganti@11 641 split m/\t/;
kkonganti@11 642 $anno{$id} = {
kkonganti@11 643 ID => $gene,
kkonganti@11 644 DESC => $product,
kkonganti@11 645 ACC => "$acc:$begin-$end",
kkonganti@11 646 };
kkonganti@11 647 }
kkonganti@11 648
kkonganti@11 649 # print Dumper(\%anno);
kkonganti@11 650
kkonganti@11 651 for my $s (@$seqs) {
kkonganti@11 652 my $id = $s->{ID};
kkonganti@11 653 next unless exists $anno{$id};
kkonganti@11 654 $s->{ID} = $anno{$id}{ID};
kkonganti@11 655 $s->{ACC} = $anno{$id}{ACC};
kkonganti@11 656 $s->{DESC} = $anno{$id}{DESC};
kkonganti@11 657 }
kkonganti@11 658
kkonganti@11 659 # print Dumper($seqs);
kkonganti@11 660
kkonganti@11 661 return $seqs;
kkonganti@11 662 }
kkonganti@11 663
kkonganti@11 664 #..............................................................................
kkonganti@11 665 sub get_ecoh {
kkonganti@11 666 my $fasta = "EcOH.fa";
kkonganti@11 667 download(
kkonganti@11 668 'https://raw.githubusercontent.com/katholt/srst2/master/data/EcOH.fasta',
kkonganti@11 669 $fasta
kkonganti@11 670 );
kkonganti@11 671
kkonganti@11 672 # https://github.com/katholt/srst2#generating-srst2-compatible-clustered-database-from-raw-sequences
kkonganti@11 673 # [clusterID]__[gene]__[allele]__[seqID] [other stuff]
kkonganti@11 674 # >1__fliC__fliC-H1__1 AB028471.1;flagellin;H1
kkonganti@11 675 # >8__wzx__wzx-O41__246 AB811617.1;O antigen flippase;O41
kkonganti@11 676 # >9__wzy__wzy-OgN31__597 LC125932.1;O antigen polyermase;OgN31
kkonganti@11 677 my $seqs = load_fasta($fasta);
kkonganti@11 678
kkonganti@11 679 for my $s (@$seqs) {
kkonganti@11 680 my @id = split m/__/, $s->{ID};
kkonganti@11 681 my @desc = split m';', $s->{DESC};
kkonganti@11 682 $s->{ID} = $id[2];
kkonganti@11 683 $s->{ACC} = shift(@desc);
kkonganti@11 684 $s->{DESC} = join( ' ', @desc );
kkonganti@11 685 }
kkonganti@11 686
kkonganti@11 687 # print Dumper($seqs);
kkonganti@11 688 return $seqs;
kkonganti@11 689 }
kkonganti@11 690
kkonganti@11 691 #..............................................................................
kkonganti@11 692 sub get_ecoli_vf {
kkonganti@11 693 my $fasta = "ecoli_vf.ffn";
kkonganti@11 694 download(
kkonganti@11 695 'https://github.com/phac-nml/ecoli_vf/raw/master/data/repaired_ecoli_vfs_shortnames.ffn',
kkonganti@11 696 $fasta
kkonganti@11 697 );
kkonganti@11 698 my $seqs = load_fasta($fasta);
kkonganti@11 699
kkonganti@11 700 # >VFG000748(gi:2865308) (espF) EspF [EspF (VF0182)] [Escherichia coli O127:H6 str. E2348/69]
kkonganti@11 701 # >VFG000749(gi:6009379) (bfpA) Bundlin [BFP (VF0174)] [Escherichia coli B171]
kkonganti@11 702 # >SPG000142 (cvac) Escherichia coli cvi cvaC operon. [X57525 434-745]
kkonganti@11 703 # >SPG000143 (iss2) Escherichia coli Iss (iss) gene, complete cds. [AF042279 292-600]
kkonganti@11 704
kkonganti@11 705 for my $s (@$seqs) {
kkonganti@11 706
kkonganti@11 707 #print STDERR Dumper("IN", $s);
kkonganti@11 708 $s->{ID} =~ m/ ^ (\w+) (?: \( (.*?) \) )? $ /x
kkonganti@11 709 or die "Can't parse $fasta at " . Dumper($s);
kkonganti@11 710 $s->{ID} = $1 if $1;
kkonganti@11 711 $s->{ACC} = $2 || $1;
kkonganti@11 712 $s->{DESC} =~ s/\s\[.*?\]$//g; # remove strain name
kkonganti@11 713 $s->{DESC} =~ m/^(?:\((.*?)\)\s+)?(.*?)$/;
kkonganti@11 714 $s->{ID} = $1 if $1;
kkonganti@11 715 $s->{DESC} = $2;
kkonganti@11 716
kkonganti@11 717 #print STDERR Dumper("OUT", $s);
kkonganti@11 718 #print STDERR "="x60, "\n";
kkonganti@11 719 }
kkonganti@11 720
kkonganti@11 721 # print Dumper($seqs);
kkonganti@11 722 return $seqs;
kkonganti@11 723 }
kkonganti@11 724
kkonganti@11 725 #..............................................................................
kkonganti@11 726 sub is_full_gene {
kkonganti@11 727 my ($s) = @_;
kkonganti@11 728 my $has_ambig = 0;
kkonganti@11 729
kkonganti@11 730 my $id = $s->{ID};
kkonganti@11 731 my $L = length( $s->{SEQ} );
kkonganti@11 732 if ( $L % 3 != 0 ) {
kkonganti@11 733 wrn("$id - length $L bp is not multiple of 3");
kkonganti@11 734 return;
kkonganti@11 735 }
kkonganti@11 736 if ( $s->{SEQ} !~ m/^[AGCT]+$/ ) {
kkonganti@11 737 wrn("$id - has non-AGTC bases");
kkonganti@11 738 return;
kkonganti@11 739 }
kkonganti@11 740
kkonganti@11 741 my $seq = Bio::Seq->new( -id => $s->{ID}, -seq => $s->{SEQ} );
kkonganti@11 742
kkonganti@11 743 my $aa = $seq->translate->seq;
kkonganti@11 744
kkonganti@11 745 if ( $aa =~ m/\*./ ) {
kkonganti@11 746 wrn("$id - has internal stop codons, trying revcom");
kkonganti@11 747 $aa = $seq->revcom->translate->seq;
kkonganti@11 748 if ( $aa =~ m/\*./ ) {
kkonganti@11 749 wrn("$id - revcom has internal stop codons too");
kkonganti@11 750 return;
kkonganti@11 751 }
kkonganti@11 752 else {
kkonganti@11 753 msg("$id - revcom resolves problem, hooray!");
kkonganti@11 754 $s->{SEQ} = $seq->revcom->seq;
kkonganti@11 755 }
kkonganti@11 756 }
kkonganti@11 757
kkonganti@11 758 return $L;
kkonganti@11 759 }
kkonganti@11 760
kkonganti@11 761 #..............................................................................
kkonganti@11 762 sub dedupe_seq {
kkonganti@11 763 my ($seq) = @_;
kkonganti@11 764 my %seen;
kkonganti@11 765 my $good = [];
kkonganti@11 766 for my $s (@$seq) {
kkonganti@11 767 if ( $seen{ $s->{SEQ} } ) {
kkonganti@11 768 wrn( "duplicate", length( $s->{SEQ} ),
kkonganti@11 769 "bp sequence:", $s->{ID}, '~', $seen{ $s->{SEQ} } );
kkonganti@11 770 }
kkonganti@11 771 else {
kkonganti@11 772 push @$good, $s;
kkonganti@11 773 }
kkonganti@11 774 $seen{ $s->{SEQ} } .= ' ' . $s->{ID};
kkonganti@11 775 }
kkonganti@11 776 msg( "dedupe_seq: read", scalar(@$seq), "/ kept", scalar(@$good) );
kkonganti@11 777 return $good;
kkonganti@11 778 }
kkonganti@11 779
kkonganti@11 780 #..............................................................................
kkonganti@11 781 sub load_tabular {
kkonganti@11 782 my ( $fname, $keycol, $sep ) = @_;
kkonganti@11 783 $keycol //= 0;
kkonganti@11 784 $sep //= "\t";
kkonganti@11 785 my $hash = {};
kkonganti@11 786 my @hdr;
kkonganti@11 787 my $row = 0;
kkonganti@11 788 open my $TSV, '<', $fname or err("Can't read TSV file: $fname");
kkonganti@11 789 while (<$TSV>) {
kkonganti@11 790 chomp;
kkonganti@11 791 my @col = split m/$sep/;
kkonganti@11 792 $row++;
kkonganti@11 793 if (@hdr) {
kkonganti@11 794 @hdr == @col or err("Header and row $row differ in column count");
kkonganti@11 795
kkonganti@11 796 #my $key = $col[$keycol] or wrn("Empty key column $keycol: $_");
kkonganti@11 797 #exists{$hash->{$col[$key}} and wrn("WARNING: dupe key $key at row: $_");
kkonganti@11 798 my $key = $col[$keycol];
kkonganti@11 799 $hash->{$key} ||=
kkonganti@11 800 { map { ( $hdr[$_] => $col[$_] ) } ( 0 .. $#hdr ) };
kkonganti@11 801 }
kkonganti@11 802 else {
kkonganti@11 803 @hdr = @col;
kkonganti@11 804 }
kkonganti@11 805 }
kkonganti@11 806 close $TSV;
kkonganti@11 807 return $hash;
kkonganti@11 808 }
kkonganti@11 809
kkonganti@11 810 #..............................................................................
kkonganti@11 811 sub load_fasta {
kkonganti@11 812 my ($fasta) = @_;
kkonganti@11 813 my %seen;
kkonganti@11 814 my $list;
kkonganti@11 815 my $dbtype = 'unknown';
kkonganti@11 816 msg("load_fasta: $fasta");
kkonganti@11 817 my $in = Bio::SeqIO->new( -file => $fasta, -format => 'fasta' );
kkonganti@11 818 while ( my $seq = $in->next_seq ) {
kkonganti@11 819 my $id = $seq->id or err("Empty ID in $fasta");
kkonganti@11 820 if ( $seen{$id} ) {
kkonganti@11 821 wrn("Duplicate ID '$id' in $fasta");
kkonganti@11 822 $id = $id . '_dupe';
kkonganti@11 823 }
kkonganti@11 824 $seen{$id}++;
kkonganti@11 825 my $s = uc( $seq->seq );
kkonganti@11 826 $dbtype = $seq->alphabet eq 'dna' ? 'nucl' : 'prot';
kkonganti@11 827 $dbtype eq 'nucl' ? $s =~ s/[^AGTC]/N/g : $s =~ s/[^A-Z]/X/g;
kkonganti@11 828 push @$list,
kkonganti@11 829 {
kkonganti@11 830 ID => $id,
kkonganti@11 831 ACC => '',
kkonganti@11 832 DESC => $seq->desc,
kkonganti@11 833 SEQ => $s,
kkonganti@11 834 TYPE => $dbtype,
kkonganti@11 835 };
kkonganti@11 836 }
kkonganti@11 837 msg( "load_fasta: read", scalar(@$list), "$dbtype sequences" );
kkonganti@11 838 return $list;
kkonganti@11 839 }
kkonganti@11 840
kkonganti@11 841 #..............................................................................
kkonganti@11 842 sub save_fasta {
kkonganti@11 843 my ( $fasta, $seq ) = @_;
kkonganti@11 844 msg("save_fasta: $fasta");
kkonganti@11 845 my %seen;
kkonganti@11 846 my $out = Bio::SeqIO->new( -file => ">$fasta", -format => 'fasta' );
kkonganti@11 847 for my $s (@$seq) {
kkonganti@11 848 $seen{ $s->{ID} }++;
kkonganti@11 849 my $freq = $seen{ $s->{ID} };
kkonganti@11 850
kkonganti@11 851 #wrn("seen $s->{ID} now $freq times") if $freq > 1;
kkonganti@11 852 # print Dumper($s);
kkonganti@11 853 my $ABX =
kkonganti@11 854 defined( $s->{ABX} ) ? join( $ABX_SEP, sort @{ $s->{ABX} } ) : '';
kkonganti@11 855 $ABX =~ s/\s+/_/g; # remove spaces!
kkonganti@11 856 my $obj = Bio::Seq->new(
kkonganti@11 857 -id => join( '~~~', $db, $s->{ID}, $s->{ACC}, $ABX ),
kkonganti@11 858 -desc => ( $s->{DESC} || $s->{ID} ),
kkonganti@11 859 -seq => $s->{SEQ},
kkonganti@11 860 );
kkonganti@11 861
kkonganti@11 862 # $obj->desc( hash_encode($s) );
kkonganti@11 863 $out->write_seq($obj);
kkonganti@11 864
kkonganti@11 865 # $seen{ $s->{ID} }++;
kkonganti@11 866 }
kkonganti@11 867 msg( "save_fasta: wrote", scalar(@$seq), "sequences" );
kkonganti@11 868 }
kkonganti@11 869
kkonganti@11 870 #----------------------------------------------------------------------
kkonganti@11 871 sub msg {
kkonganti@11 872 print STDERR "@_\n";
kkonganti@11 873 }
kkonganti@11 874
kkonganti@11 875 #----------------------------------------------------------------------
kkonganti@11 876 sub wrn {
kkonganti@11 877 msg( "WARNING:", @_ ) if $debug;
kkonganti@11 878 }
kkonganti@11 879
kkonganti@11 880 #----------------------------------------------------------------------
kkonganti@11 881 sub err {
kkonganti@11 882 msg( "ERROR:", @_ );
kkonganti@11 883 exit(1);
kkonganti@11 884 }
kkonganti@11 885
kkonganti@11 886 #----------------------------------------------------------------------
kkonganti@11 887 # Option setting routines
kkonganti@11 888
kkonganti@11 889 sub setOptions {
kkonganti@11 890 use Getopt::Long;
kkonganti@11 891
kkonganti@11 892 @Options = (
kkonganti@11 893 { OPT => "help", VAR => \&usage, DESC => "This help" },
kkonganti@11 894 {
kkonganti@11 895 OPT => "debug!",
kkonganti@11 896 VAR => \$debug,
kkonganti@11 897 DEFAULT => 0,
kkonganti@11 898 DESC => "Verbose debug output"
kkonganti@11 899 },
kkonganti@11 900 {
kkonganti@11 901 OPT => "dbdir=s",
kkonganti@11 902 VAR => \$outdir,
kkonganti@11 903 DEFAULT => abs_path("$FindBin::RealBin/../db"),
kkonganti@11 904 DESC => "Parent folder"
kkonganti@11 905 },
kkonganti@11 906 {
kkonganti@11 907 OPT => "db=s",
kkonganti@11 908 VAR => \$db,
kkonganti@11 909 DEFAULT => "",
kkonganti@11 910 DESC => "Choices: $DATABASES"
kkonganti@11 911 },
kkonganti@11 912 {
kkonganti@11 913 OPT => "force!",
kkonganti@11 914 VAR => \$force,
kkonganti@11 915 DEFAULT => 0,
kkonganti@11 916 DESC => "Force download even if exists"
kkonganti@11 917 },
kkonganti@11 918 );
kkonganti@11 919
kkonganti@11 920 &GetOptions( map { $_->{OPT}, $_->{VAR} } @Options ) || usage(1);
kkonganti@11 921
kkonganti@11 922 # Now setup default values.
kkonganti@11 923 foreach (@Options) {
kkonganti@11 924 if ( defined( $_->{DEFAULT} ) && !defined( ${ $_->{VAR} } ) ) {
kkonganti@11 925 ${ $_->{VAR} } = $_->{DEFAULT};
kkonganti@11 926 }
kkonganti@11 927 }
kkonganti@11 928 }
kkonganti@11 929
kkonganti@11 930 sub usage {
kkonganti@11 931 my ($exitcode) = @_;
kkonganti@11 932 $exitcode = 0 if $exitcode eq 'help'; # what gets passed by getopt func ref
kkonganti@11 933 $exitcode ||= 0;
kkonganti@11 934 select STDERR if $exitcode; # write to STDERR if exitcode is error
kkonganti@11 935
kkonganti@11 936 print "SYNOPIS\n Download databases for abricate to use\n";
kkonganti@11 937 print "USAGE\n $EXE [options] --db DATABASE\n";
kkonganti@11 938 print "OPTIONS\n";
kkonganti@11 939 foreach (@Options) {
kkonganti@11 940 printf " --%-13s %s%s.\n", $_->{OPT}, $_->{DESC},
kkonganti@11 941 defined( $_->{DEFAULT} ) ? " (default '$_->{DEFAULT}')" : "";
kkonganti@11 942 }
kkonganti@11 943 exit($exitcode);
kkonganti@11 944 }
kkonganti@11 945
kkonganti@11 946 #----------------------------------------------------------------------
kkonganti@11 947