shared genome threshold updated

Author: cjain7

src/cgi/core_genome_identity.cpp

@@ -115,10 +115,13 @@ int main(int argc, char** argv)
 
   std::cerr << "INFO, skch::main, parallel_for execution finished" << std::endl;
 
+  std::unordered_map <std::string, uint64_t> genomeLengths;    // name of genome -> length
+  cgi::computeGenomeLengths(parameters, genomeLengths);
+
   //report output in file
-  cgi::outputCGI (parameters, finalResults, fileName);
+  cgi::outputCGI (parameters, genomeLengths, finalResults, fileName);
 
   //report output as matrix
   if (parameters.matrixOutput)
-    cgi::outputPhylip (parameters, finalResults, fileName);
+    cgi::outputPhylip (parameters, genomeLengths, finalResults, fileName);
 }

src/cgi/include/computeCoreIdentity.hpp

@@ -11,12 +11,14 @@
 #include <unordered_map>
 #include <fstream>
 #include <omp.h>
+#include <zlib.h>  
 
 //Own includes
 #include "map/include/base_types.hpp"
 #include "cgi/include/cgid_types.hpp"
 
 //External includes
+#include "common/kseq.h"
 #include "common/prettyprint.hpp"
 
 namespace cgi
@@ -39,6 +41,52 @@ namespace cgi
     }
   }
 
