Added new argument -G to count gaps towards depth for variants.

Author: gkarthik

src/call_variants.cpp

@@ -32,7 +32,7 @@ double *get_frequency_depth(
 int call_variants_from_plup(std::istream &cin, std::string out_file,
                             uint8_t min_qual, double min_threshold,
                             uint32_t min_depth, std::string ref_path,
-                            std::string gff_path) {
+                            std::string gff_path, bool gapped_depth) {
   std::string line, cell, bases, qualities, region;
   ref_antd refantd(ref_path, gff_path);
   std::ostringstream out_str;
@@ -108,7 +108,13 @@ int call_variants_from_plup(std::istream &cin, std::string out_file,
     // Get ungapped depth
     pdepth = 0;
     for (std::vector<allele>::iterator it = ad.begin(); it != ad.end(); ++it) {
-      if (it->nuc[0] == '*' || it->nuc[0] == '+' || it->nuc[0] == '-') continue;
+      // if ((gapped_depth && (it->nuc[0] == '+' || it->nuc[0] == '-')) || (!gapped_depth && (it->nuc[0] == '+' || it->nuc[0] == '*' || it->nuc[0] == '-'))) continue;
+      // Split into two if statements for readability.
+      if (gapped_depth) {
+        if (it->nuc[0] == '+' || it->nuc[0] == '-') continue; // Count gaps ('*') towards depth
+      } else {
+        if (it->nuc[0] == '+' || it->nuc[0] == '*' || it->nuc[0] == '-') continue;
+      }
       pdepth += it->depth;
     }
     if (pdepth < min_depth) {  // Check for minimum depth

src/call_variants.h

@@ -19,7 +19,7 @@
 int call_variants_from_plup(std::istream &cin, std::string out_file,
                             uint8_t min_qual, double min_threshold,
                             uint32_t min_depth, std::string ref_path,
-                            std::string gff_path);
+                            std::string gff_path, bool gapped_depth);
 std::vector<allele>::iterator get_ref_allele(std::vector<allele> &ad, char ref);
 
 #endif

src/ivar.cpp

@@ -44,6 +44,7 @@ struct args_t {
   bool write_no_primers_flag;    // -e
   std::string gff;               // -g
   bool keep_for_reanalysis;      // -k
+  bool gapped_depth;             // -G
 } g_args;
 
 void print_usage() {
@@ -110,6 +111,7 @@ void print_variants_usage() {
          "           -t    Minimum frequency threshold(0 - 1) to call variants "
          "(Default: 0.03)\n"
          "           -m    Minimum read depth to call variants (Default: 0)\n"
+         "           -G    Count gaps towards depth. By default, gaps are not counted\n"
          "           -r    Reference file used for alignment. This is used to "
          "translate the nucleotide sequences and identify intra host single "
          "nucleotide variants\n"
@@ -241,7 +243,7 @@ void print_version_info() {
 }
 
 static const char *trim_opt_str = "i:b:f:x:p:m:q:s:ekh?";
-static const char *variants_opt_str = "p:t:q:m:r:g:h?";
+static const char *variants_opt_str = "p:t:q:m:r:g:Gh?";
 static const char *consensus_opt_str = "i:p:q:t:c:m:n:kh?";
 static const char *removereads_opt_str = "i:p:t:b:h?";
 static const char *filtervariants_opt_str = "p:t:f:h?";
@@ -361,6 +363,7 @@ int main(int argc, char *argv[]) {
     g_args.min_depth = 0;
     g_args.ref = "";
     g_args.gff = "";
+    g_args.gapped_depth = false;
     opt = getopt(argc, argv, variants_opt_str);
     while (opt != -1) {
       switch (opt) {
@@ -382,6 +385,9 @@ int main(int argc, char *argv[]) {
         case 'g':
           g_args.gff = optarg;
           break;
+        case 'G':
+          g_args.gapped_depth = true;
+          break;
         case 'h':
         case '?':
           print_variants_usage();
@@ -422,7 +428,7 @@ int main(int argc, char *argv[]) {
     }
     res = call_variants_from_plup(std::cin, g_args.prefix, g_args.min_qual,
                                   g_args.min_threshold, g_args.min_depth,
-                                  g_args.ref, g_args.gff);
+                                  g_args.ref, g_args.gff, g_args.gapped_depth);
   }
   // ivar consensus
   else if (cmd.compare("consensus") == 0) {

tests/test_variants.cpp

@@ -11,7 +11,7 @@ int call_var_check_outfile(std::string prefix, uint8_t min_qual,
   std::string path = "../data/test.indel.mpileup";
   std::ifstream mplp(path);
   call_variants_from_plup(mplp, prefix, min_qual, min_threshold, min_depth,
-                          "../data/db/test_ref.fa", "../data/test.gff");
+                          "../data/db/test_ref.fa", "../data/test.gff", false);
   std::ifstream outFile(prefix + ".tsv");
   std::string l;
   getline(outFile, l);  // Ignore first line

tests/test_variants_reverse_strand.cpp

@@ -11,7 +11,7 @@ int call_var_check_outfile(std::string prefix, uint8_t min_qual,
   std::string path = "../data/test.strand.pileup";
   std::ifstream mplp(path);
   call_variants_from_plup(mplp, prefix, min_qual, min_threshold, min_depth,
-                          "../data/db/test_ref.fa", "../data/test_strand.gff");
+                          "../data/db/test_ref.fa", "../data/test_strand.gff", false);
   std::ifstream outFile(prefix + ".tsv");
   std::string l;
   getline(outFile, l);  // Ignore first line