How to read .bwt file .sa file?

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:
"_bwt_restore_bwt", referenced from:
clang: error: linker command failed with exit code 1 (use -v to see invocation)


Any help ?

Thanks.

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_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
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;
}


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.

• Thanks. But there is no main() in bwa.c by the way I want to do something ppl did before. like this call to bwt_dump_bwt in this link: github.com/nsoranzo/bwa/blob/master/bwtindex.c The problem is I'm new to c programming. That why I need a help and still I can't understand how to do it. Feb 12 '18 at 8:42
• There's a main in a different file then. You linked to a fork of the bwa repository, which apparently used to include that file. Feb 12 '18 at 8:44
• I ran clang -o main bwa.o main.c -lz -lpthread -lm -lrt and I got: ld: library not found for -lrt and when I removed -lrt I got many other errors. Any idea? I'm a Mac user maybe this is the problem. Feb 12 '18 at 9:18
• You'll have to figure out the equivalent library (or install it) on OS X, or whatever system you're using. Feb 12 '18 at 9:19
• I ran BWA successfully on my system. But I can't read .bwt with the program that I wrote. Feb 12 '18 at 9:21