+  /**
+   * @brief                       compute genome lengths in reference and query genome set
+   * @param[out] genomeLengths
+   */
+  void computeGenomeLengths(skch::Parameters &parameters, std::unordered_map <std::string, uint64_t> &genomeLengths) 
+  { 
+    for(auto &e : parameters.querySequences)
+    {
+      //Open the file using kseq
+      FILE *file = fopen(e.c_str(), "r");
+      gzFile fp = gzdopen(fileno(file), "r");
+      kseq_t *seq = kseq_init(fp);
+      int l; uint64_t genomeLen = 0;
+
+      while ((l = kseq_read(seq)) >= 0) 
+        genomeLen = genomeLen + (uint64_t)strlen(seq->seq.s);
+
+      genomeLengths[e] = genomeLen;
+
+      kseq_destroy(seq);  
+      gzclose(fp); //close the file handler 
+      fclose(file);
+    }
+
+    for(auto &e : parameters.refSequences)
+    {
+      if( genomeLengths.find(e) == genomeLengths.end() )
+      {
+        //Open the file using kseq
+        FILE *file = fopen(e.c_str(), "r");
+        gzFile fp = gzdopen(fileno(file), "r");
+        kseq_t *seq = kseq_init(fp);
+        int l; uint64_t genomeLen = 0;
+
+        while ((l = kseq_read(seq)) >= 0) 
+          genomeLen = genomeLen + (uint64_t)strlen(seq->seq.s);
+
+        genomeLengths[e] = genomeLen;
+
+        kseq_destroy(seq);  
+        gzclose(fp); //close the file handler 
+        fclose(file);
+      }
+    }
+  }
+
   /**
    * @brief                             output blast tabular mappings for visualization 
    * @param[in]   parameters            algorithm parameters
@@ -242,10 +290,12 @@ namespace cgi
   /**
    * @brief                             output FastANI results to file
    * @param[in]   parameters            algorithm parameters
+   * @param[in]   genomeLengths
    * @param[in]   CGI_ResultsVector     results
    * @param[in]   fileName              file name where results will be reported
    */
   void outputCGI(skch::Parameters &parameters,
+      std::unordered_map <std::string, uint64_t> &genomeLengths,
       std::vector<cgi::CGI_Results> &CGI_ResultsVector,
       std::string &fileName)
   {
@@ -257,10 +307,22 @@ namespace cgi
     //Report results
     for(auto &e : CGI_ResultsVector)
     {
-      if(e.countSeq >= parameters.minFragments)
+      std::string qryGenome = parameters.querySequences[e.qryGenomeId];
+      std::string refGenome = parameters.refSequences[e.refGenomeId];
+
+      assert(genomeLengths.find(qryGenome) != genomeLengths.end());
+      assert(genomeLengths.find(refGenome) != genomeLengths.end());
+
+      uint64_t queryGenomeLength = genomeLengths[qryGenome];
+      uint64_t refGenomeLength = genomeLengths[refGenome]; 
+      uint64_t minGenomeLength = std::min(queryGenomeLength, refGenomeLength);
+      uint64_t sharedLength = e.countSeq * parameters.minReadLength;
+
+      //Checking if shared genome is above a certain fraction of genome length
+      if(sharedLength >= minGenomeLength * parameters.minFraction)
       {
-        outstrm << parameters.querySequences[e.qryGenomeId]
-          << "\t" << parameters.refSequences[e.refGenomeId]
+        outstrm << qryGenome
+          << "\t" << refGenome
           << "\t" << e.identity 
           << "\t" << e.countSeq
           << "\t" << e.totalQueryFragments
@@ -274,10 +336,12 @@ namespace cgi
   /**
    * @brief                             output FastANI results as lower triangular matrix
    * @param[in]   parameters            algorithm parameters
+   * @param[in]   genomeLengths
    * @param[in]   CGI_ResultsVector     results
    * @param[in]   fileName              file name where results will be reported
    */
   void outputPhylip(skch::Parameters &parameters,
+      std::unordered_map <std::string, uint64_t> &genomeLengths,
       std::vector<cgi::CGI_Results> &CGI_ResultsVector,
       std::string &fileName)
   {
@@ -313,10 +377,22 @@ namespace cgi
     //transform FastANI results into 3-tuples
     for(auto &e : CGI_ResultsVector)
     {
-      if(e.countSeq >= parameters.minFragments)
+      std::string qryGenome = parameters.querySequences[e.qryGenomeId];
+      std::string refGenome = parameters.refSequences[e.refGenomeId];
+
+      assert(genomeLengths.find(qryGenome) != genomeLengths.end());
+      assert(genomeLengths.find(refGenome) != genomeLengths.end());
+
+      uint64_t queryGenomeLength = genomeLengths[qryGenome];
+      uint64_t refGenomeLength = genomeLengths[refGenome]; 
+      uint64_t minGenomeLength = std::min(queryGenomeLength, refGenomeLength);
+      uint64_t sharedLength = e.countSeq * parameters.minReadLength;
+
+      //Checking if shared genome is above a certain fraction of genome length
+      if(sharedLength >= minGenomeLength * parameters.minFraction)
       {
-        int qGenome = genome2Int [ parameters.querySequences[e.qryGenomeId] ];
-        int rGenome = genome2Int [ parameters.refSequences[e.refGenomeId] ];
+        int qGenome = genome2Int [ qryGenome ];
+        int rGenome = genome2Int [ refGenome ];
 
         if (qGenome != rGenome)   //ignore if both genomes are same
         {

src/map/include/map_parameters.hpp

@@ -23,7 +23,7 @@ namespace skch
     int kmerSize;                                     //kmer size for sketching
     int windowSize;                                   //window size used for sketching 
     int minReadLength;                                //minimum read length which code maps
-    int minFragments;                                 //minimum fragment mappings for trusting ANI value
+    float minFraction;                                //minimum genome fraction for trusting ANI value
     int threads;                                      //thread count
     int alphabetSize;                                 //alphabet size
     uint64_t referenceSize;                           //Approximate reference size

src/map/include/parseCmdArgs.hpp

@@ -10,6 +10,7 @@
 #include <iostream>
 #include <string>
 #include <fstream>
+#include <cassert>
 
 //Own includes
 #include "map/include/map_parameters.hpp"
@@ -57,7 +58,7 @@ Example usage: \n\
 
     cmd.defineOption("fragLen", "fragment length [default : 3,000]", ArgvParser::OptionRequiresValue);
 
-    cmd.defineOption("minFrag", "minimum matched fragments for trusting ANI [default : 50]", ArgvParser::OptionRequiresValue);
+    cmd.defineOption("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]", ArgvParser::OptionRequiresValue);
 
     cmd.defineOption("visualize", "output mappings for visualization, can be enabled for single genome to single genome comparison only [disabled by default]");
 
@@ -167,20 +168,24 @@ Example usage: \n\
     if (cmd.foundOption("version"))
     {
       std::cerr << "version 1.2\n\n";
-      exit(1);
+      exit(0);
     }
     if (result != ArgvParser::NoParserError)
     {
       std::cerr << cmd.parseErrorDescription(result) << "\n";
-      exit(1);
+
+      if (result == ArgvParser::ParserHelpRequested)
+        exit(0);
+      else
+        exit(1);
     }
     else if (!cmd.foundOption("ref") && !cmd.foundOption("refList"))
-    {
+    { 
       std::cerr << "Provide reference file (s)\n";
       exit(1);
     }
     else if (!cmd.foundOption("query") && !cmd.foundOption("queryList"))
-    {
+    { 
       std::cerr << "Provide query file (s)\n";
       exit(1);
     }
@@ -296,14 +301,15 @@ Example usage: \n\
     else
       parameters.minReadLength = 3000;
 
-    if(cmd.foundOption("minFrag"))
+    if(cmd.foundOption("minFraction"))
     {
-      str << cmd.optionValue("minFrag");
-      str >> parameters.minFragments;
+      str << cmd.optionValue("minFraction");
+      str >> parameters.minFraction;
       str.clear();
+      assert(parameters.minFraction >= 0.0 && parameters.minFraction <= 1.0);
     }
     else
-      parameters.minFragments = 50;
+      parameters.minFraction = 0.2;
 
     parameters.percentageIdentity = 80;