Question: I want to write some c code to convert fastq to unmapped bam file. I use htslib to read and write bam file and have some problem. The bam file I generated seems not validated. I guess I missed some part but I dont know what it is. This is my toy code:

// trim_fq_to_bam.c
#include // for reading compressed .fq file
#include
#include
#include
#include
#include
#include
KSEQ_INIT(gzFile, gzread)
#define bam1_seq_seti(s, i, c) ( (s)[(i)>>1] = ((s)[(i)>>1] & 0xf<<(((i)&1)<<2)) | (c)<<((~(i)&1)<<2) )

static const char init_header[] = "@HD\tVN:1.4\tSO:unknown\n";

void dump_read(bam1_t* b) {
printf("->core.tid:(%d)\n", b->core.tid);
printf("->core.pos:(%d)\n", b->core.pos);
printf("->core.bin:(%d)\n", b->core.bin);
printf("->core.qual:(%d)\n", b->core.qual);
printf("->core.l_qname:(%d)\n", b->core.l_qname);
printf("->core.flag:(%d)\n", b->core.flag);
printf("->core.n_cigar:(%d)\n", b->core.n_cigar);
printf("->core.l_qseq:(%d)\n", b->core.l_qseq);
printf("->core.mtid:(%d)\n", b->core.mtid);
printf("->core.mpos:(%d)\n", b->core.mpos);
printf("->core.isize:(%d)\n", b->core.isize);
if (b->data) {
printf("->data:");
int i;
for (i = 0; i < b->l_data; ++i) {
printf("%x ", b->data[i]);
}
printf("\n");
}
if (b->core.l_qname) {
printf("qname: %s\n",bam_get_qname(b));
}
if (b->core.l_qseq) {
printf("qseq:");
int i;
for (i = 0; i < b->core.l_qseq; ++i) {
printf("%c", seq_nt16_str[bam_seqi(bam_get_seq(b),i)]);
}
printf("\n");
printf("qual:");
for (i = 0; i < b->core.l_qseq; ++i) {
printf("%c",bam_get_qual(b)[i]);
}
printf("\n");

}

if (bam_get_l_aux(b)) {
int i = 0;
uint8_t* aux = bam_get_aux(b);

while (i < bam_get_l_aux(b)) {
printf("%.2s:%c:",aux+i,*(aux+i+2));
i += 2;
switch (*(aux+i)) {
case 'Z':
while (*(aux+1+i) != '\0') { putc(*(aux+1+i), stdout); ++i; }
break;
}
putc('\n',stdout);
++i;++i;
}
}
printf("\n");
}

int main(int argc, char const *argv[])
{
char *fn_in = "test_data/rd1.fq.gz";
char *fn_out = "test_data/new.bam";
int l = 0;
gzFile *f_rd = gzopen(fn_in, "r");
if (!f_rd){fprintf(stderr, "cant open file: %s\n", fn_in); exit(EXIT_FAILURE);}
kseq_t *seq;
seq = kseq_init(f_rd);
BGZF *fp = bgzf_open(fn_out,"w");

// write header
bam_hdr_t *hdr = bam_hdr_init();
hdr->l_text = strlen(init_header);
hdr->text = strdup(init_header);
hdr->n_targets = 0;
bam_hdr_write(fp, hdr);

// main loop, iter through each fastq records
while ((l = kseq_read(seq)) >= 0){
printf("%s, %s\n", seq->seq.s, seq->qual.s);
bam1_t *q = bam_init1();
//`q->data` structure: qname-cigar-seq-qual-aux
q->l_data = seq->name.l+1+(int)(1.5*seq->seq.l+(seq->seq.l % 2 != 0)); // +1 includes the tailing '\0'
if (q->m_data < q->l_data) {
q->m_data = q->l_data;
kroundup32(q->m_data);
q->data = (uint8_t*)realloc(q->data, q->m_data);
}
q->core.flag = BAM_FMUNMAP;
q->core.l_qname = seq->name.l+1; // +1 includes the tailing '\0'
q->core.l_qseq = seq->seq.l;
q->core.n_cigar = 0; // we have no cigar sequence
memcpy(q->data, seq->name.s, q->core.l_qname); // first set qname
uint8_t *s = bam_get_seq(q);
for (int i = 0; i < q->core.l_qseq; ++i){

bam1_seq_seti(s, i, seq_nt16_table[seq->seq.s[i]]);
}
s = bam_get_qual(q);
for (int i = 0; i < seq->qual.l; ++i) {
s[i] = seq->qual.s[i];
}
dump_read(q);
int ret = bam_write1(fp, q);
printf("return value: %d\n", ret);
bam_destroy1(q);
}
kseq_destroy(seq); // free seq
gzclose(f_rd); // close fastq file
bgzf_close(fp); // close bam file

// validate bam file
BGZF *fp_new = bgzf_open(fn_out,"r");
bam1_t *q = bam_init1();
int ret = bam_read1(fp_new, q);
printf("%d\n", ret);
printf("%d, %d\n", (char *)bam_get_aux(q) - (char *)q->data , q->l_data);
dump_read(q);
return 0;
}

the function `dump_read()` is used to print a bam record to see if I have successfully create a valid `bam1_t` structure, and the code is copied from here:
https://github.com/samtools/samtools/blob/74a3e0a744f91bc2b1f65a4144a509d6dd791631/test/merge/test_bam_translate.c
with some minor modification. For printing the sequence the original code is `printf("%c",seq_nt16_str[seq_nt16_table[bam_seqi(bam_get_seq(b),i)]]);` . but I didnt get the point of jointly using `seq_nt16_table[]` and `seq_nt16_str[]`, since `bam_seqi()` returns the 4 bit representation of a base and `seq_nt16_str` convert it to str so these two would be sufficient.

The basic structure of the code is:

open fq file and bam file ->

write bam header (I am not sure if I have done it correctly) ->

iterate through each fastq reads, convert it to a `bam1_t` structure ->

use `dump_read()` to check if the `bam1_t` is correct. ->

write `bam1_t` to the bam file. check the return to see if the writing is successful.

So far everything looks fine. the `dump_read()` print the same seq and quality strings as the fastq reads. `bam_write1()` gives positive returns. And I generated a bam file. But when I try `samtools view -h new.bam`, it gives an error `[main_samview] truncated file.` Then I hacked the code and find when `sam_read1() < -1` it will print this err msg. So I reopened my bam file and calls `bam_read1(fp_new, q);`, the return is -4, I checked the source code to see when it will returns -4 and stuck. I don't know which part of my code is wrong.

Also I wonder if there are any toturial on how to use htslib, or any examples so I can start from there. I usually work on Python and I am not familiar with C. Now I just read the source code of htslib but it is a bit messy, with a lot of uncomment blocks.

Answer: finally I solved the problem. Yes `bgzf_open(fn_out,"w")` is a lower level functions and one should use `sam_open(fn_out,"wb")` to write to a binary sam file (which is bam). But just change the open function does not fix everything. I also need to add these lines to set some variables to "NULL".

q->core.tid = -1;
q->core.pos = -1;
q->core.mtid = -1;
q->core.mpos = -1;

as an example, `tid` is the chromosome index in a char array which start at zero. So for unaligned bam where we have no chromosome, we need to set it to -1.

Thank you @Devon Ryan.