const void* rowWorse)

{

ASSERT (subtype_col != no_index);

ASSERT (start_col != no_index);

ASSERT (stop_col != no_index);

ASSERT (strand_col != no_index);

ASSERT (rowBetter);

ASSERT (rowWorse);

ASSERT (rowBetter != rowWorse);

const StringVector& rowBetter_ = * static_cast <const StringVector*> (rowBetter);

const StringVector& rowWorse_ = * static_cast <const StringVector*> (rowWorse);

ASSERT (rowBetter_ [strand_col] == rowWorse_ [strand_col]);

if ( rowBetter_ [subtype_col] == "POINT"

&& rowWorse_ [subtype_col] == "POINT_DISRUPT"

&& ( (rowBetter_ [strand_col] == "+" && rowBetter_ [start_col] == rowWorse_ [ start_col])

|| (rowBetter_ [strand_col] == "-" && rowBetter_ [stop_col] == rowWorse_ [ stop_col])

)

)

return 1;

return 0;

{

ThisApplication ()

: ShellApplication ("Identify AMR and virulence genes in proteins and/or contigs and print a report", true, true, true, true)

{

addFlag ("update", "Update the AMRFinder database", 'u'); // PD-2379

addFlag ("force_update", "Force updating the AMRFinder database", 'U'); // PD-3469

#ifdef DIR

addKey ("dir", "Common directory of the --protein, --nucleotide and --gff files", "", '\0', "DIRECTORY");

#endif

addKey ("protein", "Input protein FASTA file (can be gzipped)", "", 'p', "PROT_FASTA");

addKey ("nucleotide", "Input nucleotide FASTA file (can be gzipped)", "", 'n', "NUC_FASTA");

addKey ("gff", "GFF file for protein locations (can be gzipped). Locations are in the --nucleotide file. Protein ids should be in the attribute 'Name=<id>' (9th field) of the rows with type 'CDS' or 'gene' (3rd field).", "", 'g', "GFF_FILE");

{

const string annots (Gff::names. toString (", "));

addKey ("annotation_format", "Type of GFF file: " + annots, "genbank", 'a', "ANNOTATION_FORMAT");

}

addKey ("database", "Alternative directory with AMRFinder database. Default: $AMRFINDER_DB", "", 'd', "DATABASE_DIR");

addFlag ("database_version", "Print database version and exit", 'V');

addKey ("ident_min", "Minimum proportion of identical amino acids in alignment for hit (0..1). -1 means use a curated threshold if it exists and " + toString (ident_min_def) + " otherwise", "-1", 'i', "MIN_IDENT");

// "PD-3482

addKey ("coverage_min", "Minimum coverage of the reference protein (0..1)", toString (partial_coverage_min_def), 'c', "MIN_COV");

addKey ("organism", "Taxonomy group. To see all possible taxonomy groups use the --list_organisms flag", "", 'O', "ORGANISM");

addFlag ("list_organisms", "Print the list of all possible taxonomy groups for mutations identification and exit", 'l');

addKey ("translation_table", "NCBI genetic code for translated BLAST", "11", 't', "TRANSLATION_TABLE");

addFlag ("plus", "Add the plus genes to the report"); // PD-2789

//addFlag ("pseudo", "Report pseudogenes (frame shifts, internal stop codons)");

addFlag ("report_common", "Report proteins common to a taxonomy group"); // PD-2756

addFlag ("report_all_equal", "Report all equally-scoring BLAST and HMM matches"); // PD-3772

addKey ("name", "Text to be added as the first column \"name\" to all rows of the report, for example it can be an assembly name", "", '\0', "NAME");

addFlag ("print_node", "Print hierarchy node (family)"); // PD-4394

addKey ("mutation_all", "File to report all mutations", "", '\0', "MUT_ALL_FILE");

addKey ("output", "Write output to OUTPUT_FILE instead of STDOUT", "", 'o', "OUTPUT_FILE");

addKey ("protein_output", "Output protein FASTA file of reported proteins", "", '\0', "PROT_FASTA_OUT");

addKey ("nucleotide_output", "Output nucleotide FASTA file of reported nucleotide sequences", "", '\0', "NUC_FASTA_OUT");

addKey ("nucleotide_flank5_output", "Output nucleotide FASTA file of reported nucleotide sequences with 5' flanking sequences", "", '\0', "NUC_FLANK5_FASTA_OUT");

addKey ("nucleotide_flank5_size", "5' flanking sequence size for NUC_FLANK5_FASTA_OUT", "0", '\0', "NUC_FLANK5_SIZE");

addKey ("blast_bin", "Directory for BLAST. Deafult: $BLAST_BIN", "", '\0', "BLAST_DIR");

addKey ("hmmer_bin", "Directory for HMMer", "", '\0', "HMMER_DIR");

addFlag ("pgap", "Input files PROT_FASTA, NUC_FASTA and GFF_FILE are created by the NCBI PGAP"); // = --annotation_format pgap

addFlag ("gpipe_org", "NCBI internal GPipe organism names");

addKey ("parm", "amr_report parameters for testing: -nosame -noblast -skip_hmm_check -bed", "", '\0', "PARM");

version = SVN_REV;

documentationUrl = "https://github.com/ncbi/amr/wiki";

updatesUrl = "https://www.ncbi.nlm.nih.gov/mailman/listinfo/amrfinder-announce";

updatesDoc = "\

Subscribe to the amrfinder-announce mailing list for database and software update notifications";

}

void initEnvironment () final

{

ShellApplication::initEnvironment ();

var_cast (name2arg ["threads"] -> asKey ()) -> defaultValue = to_string (threads_def);

}

void fastaCheck (const string &fName,

bool prot,

const string &qcS,

const string &logFName,

size_t &nSeq,

size_t &len_max,

size_t &len_total,

const string &outFName) const

// Input: fName, outFName: quoted

{

ASSERT (fName != logFName);

if (! outFName. empty ())

{

ASSERT (outFName != fName);

ASSERT (outFName != logFName);

}

exec (fullProg ("fasta_check") + fName + " " + (prot ? "-aa -stop_codon -ambig_max " + ambigS + prependS (outFName, " -out ") : "-len " + tmp + "/len -hyphen -ambig") + qcS + " -log " + logFName + " > " + tmp + "/nseq", logFName);

// "-stop_codon" PD-4771

const StringVector vec (tmp + "/nseq", (size_t) 10, true);

if (vec. size () != 3)

throw runtime_error (string (prot ? "Protein" : "DNA") + " fasta_check failed: " + vec. toString ("\n"));

nSeq = str2<size_t> (vec [0]);

len_max = str2<size_t> (vec [1]);

len_total = str2<size_t> (vec [2]);

QC_ASSERT (nSeq);

QC_ASSERT (len_max);

QC_ASSERT (len_total);

}

StringVector db2organisms () const

{

const TextTable taxgroup (tmp + "/db/taxgroup.tsv");

const TextTable AMRProt_mutation (tmp + "/db/AMRProt-mutation.tsv");

const TextTable AMRProt_susceptible (tmp + "/db/AMRProt-susceptible.tsv");

taxgroup. qc ();

AMRProt_mutation. qc ();

AMRProt_susceptible. qc ();

StringVector vec (taxgroup. col2values (0));

vec << AMRProt_mutation. col2values (0);

vec << AMRProt_susceptible. col2values (0);

vec. sort ();

vec. uniq ();

return vec;

}

void prepare_fasta_extract (const StringVector &columns,

const string &tmpSuf,

bool saveHeader) const

// Input: tmp + "/amr"

{

TextTable t (tmp + "/amr");

t. qc ();

t. filterColumns (columns);

t. rows. filterValue ([] (const StringVector& row) { return row [0] == na; });

t. rows. sort ();

t. rows. uniq ();

t. saveHeader = saveHeader;

t. qc ();

t. saveFile (tmp + "/" + tmpSuf);

}

void amrTab_disruptions (TextTable &amrTab,

const string &db,

const string &dna_flat,

uint gencode,

const string &qcS) const

// PD-5301

{

const TextTable::ColNum contig_col = amrTab. col2num (contig_colName);

const TextTable::ColNum prot_col = amrTab. col2num (closestRefAccession_colName);

const TextTable::ColNum genesymbol_col = amrTab. col2num (genesymbol_colName);

TextTable disrRawTab;

disrRawTab. header << TextTable::Header ("contig")

<< TextTable::Header ("prot")

<< TextTable::Header ("disruption"); // Contains Disruption::genesymbol_raw()

disrRawTab. saveHeader = false;

for (const StringVector& row : amrTab. rows)

{

const string& contig = row [contig_col];

const string& prot = row [prot_col];

const string& genesymbol = row [genesymbol_col];

const size_t pos = genesymbol. find (disruption_delim);

if ( contig != na

&& prot != na

&& pos != string::npos

)

{

StringVector disrRow; disrRow. reserve (3);

disrRow << contig << prot << genesymbol. substr (pos + string (disruption_delim). size ());

disrRawTab. rows << std::move (disrRow);

}

}

disrRawTab. qc ();

if (disrRawTab. rows. empty ())

return;

disrRawTab. saveFile (tmp + "/disr_raw");

{

OFStream f (tmp + "/disr");

f << "#contig\tprot\tdisr\tdisr_raw\n";

// 0 1 2 3

}

exec (fullProg ("disruption2genesymbol") + dna_flat + " " + shellQuote (db + "/AMRProt-susceptible.fa")

+ " " + tmp + "/disr_raw -prot_id_pos 1 -gencode " + to_string (gencode) + qcS + " -noprogress >> " + tmp + "/disr");

const TextTable disrTab (tmp + "/disr");

disrTab. qc ();

for (StringVector& row : amrTab. rows)

{

const string& contig = row [contig_col];

const string& prot = row [prot_col];

string& genesymbol = row [genesymbol_col];

const size_t pos = genesymbol. find (disruption_delim);

if ( contig != na

&& prot != na

&& pos != string::npos

)

{

bool found = false;

const string disrS (genesymbol. substr (pos + string (disruption_delim). size ()));

for (const StringVector& disrRow : disrTab. rows)

if ( disrRow [0] == contig

&& disrRow [1] == prot

&& disrRow [3] == disrS

)

{

genesymbol. erase (pos);

genesymbol += "_" + disrRow [2];

found = true;

break;

}

if (! found)

throw runtime_error ("Disruption is not replaced by gene symbol:\n" + contig + " " + prot + " " + disrS);

}

}

}

#if 0

void execSpeed (const string &cmd,

const string &speedParam,

const string &coutFile,

const string &cerrFile,

const string &logFile) const

{

try { exec (cmd + " " + speedParam + " > " + coutFile + " 2> " + cerrFile, logFile); }

catch (...)

{ exec (cmd + " > " + coutFile + " 2> " + cerrFile, logFile); }

}

#endif

void shellBody () const final

{

const bool force_update = getFlag ("force_update");

const bool update = getFlag ("update") || force_update;

const string dir =

#ifdef DIR

appendS (getArg ("dir"), "/");

#else

"";

#endif

const string prot = shellQuote (prependS (getArg ("protein"), dir));

const string dna = shellQuote (prependS (getArg ("nucleotide"), dir));

const string gff = shellQuote (prependS (getArg ("gff"), dir));

string db = getArg ("database");

const bool pgap = getFlag ("pgap");

Gff::Type gffType = Gff::name2type (getArg ("annotation_format"));

const double ident = arg2double ("ident_min");

const double cov = arg2double ("coverage_min");

const string organism = shellQuote (getArg ("organism"));

const bool list_organisms = getFlag ("list_organisms");

const uint gencode = arg2uint ("translation_table");

const bool add_plus = getFlag ("plus");

const bool report_common = getFlag ("report_common");

const string mutation_all = getArg ("mutation_all");

string blast_bin = getArg ("blast_bin");

string hmmer_bin = getArg ("hmmer_bin");

const bool equidistant = getFlag ("report_all_equal");

const bool print_node = getFlag ("print_node");

const string input_name = shellQuote (getArg ("name"));

const string parm = getArg ("parm");

const string output = getArg ("output");

const string prot_out = shellQuote (getArg ("protein_output"));

const string dna_out = shellQuote (getArg ("nucleotide_output"));

const string dnaFlank5_out = shellQuote (getArg ("nucleotide_flank5_output"));

const uint dnaFlank5_size = arg2uint ("nucleotide_flank5_size");

const bool gpipe_org = getFlag ("gpipe_org");

const bool database_version = getFlag ("database_version");

const string logFName (tmp + "/log"); // Command-local log file

const bool screen = ! isRedirected (cerr);

if (database_version)

cout << "Software directory: " << shellQuote (execDir) << endl;

else

stderr << "Software directory: " << colorizeDir (execDir, screen) << '\n';

if (database_version)

cout << "Software version: " << version << endl;

else

stderr << "Software version: " << version << '\n';

if (contains (input_name, '\t'))

throw runtime_error ("NAME cannot contain a tab character");

#if 0

if (threads_max < threads_max_min)

throw runtime_error ("Number of threads cannot be less than " + to_string (threads_max_min));

#endif

if (ident != -1.0 && (ident < 0.0 || ident > 1.0))

throw runtime_error ("ident_min must be between 0 and 1");

if (cov < 0.0 || cov > 1.0)

throw runtime_error ("coverage_min must be between 0 and 1");

if (report_common && emptyArg (organism))

throw runtime_error ("--report_common requires --organism");

if (report_common && ! add_plus)

throw runtime_error ("--report_common requires --plus");

// PD-3437

if (! mutation_all. empty () && emptyArg (organism))

{

const Warning warning (stderr);

stderr << "--mutation_all option used without -O/--organism option. No point mutations will be screened";

}

OFStream::prepare (output);

string defaultDb;

#ifdef CONDA_DB_DIR

// we're in condaland

if (const char* s = getenv("CONDA_PREFIX")) {

defaultDb = string (s) + "/share/amrfinderplus/data/latest";

} else if (const char* s = getenv("PREFIX")) {

const Warning warning (stderr);

stderr << "This was compiled for running under bioconda, but $CONDA_PREFIX was not found" << '\n';

defaultDb = string (s) + "/share/amrfinderplus/data/latest";

stderr << "Reverting to $PREFIX: " << colorizeDir (defaultDb, screen);

} else {

const Warning warning (stderr);

stderr << "This was compiled for running under bioconda, but $CONDA_PREFIX was not found" << '\n';

stderr << "Reverting to hard coded directory: " << colorizeDir (CONDA_DB_DIR "/latest", screen);

defaultDb = CONDA_DB_DIR "/latest";

}

#else

// not in condaland

defaultDb = execDir + "data/latest";

#endif

ASSERT (isRight (defaultDb, "/latest"));

// db

if (db. empty ())

{

if (const char* s = getenv ("AMRFINDER_DB"))

db = s;

else

db = defaultDb;

}

ASSERT (! db. empty ());

// blast_bin

if (blast_bin. empty ())

if (const char* s = getenv ("BLAST_BIN"))

blast_bin = string (s);

if (! blast_bin. empty ())

{

addDirSlash (blast_bin);

prog2dir ["blastp"] = blast_bin;

prog2dir ["blastx"] = blast_bin;

prog2dir ["tblastn"] = blast_bin;

prog2dir ["blastn"] = blast_bin;

//prog2dir ["makeblastdb"] = blast_bin;

}

if (! hmmer_bin. empty ())

{

addDirSlash (hmmer_bin);

prog2dir ["hmmsearch"] = hmmer_bin;

}

if (update)

{

// PD-2447

if (! emptyArg (prot) || ! emptyArg (dna))

throw runtime_error ("AMRFinder -u/--update option cannot be run with -n/--nucleotide or -p/--protein options");

if (! getArg ("database"). empty ())

throw runtime_error ("AMRFinder update option (-u/--update) only operates on the default database directory. The -d/--database option is not permitted");

if (getenv ("AMRFINDER_DB"))

{

const Warning warning (stderr);

stderr << "AMRFINDER_DB is set, but AMRFinder auto-update only downloads to the default database directory";

db = defaultDb;

}

const Dir dbDir (db);

if (! dbDir. items. empty () && dbDir. items. back () == "latest")

{

prog2dir ["amrfinder_update"] = execDir;

exec (fullProg ("amrfinder_update") + " -d " + shellQuote (dbDir. getParent ()) + ifS (force_update, " --force_update")

+ makeKey ("blast_bin", blast_bin)

+ makeKey ("hmmer_bin", hmmer_bin)

+ ifS (getQuiet (), " -q") + ifS (qc_on, " --debug") + " > " + logFName, logFName);

}

else

{

const Warning warning (stderr);

stderr << "Updating database directory works only for databases with the default data directory format." << '\n'

<< " Please see " + colorizeUrl ("https://github.com/ncbi/amr/wiki", screen) + " for details." << "\n"

<< " Current database directory is: " << colorizeDir (dbDir. get (), screen) << "\n"

<< " New database directories will be created as subdirectories of: " << colorizeDir (dbDir. getParent (), screen);

}

}

const string downloadLatestInstr ("\nTo download the latest version to the default directory run: amrfinder -u");

if (! directoryExists (db)) // PD-2447

//throw runtime_error ("No valid AMRFinder database found: " + db + ifS (! update, downloadLatestInstr));

throw runtime_error ("No valid AMRFinder database is found.\nThis directory (or symbolic link to directory) is not found: " + db + ifS (! update, downloadLatestInstr));

if (database_version)

cout << "Database directory: " << shellQuote (path2canonical (db)) << endl;

else

stderr << "Database directory: " << colorizeDir (path2canonical (db), screen) << '\n';

setSymlink (db, tmp + "/db", true);

{

// PD-4925

const string dbTest (db + "/AMRProt.fa.phr");

if (! fileExists (dbTest))

throw runtime_error ("The BLAST database for AMRProt.fa (" + dbTest + ") was not found.\nUse amrfinder -u or amrfinder --force_update to download and prepare database for AMRFinderPlus");

}

// PD-3051

{

istringstream versionIss (version);

const SoftwareVersion softwareVersion (versionIss);

const SoftwareVersion softwareVersion_min (db + "/database_format_version.txt");

//stderr << "Software version: " << softwareVersion. str () << '\n';

const DataVersion dataVersion (db + "/version.txt");

istringstream dataVersionIss (dataVer_min);

const DataVersion dataVersion_min (dataVersionIss);

if (database_version)

cout << "Database version: " << dataVersion. str () << endl;

else

stderr << "Database version: " << dataVersion. str () << '\n';

if (softwareVersion < softwareVersion_min)

throw runtime_error ("Database requires software version at least " + softwareVersion_min. str ());

if (dataVersion < dataVersion_min)

throw runtime_error ("Software requires database version at least " + dataVersion_min. str () + downloadLatestInstr);

if (database_version)

{

if (! emptyArg (prot))

throw runtime_error ("No processing of " + prot + " is done");

if (! emptyArg (dna))

throw runtime_error ("No processing of " + dna + " is done");

return;

}

}

if (list_organisms)

{

const StringVector organisms (db2organisms ());

cout << endl << "Available --organism options: " + organisms. toString (", ") << endl;

return;

}

{

string searchMode;

StringVector includes;

if (emptyArg (prot))

{

if (emptyArg (dna))

{

if (update)

return;

throw runtime_error ("Parameter --protein or --nucleotide must be present");

}

else

{

if (! emptyArg (gff))

throw runtime_error ("Parameter --gff is redundant");

searchMode = "translated nucleotide";

}

}

else

{

searchMode = "protein";

if (emptyArg (dna))

{

searchMode += "-only";

includes << key2shortHelp ("nucleotide") + " and " + key2shortHelp ("gff") + " options to add translated searches";

}

else

{

if (emptyArg (gff))

throw runtime_error ("If parameters --protein and --nucleotide are present then parameter --gff must be present");

searchMode = "combined translated and protein";

}

}

if (emptyArg (prot) && ! emptyArg (prot_out))

throw runtime_error ("Parameter --protein must be present for --protein_output");

if (emptyArg (dna) && ! emptyArg (dna_out))

throw runtime_error ("Parameter --nucleotide must be present for --nucleotide_output");

if (emptyArg (dna) && ! emptyArg (dnaFlank5_out))

throw runtime_error ("Parameter --nucleotide must be present for --nucleotide_flank5_output");

if (emptyArg (dnaFlank5_out) && dnaFlank5_size > 0)

throw runtime_error ("Parameter --nucleotide_flank5_output must be present for --nucleotide_flank5_size");

if (! emptyArg (dnaFlank5_out) && dnaFlank5_size == 0)

throw runtime_error ("Parameter --nucleotide_flank5_size must be present with a positive value for --nucleotide_flank5_output");

ASSERT (! searchMode. empty ());

if (emptyArg (organism))

includes << key2shortHelp ("organism") + " option to add mutation searches and suppress common proteins";

else

searchMode += " and mutation";

stderr << "AMRFinder " << searchMode << " search\n";

for (const string& include : includes)

stderr << colorize (" - include " + include, screen) << '\n';

}

// Quoted names

const string prot_flat = uncompress (prot, "prot_flat");

const string dna_flat = uncompress (dna, "dna_flat");

const string gff_flat = uncompress (gff, "gff_flat");

{

StringVector emptyFiles;

if (! emptyArg (prot) && ! getFileSize (unQuote (prot_flat))) emptyFiles << prot;

if (! emptyArg (dna) && ! getFileSize (unQuote (dna_flat))) emptyFiles << dna;

if (! emptyArg (gff) && ! getFileSize (unQuote (gff_flat))) emptyFiles << gff;

for (const string& emptyFile : emptyFiles)

{

const Warning warning (stderr);

stderr << "Empty file: " << emptyFile;

}

}

// organism --> organism1

string organism1;

bool suppress_common = false;

if (! emptyArg (organism))

{

organism1 = unQuote (organism);

replace (organism1, ' ', '_');

ASSERT (! organism1. empty ());

if (gpipe_org)

{

const TextTable tab (db + "/taxgroup.tsv");

tab. qc ();

const TextTable::ColNum taxgroup_col = tab. col2num ("taxgroup");

const TextTable::ColNum gpipeTaxgroup_col = tab. col2num ("gpipe_taxgroup");

bool found = false;

for (const StringVector& row : tab. rows)

{

QC_ASSERT (! row [taxgroup_col]. empty ());

const StringVector gpipeOrgVec (row [gpipeTaxgroup_col], ',', true);

//QC_ASSERT (gpipeOrgVec. size () >= 1);

if (gpipeOrgVec. contains (organism1))

{

organism1 = row [taxgroup_col];

found = true;

break;

}

}