How to read .bwt file .sa file?

Question

In BWA project in bwt.c file I found bwt_dump_bwt method and bwt_dump_sa method. I want to use them to read .bwt file and .sa file. I wrote the following program:

int main(void)
{
  char *tmp, *prefix;
    bwt_t *bwt;
    prefix = bwa_idx_infer_prefix("hg19bwaidx");

    tmp = calloc(strlen(prefix) + 5, 1);
    strcat(strcpy(tmp, prefix), ".bwt");
    bwt = bwt_restore_bwt(tmp);

  return 0;
}

I got the following error:

Undefined symbols for architecture x86_64:
  "_bwa_idx_infer_prefix", referenced from:
      _main in readbwt-3022ec.o
  "_bwt_restore_bwt", referenced from:
      _main in readbwt-3022ec.o
ld: symbol(s) not found for architecture x86_64
clang: error: linker command failed with exit code 1 (use -v to see invocation)

Any help ?

Thanks.

Answer 1

In the bwa source code, file bwt.{c,h} implement basic operations on .bwt and .sa files. They depend on utils.{c,h}, but you can create a fake utils.h to cut this dependency as is shown below. Then you only need three files: kvec.h, bwt.h and bwt.c to access .bwt and .sa.

The following example demonstrates how to find super-maximal exact matches (SMEMs). Note that internally bwa concatenates all database sequences into one long string. If you use bwt.{c,h} only, you get coordinates on this long string but not on each chromosome. You have to call APIs in bntseq.{c,h}, which are fairly independent, too. You can go on to implement a mini bwa-mem with this example, but you will need all the bwa-mem source code, except the main function.

Compiling and running

# compile the example
git clone https://github.com/lh3/bwa
mkdir -p tiny-smem && cd tiny-smem
cp ../bwa/{bwt.c,bwt.h,kvec.h} .              # you need at least these three files
### then create utils.h and test.c shown below ###
gcc -g -O2 -Wall bwt.c test.c -o tiny-smem    # compile tiny-smem

# create a bwa index if you don't have it. We use chr20 from GRCh38 as an example:
(cd ../bwa && make)                           # compile bwa for indexing
../bwa/bwa index -p chr20 path/to/human/chr20.fa

# run the example
./tiny-smem chr20.bwt chr20.sa TTTTCTTCTTTTGTAATATAAGCATTgGTTGTTCTAAATTTCTCTCTTGGCACTGCTTTA

The last command line above outputs:

0       0       26      1       -125734
0       14      27      4       -3657204        +16021887       +11802682       +54692983
0       15      29      1       +8344885
0       17      31      1       -39545379
0       19      37      1       +55060482
0       25      39      2       +1198672        -26409692
0       27      60      1       -125700

where column 2 and 3 give the coordinate of SMEM on the query, col4 the number of exact hits and the following columns are the list of hits in the concatenated coordinate.

File: utils.h

#ifndef FAKE_UTILS_H
#define FAKE_UTILS_H

#define xassert(cond, msg) assert(cond)
#define xopen(fn, mode) fopen((fn), (mode))
#define err_fclose(fp) fclose(fp)
#define err_fread_noeof(ptr, size, n, fp) fread((ptr), (size), (n), (fp))
#define err_fwrite(ptr, size, n, fp) fwrite((ptr), (size), (n), (fp))
#define err_fflush(fp) fflush(fp)
#define err_fseek(fp, off, whence) fseek((fp), (off), (whence))
#define err_ftell(fp) ftell(fp)

#endif // ~FAKE_UTILS_H

File: test.c

#include <stdlib.h>
#include <string.h>
#include <ctype.h>
#include <stdio.h>
#include "bwt.h"

int main(int argc, char *argv[]) {
  bwt_t *bwt;
  int j, max_intv = 100; // ignore SA intervals larger than this value
  bwtintv_v r = {0,0,0};

  if (argc < 4) { fprintf(stderr, "Usage: ./test <.bwt> <.sa> <seq>\n"); return 1; }
  if ((bwt = bwt_restore_bwt(argv[1])) == NULL) abort(); // read .bwt
  bwt_restore_sa(argv[2], bwt);                          // read .sa
  for (j = 3; j < argc; ++j) {
    int i, k, x, len = strlen(argv[j]);
    uint8_t *seq = (uint8_t*)malloc(len);
    for (i = 0; i < len; ++i) { // convert A/C/G/T/N => 0/1/2/3/4
      int c = tolower(argv[j][i]);
      if (c == 'a') seq[i] = 0;
      else if (c == 'c') seq[i] = 1;
      else if (c == 'g') seq[i] = 2;
      else if (c == 't') seq[i] = 3;
      else seq[i] = 4;
    }
    x = 0, r.n = 0;
    while (x < len) {
      x = bwt_smem1(bwt, len, seq, x, 1, &r, 0); // find SMEMs
      for (k = 0; k < r.n; ++k) { // traverse each SMEM
        bwtintv_t *p = &r.a[k]; // p->info encodes start and end on the query
        printf("%d\t%d\t%d\t%d", j-3, (int)(p->info>>32), (int)p->info, (int)p->x[2]);
        if (r.a[k].x[2] <= max_intv) { // don't retrieve coordinates for huge intervals
          int l_mem = (int)p->info - (int)(p->info>>32);
          for (i = 0; i < p->x[2]; ++i) { // traverse each hit
            uint64_t y = bwt_sa(bwt, p->x[0] + i); // get suffix array coordinate
            if (y < bwt->seq_len>>1) printf("\t+%ld", (long)y);
            else printf("\t-%ld", (long)(bwt->seq_len - y - l_mem));
          }
        }
        putchar('\n');
      }
    }
    free(seq);
  }
  free(r.a);
  bwt_destroy(bwt);
  return 0;
}

Answer 2

You need to include bwa.o in your call to clang, since it holds the definitions of those functions. This will cause a good number of object files, so you might as well just clang -o foo *.o foo.c -lz -lpthread -lm -lrt and get it over with. Since this was quite obviously never meant to be done, you'll have to remove at least one main() function from bwa.c.

Note that further questions regarding this will be considered off-topic and closed, since it's purely about compiling something that you want to write yourself.