Update bustools inspect to account for multi-barcode onlists (e.g for SPLiT-seq)

src/bustools_inspect.cpp

@@ -50,22 +50,83 @@ void bustools_inspect(Bustools_opt &opt) {
     }
   }
 
-  /* Load whitelist. */
-  std::unordered_set<uint64_t> whitelist;
+  /* Load on-list. A "complex" on-list has one column per barcode subsequence, and any
+     combination of entries across the columns forms a valid barcode, so each column gets
+     its own set (as in bustools correct) rather than concatenating each line into one
+     barcode. Ragged columns are padded with "-". */
+  std::vector<std::unordered_set<uint64_t>> wbc; // Each set contains the on-list of one column
+  std::vector<uint32_t> wc_bclen; // Barcode length of each column
   if (opt.whitelist.size()) {
     std::ifstream wl(opt.whitelist);
-    std::string inp;
+    std::string line;
     uint32_t flag; // Unused
-    while (std::getline(wl, inp)) {
-      std::string str = inp;
-      str.erase(std::remove_if(str.begin(), str.end(), [](unsigned char c) {
-        return c == ' ' || c == '\t'; // Remove spaces and tabs (e.g. if we have split barcodes in our list)
-      }), str.end());
-      whitelist.insert(stringToBinary(str, flag));
+    bool first_line = true;
+    while (std::getline(wl, line)) {
+      std::stringstream ss(line);
+      std::string barcode;
+      size_t i = 0;
+      while (ss >> barcode) {
+        std::transform(barcode.begin(), barcode.end(), barcode.begin(), ::toupper);
+        if (barcode == "-") { // Empty barcode (padding for a ragged column)
+          if (first_line) { // The first line establishes each column's length, so it can't be padded
+            std::cerr << "Error: on-list file malformed; the first line cannot contain empty barcodes"
+                      << std::endl;
+            exit(1);
+          }
+          ++i;
+          continue;
+        }
+        uint64_t bc = stringToBinary(barcode, flag);
+        if (first_line) { // First line establishes all the barcode sets
+          std::unordered_set<uint64_t> bc_set;
+          bc_set.reserve(100000);
+          bc_set.insert(bc);
+          wbc.push_back(std::move(bc_set));
+          wc_bclen.push_back(barcode.size());
+        } else if (i >= wbc.size()) { // Too many barcodes in this line
+          std::cerr << "Error: on-list file malformed; encountered " << (i+1)
+                    << " barcodes on a line while " << wbc.size() << " barcodes on a previous line"
+                    << std::endl;
+          exit(1);
+        } else if (barcode.length() != wc_bclen[i]) {
+          std::cerr << "Error: on-list file malformed; encountered barcode length " << wc_bclen[i]
+                    << " on a line but barcode length " << barcode.length() << " on another line"
+                    << std::endl;
+          exit(1);
+        } else {
+          wbc[i].insert(bc);
+        }
+        ++i;
+      }
+      if (i == 0) continue; // Empty line
+      first_line = false;
     }
     wl.close();
+
+    if (wbc.size() == 0) {
+      std::cerr << "Error: on-list file malformed; no barcodes found" << std::endl;
+      exit(1);
+    }
   }
 
+  /* A barcode is on the on-list if each of its subsequences is found in the set of the
+     column it comes from; the last column occupies the least significant bits. */
+  auto in_onlist = [&wbc, &wc_bclen](uint64_t bc) -> bool {
+    if (wbc.empty()) {
+      return false;
+    }
+    uint64_t shift = 0;
+    for (size_t j = wbc.size(); j-- > 0; ) {
+      uint32_t nbits = 2 * wc_bclen[j];
+      uint64_t mask = nbits >= 64 ? ~0ULL : ((1ULL << nbits) - 1);
+      if (wbc[j].find((bc >> shift) & mask) == wbc[j].end()) {
+        return false;
+      }
+      shift += nbits;
+    }
+    return true;
+  };
+
   /* Inspect. */
   size_t N = 100000;
   BUSData *p = new BUSData[N];
@@ -83,6 +144,19 @@ void bustools_inspect(Bustools_opt &opt) {
   uint32_t bclen = h.bclen;
   uint64_t len_mask = ((1ULL << (2*bclen)) - 1); // Only include n least significant bits where n=2*bclen
 
+  if (wbc.size()) { // Warn (rather than exit) so the remaining stats are still reported
+    uint32_t wl_bclen = 0;
+    for (auto l : wc_bclen) {
+      wl_bclen += l;
+    }
+    if (bclen != wl_bclen) {
+      std::cerr << "Warning: barcode length and on-list length differ, barcodes = " << bclen
+                << ", on-list = " << wl_bclen << std::endl
+                << "         check that your on-list matches the technology used;"
+                << " on-list statistics will be unreliable" << std::endl;
+    }
+  }
+
   /* Number of records. */
   size_t nr = 0;
 
@@ -139,7 +213,7 @@ void bustools_inspect(Bustools_opt &opt) {
 
     for (size_t i = 0; i < rc; i++) {
       if (curr_bc != p[i].barcode) {
-        if (whitelist.find(curr_bc & len_mask) != whitelist.end()) {
+        if (in_onlist(curr_bc & len_mask)) {
           reads_wl += readsPerBc_count;
         }
 
@@ -154,7 +228,7 @@ void bustools_inspect(Bustools_opt &opt) {
         umisPerBc.push_back(umisPerBc_count);
         umisPerBc_count = 1;
 
-        if (whitelist.find(p[i].barcode & len_mask) != whitelist.end()) {
+        if (in_onlist(p[i].barcode & len_mask)) {
           ++bc_wl;
         }
       } else if (curr_umi != p[i].UMI) {
@@ -204,7 +278,7 @@ void bustools_inspect(Bustools_opt &opt) {
   /* Done reading BUS file. */
 
   // Some stats have stragglers
-  if (whitelist.find(curr_bc & len_mask) != whitelist.end()) {
+  if (in_onlist(curr_bc & len_mask)) {
     reads_wl += readsPerBc_count;
   }
   umisPerBc.push_back(umisPerBc_count);