corrected ANI using orthoANI_usearch/blastn

README.md

@@ -42,6 +42,10 @@ For above use case, REFERENCE\_LIST should be a file containing directory paths
 $ ./fastANI --ql [QUERY_LIST] --rl [REFERENCE_LIST] -o [OUTPUT_FILE]
 ```
 Again, QUERY\_LIST and REFERENCE\_LIST are files containing paths to genomes, one per line.
+* **Corrected ANI values.** the final ANI values can be corrected by learning from alignment-based ANI (e.g., orthoANI), expecially for ANI range 70% to 85%.
+```sh
+$ ./fastANI --ql [QUERY_LIST] --rl [REFERENCE_LIST] -o [OUTPUT_FILE] --correct
+```
 
 **Output format.** In all above use cases, OUTPUT\_FILE will contain tab delimited row(s) with query genome, reference genome, ANI value, count of bidirectional fragment mappings, and total query fragments. Alignment fraction (wrt. the query genome) is simply the ratio of mappings and total fragments. Optionally, users can also get a second `.matrix` file with identity values arranged in a [phylip-formatted lower triangular matrix](https://www.mothur.org/wiki/Phylip-formatted_distance_matrix) by supplying `--matrix` parameter. **NOTE:** No ANI output is reported for a genome pair if ANI value is much below 80%. Such case should be computed at [amino acid level](http://enve-omics.ce.gatech.edu/aai/).
 

src/cgi/include/computeCoreIdentity.hpp

@@ -289,7 +289,11 @@ namespace cgi
       currentResult.countSeq = std::distance(it, rangeEndIter);
       currentResult.totalQueryFragments = totalQueryFragments;
       currentResult.identity = sumIdentity/currentResult.countSeq;
-
+      if(parameters.correct)
+      {
+      // linear correction of final identity values according to large scale orthoANI(ANI_u) computations
+        currentResult.identity = currentResult.identity * 41.00/36.00 - 500.00/36.00;
+      }
       CGI_ResultsVector.push_back(currentResult);
 
       //Advance the iterator it

src/map/include/map_parameters.hpp

@@ -34,7 +34,8 @@ namespace skch
     std::vector<std::string> querySequences;          //query sequence(s)
     std::string outFileName;                          //output file name
     bool reportAll;                                   //Report all alignments if this is true
-    bool visualize;                                   //Visualize the conserved regions of two genomes
+    bool visualize;
+    bool correct;                                 //Visualize the conserved regions of two genomes
     bool matrixOutput;                                //report fastani results as lower triangular matrix
     float maxRatioDiff;                               //max Ratio for sanity check
     bool sanityCheck;                                 // Sanity check for extreme Cases

src/map/include/parseCmdArgs.hpp

@@ -123,6 +123,7 @@ namespace skch
     parameters.p_value = 1e-03;
     parameters.percentageIdentity = 80;
     parameters.visualize = false;
+    parameters.correct = false;
     parameters.matrixOutput = false;
     parameters.referenceSize = 5000000;
     parameters.maxRatioDiff = 100.0;
@@ -146,6 +147,7 @@ namespace skch
     auto minfraction_cmd = (clipp::option("--minFraction") & clipp::value("value", parameters.minFraction)) % "minimum fraction of genome that must be shared for trusting ANI. If reference and query genome size differ, smaller one among the two is considered. [default : 0.2]";
     auto maxratio_cmd = (clipp::option("--maxRatioDiff") & clipp::value("value", parameters.maxRatioDiff)) % "maximum difference between (Total Ref. Length/Total Occ. Hashes) and (Total Ref. Length/Total No. Hashes). [default : 10.0]";
     auto visualize_cmd = clipp::option("--visualize").set(parameters.visualize).doc("output mappings for visualization, can be enabled for single genome to single genome comparison only [disabled by default]");
+    auto correct_cmd = clipp::option("--correct").set(parameters.correct).doc("correct final ANI values by learning alignment-based ANI values");
     auto matrix_cmd = clipp::option("--matrix").set(parameters.matrixOutput).doc("also output ANI values as lower triangular matrix (format inspired from phylip). If enabled, you should expect an output file with .matrix extension [disabled by default]");
     auto output_cmd = (clipp::option("-o", "--output") & clipp::value("value", parameters.outFileName)) % "output file name";
     auto sanitycheck_cmd = clipp::option("-s", "--sanityCheck").set(parameters.sanityCheck).doc("run sanity check");
@@ -164,6 +166,7 @@ namespace skch
        minfraction_cmd,
        maxratio_cmd,
        visualize_cmd,
+       correct_cmd,
        matrix_cmd,
        output_cmd,
        sanitycheck_cmd,