diff --git a/src/bustools_inspect.cpp b/src/bustools_inspect.cpp index 80dcd38..fc40cb8 100644 --- a/src/bustools_inspect.cpp +++ b/src/bustools_inspect.cpp @@ -50,22 +50,83 @@ void bustools_inspect(Bustools_opt &opt) { } } - /* Load whitelist. */ - std::unordered_set 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> wbc; // Each set contains the on-list of one column + std::vector 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 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);