Changeset: 7cf718a87aee for MonetDB
URL: http://dev.monetdb.org/hg/MonetDB?cmd=changeset;node=7cf718a87aee
Modified Files:
sql/backends/monet5/bam/bam_loader.c
sql/backends/monet5/bam/bam_loader.h
Branch: DVframework_bam
Log Message:
Improved _process_bam_alignment method greatly. First, Samtools was used to
transform an alignment structure into an ASCII string. The bam_loader would
then parse this string. Since this approach was way too lame, I now retrieve
all the values from the alignment structure directly and tested the approach on
a toy file.
diffs (truncated from 426 to 300 lines):
diff --git a/sql/backends/monet5/bam/bam_loader.c
b/sql/backends/monet5/bam/bam_loader.c
--- a/sql/backends/monet5/bam/bam_loader.c
+++ b/sql/backends/monet5/bam/bam_loader.c
@@ -1,4 +1,6 @@
#include "monetdb_config.h"
+#include "bam.h"
+#include "kstring.h"
#include "bam_loader.h"
/*
@@ -92,8 +94,8 @@ static bit _append_option_to_bat(_temp_c
static str _read_bam_header_line(str *header, _bam_header_line *ret_hl, bit
*eof);
static void _free_bam_header_line(_bam_header_line *hl);
static str _process_bam_alignment(sht file_id, lng virtual_offset,
bam_header_t *header, bam1_t *alignment, _temp_container *ret_tc, sht schema);
-static bit _parse_alignment_str(str *sam_alig, str *dest);
-static bit _parse_alignment_lng(str *sam_alig, lng *dest);
+/*static bit _parse_alignment_str(str *sam_alig, str *dest);
+static bit _parse_alignment_lng(str *sam_alig, lng *dest);*/
/* Generic functions */
static str _init_temp_subcontainer(_temp_subcontainer *ret_tsc, str
*col_names, int *col_types, sht num_cols);
@@ -230,7 +232,7 @@ static str
if(ret_tc->table_names == NULL || ret_tc->num_columns == NULL)
throw(MAL, "init_temp_container", MAL_MALLOC_FAIL);
- for(t = 0; t < num_tables; t++)
+ for(t = 0; t < num_tables; ++t)
{
ret_tc->table_names[t] = table_names[t];
ret_tc->num_columns[t] = num_cols[t];
@@ -301,7 +303,7 @@ static str
{
if((err = _process_bam_alignment(file_id, voffset, header, alignment,
ret_tc, dbschema)) != MAL_SUCCEED)
throw(MAL, "loadfile", "Error processing bam alignment: %s\n",
err);
- num_alignments++;
+ ++num_alignments;
voffset = bam_tell(input);
}
duration_sec = (GDKms() - start) / 1000.0f;
@@ -359,7 +361,7 @@ static str
/* read and interpret the header tag */
if(strcmp(hl->header_tag, "HD") == 0)
{
- num_hd_lines++;
+ ++num_hd_lines;
if(num_hd_lines > 1)
{
_free_bam_header_line(hl);
@@ -369,7 +371,7 @@ static str
hd_fields_found[0] = FALSE;
hd_fields_found[1] = FALSE;
- for(o = 0; o < hl->num_options; o++)
+ for(o = 0; o < hl->num_options; ++o)
{
if(_append_option_to_bat_cond_flt(ret_tc, &hl->options[o],
"VN", 0, 2, &appendErr, &hd_fields_found[0])) continue;
if(_append_option_to_bat_cond_str(ret_tc, &hl->options[o],
"SO", 0, 3, &appendErr, &hd_fields_found[1])) continue;
@@ -399,10 +401,10 @@ static str
throw(MAL, "_process_bam_header", "Appending file_id to SQ BAT
failed\n");
}
- for(i=0; i<6; i++)
+ for(i=0; i<6; ++i)
sq_fields_found[i] = FALSE;
- for(o = 0; o < hl->num_options; o++)
+ for(o = 0; o < hl->num_options; ++o)
{
if(_append_option_to_bat_cond_str(ret_tc, &hl->options[o],
"SN", 1, 0, &appendErr, &sq_fields_found[0])) continue;
if(_append_option_to_bat_cond_lng(ret_tc, &hl->options[o],
"LN", 1, 2, &appendErr, &sq_fields_found[1])) continue;
@@ -446,10 +448,10 @@ static str
throw(MAL, "_process_bam_header", "Appending file_id to RG BAT
failed\n");
}
- for(i=0; i<12; i++)
+ for(i=0; i<12; ++i)
rg_fields_found[i] = FALSE;
- for(o = 0; o < hl->num_options; o++)
+ for(o = 0; o < hl->num_options; ++o)
{
if(_append_option_to_bat_cond_str(ret_tc, &hl->options[o],
"ID", 2, 0, &appendErr, &rg_fields_found[0])) continue;
if(_append_option_to_bat_cond_str(ret_tc, &hl->options[o],
"CN", 2, 2, &appendErr, &rg_fields_found[1])) continue;
@@ -501,10 +503,10 @@ static str
throw(MAL, "_process_bam_header", "Appending file_id to PG BAT
failed\n");
}
- for(i=0; i<5; i++)
+ for(i=0; i<5; ++i)
pg_fields_found[i] = FALSE;
- for(o = 0; o < hl->num_options; o++)
+ for(o = 0; o < hl->num_options; ++o)
{
if(_append_option_to_bat_cond_str(ret_tc, &hl->options[o],
"ID", 3, 0, &appendErr, &pg_fields_found[0])) continue;
if(_append_option_to_bat_cond_str(ret_tc, &hl->options[o],
"PN", 3, 2, &appendErr, &pg_fields_found[1])) continue;
@@ -698,7 +700,7 @@ static void
if(hl == NULL)
return;
- for(o=0; o<hl->num_options; o++)
+ for(o=0; o<hl->num_options; ++o)
{
/* for some reason, freeing hl->options[o].value causes segmentation
fault, while these
strings are created with GDKmalloc... */
@@ -708,36 +710,109 @@ static void
GDKfree(hl);
}
+/*
+* Given a Samtools native structure bam1_t, retrieve all information from it
and append it to the _temp_container.
+* Retrieval code is inspired by the code found in bam.c::bam_format1_core
+*
+* Old implementation actually called bam.c::bam_format1_core and then parsed
the resulting SAM-string. This new and less lame
+* implementation has been compared with the old implementation on the file
toy.bam and resulted in the exact same database contents.
+*/
static str
_process_bam_alignment(sht file_id, lng virtual_offset, bam_header_t *header,
bam1_t *alignment, _temp_container *ret_tc, sht dbschema)
{
- str err = NULL;
str qname, rname, cigar, rnext, seq, qual;
- lng flag, pos, mapq, pnext, tlen;
+ sht flag, mapq;
+ int pos, pnext, tlen;
- str tag = NULL, type = NULL, value = NULL; /* auxiliary options */
-
- str sam_alig;
- str parse_err = NULL;
+ uint8_t *s;
+ int i;
(void)dbschema;
- sam_alig = bam_format1(header, alignment);
- if( !_parse_alignment_str(&sam_alig, &qname ))
parse_err = "QNAME";
- if(parse_err == NULL && !_parse_alignment_lng(&sam_alig, &flag ))
parse_err = "FLAG" ;
- if(parse_err == NULL && !_parse_alignment_str(&sam_alig, &rname ))
parse_err = "RNAME";
- if(parse_err == NULL && !_parse_alignment_lng(&sam_alig, &pos ))
parse_err = "POS" ;
- if(parse_err == NULL && !_parse_alignment_lng(&sam_alig, &mapq ))
parse_err = "MAPQ" ;
- if(parse_err == NULL && !_parse_alignment_str(&sam_alig, &cigar ))
parse_err = "CIGAR";
- if(parse_err == NULL && !_parse_alignment_str(&sam_alig, &rnext ))
parse_err = "RNEXT";
- if(parse_err == NULL && !_parse_alignment_lng(&sam_alig, &pnext ))
parse_err = "PNEXT";
- if(parse_err == NULL && !_parse_alignment_lng(&sam_alig, &tlen ))
parse_err = "TLEN" ;
- if(parse_err == NULL && !_parse_alignment_str(&sam_alig, &seq ))
parse_err = "SEQ" ;
- if(parse_err == NULL && !_parse_alignment_str(&sam_alig, &qual ))
parse_err = "QUAL" ;
+ //qname
+ qname = bam1_qname(alignment);
+
+ //flag
+ flag = alignment->core.flag;
+
+ //rname
+ if(alignment->core.tid < 0)
+ rname = "*";
+ else
+ rname = header->target_name[alignment->core.tid];
+
+ //pos
+ pos = alignment->core.pos + 1;
+
+ //mapq
+ mapq = alignment->core.qual;
+
+ //cigar
+ if(alignment->core.n_cigar == 0)
+ cigar = "*";
+ else
+ {
+ uint32_t *cigar_bin = bam1_cigar(alignment);
+
+ /* Use kstring library as is done in Samtools, to conveniently
reconstruct cigar string */
+ kstring_t cigar_stream;
+ cigar_stream.l = cigar_stream.m = 0; cigar_stream.s = 0;
+
+ for (i=0; i<alignment->core.n_cigar; ++i)
+ {
+ kputw(cigar_bin[i]>>BAM_CIGAR_SHIFT, &cigar_stream);
+ kputc(bam_cigar_opchr(cigar_bin[i]), &cigar_stream);
+ }
+ cigar = (str)GDKmalloc(alignment->core.n_cigar * sizeof(char));
+ if(cigar == NULL)
+ throw(MAL, "_process_bam_alignment", MAL_MALLOC_FAIL);
+ cigar = cigar_stream.s;
+ }
+
+ //rnext
+ if(alignment->core.mtid < 0)
+ rnext = "*";
+ else if(alignment->core.mtid == alignment->core.tid)
+ rnext = "=";
+ else
+ rnext = header->target_name[alignment->core.mtid];
+
+ //pnext
+ pnext = alignment->core.mpos + 1;
+
+ //tlen
+ tlen = alignment->core.isize;
+
+ //seq and qual
+ if(alignment->core.l_qseq)
+ {
+ s = bam1_seq(alignment);
+ seq = (str)GDKmalloc((alignment->core.l_qseq + 1) * sizeof(char)); /*
+1 for \0 character */
+ if(seq == NULL)
+ throw(MAL, "_process_bam_alignment", MAL_MALLOC_FAIL);
+ for (i=0; i<alignment->core.l_qseq; ++i)
+ seq[i] = bam_nt16_rev_table[bam1_seqi(s, i)];
+ seq[alignment->core.l_qseq] = '\0';
+ s = bam1_qual(alignment);
+ if (s[0] == 0xff)
+ qual = "*";
+ else
+ {
+ qual = (str)GDKmalloc((alignment->core.l_qseq + 1) *
sizeof(char)); /* +1 for \0 character */
+ if(qual == NULL)
+ throw(MAL, "_process_bam_alignment", MAL_MALLOC_FAIL);
+ for (i=0; i<alignment->core.l_qseq; ++i)
+ qual[i] = s[i] + 33;
+ qual[alignment->core.l_qseq] = '\0';
+ }
+ }
+ else
+ {
+ seq = "*";
+ qual = "*";
+ }
- if(parse_err != NULL)
- throw(MAL, "_process_bam_alignment", "Parse error on field: %s", err);
if(BUNappend(ret_tc->tables_columns[4].column_bats[0] , (ptr)
&virtual_offset , TRUE) == NULL
|| BUNappend(ret_tc->tables_columns[4].column_bats[1] , (ptr) &file_id
, TRUE) == NULL
@@ -759,48 +834,68 @@ static str
/* parse auxiliary data */
- /* tabs have been stripped already, so we are directly faced with an
optional field, a \n or a \0 */
- while(TRUE)
- {
- if(*sam_alig == '\n' || *sam_alig == '\0') break;
-
- if(_read_string_until_delim(&sam_alig, &tag, ":\t\n\0", 4) != 2)
- throw(MAL, "_process_bam_alignment", "Auxiliary data tag length !=
2 found in alignment: %s", tag);
- if(*sam_alig != ':')
- throw(MAL, "_process_bam_alignment", "Auxiliary data tag contains
one of '\\t\\n\\0'");
- sam_alig++;
-
- if(_read_string_until_delim(&sam_alig, &type, ":\t\n\0", 4) != 1)
- throw(MAL, "_process_bam_alignment", "Auxiliary data type length
!= 1 found in alignment: %s:%s", tag, type);
- if(*sam_alig != ':')
- throw(MAL, "_process_bam_alignment", "Auxiliary data type contains
one of '\\t\\n\\0'");
- sam_alig++;
-
- if(_read_string_until_delim(&sam_alig, &value, "\t\n\0", 3) == -1)
- throw(MAL, "_process_bam_alignment", MAL_MALLOC_FAIL);
-
- if(BUNappend(ret_tc->tables_columns[5].column_bats[0] , (ptr) tag
, TRUE) == NULL
- || BUNappend(ret_tc->tables_columns[5].column_bats[1] , (ptr)
&virtual_offset , TRUE) == NULL
- || BUNappend(ret_tc->tables_columns[5].column_bats[2] , (ptr) &file_id
, TRUE) == NULL
- || BUNappend(ret_tc->tables_columns[5].column_bats[3] , (ptr) type
, TRUE) == NULL
- || BUNappend(ret_tc->tables_columns[5].column_bats[4] , (ptr) value
, TRUE) == NULL
+ s = bam1_aux(alignment);
+ while (s < alignment->data + alignment->data_len) {
+ uint8_t type, tag[2];
+ kstring_t aux_value_stream;
+ aux_value_stream.l = aux_value_stream.m = 0; aux_value_stream.s = 0;
+ tag[0] = s[0];
+ tag[1] = s[1];
+ s += 2;
+ type = *s;
+ ++s;
+ if (type == 'A') { kputc(*s
, &aux_value_stream) ; ++s; }
+ else if (type == 'C') { type = 'i'; kputw(*s
, &aux_value_stream) ; ++s; }
+ else if (type == 'c') { type = 'i';
kputw(*(int8_t*)s , &aux_value_stream) ; ++s; }
+ else if (type == 'S') { type = 'i';
kputw(*(uint16_t*)s , &aux_value_stream) ; s += 2; }
+ else if (type == 's') { type = 'i';
kputw(*(int16_t*)s , &aux_value_stream) ; s += 2; }
+ else if (type == 'I') { type = 'i';
kputuw(*(uint32_t*)s, &aux_value_stream) ; s += 4; }
+ else if (type == 'i') { type = 'i';
kputw(*(int32_t*)s , &aux_value_stream) ; s += 4; }
+ else if (type == 'f') {
ksprintf(&aux_value_stream, "%g", *(float*)s) ; s += 4; }
+ else if (type == 'd') {
ksprintf(&aux_value_stream, "%lg", *(double*)s); s += 8; }
+ else if (type == 'Z' || type == 'H') { while (*s)
kputc(*s++, &aux_value_stream) ; ++s; }
+ else if (type == 'B')
+ {
+ uint8_t sub_type = *(s++);
+ int32_t n;
+ memcpy(&n, s, 4);
+ s += 4; // no point to the start of the array
+ kputc(sub_type, &aux_value_stream); // write the typing
+ for (i = 0; i < n; ++i)
+ {
+ kputc(',', &aux_value_stream);
+ if ('c' == sub_type) { kputw(*(int8_t*)s,
&aux_value_stream) ; ++s; }
+ else if ('C' == sub_type) { kputw(*(uint8_t*)s,
&aux_value_stream) ; ++s; }
+ else if ('s' == sub_type) { kputw(*(int16_t*)s,
&aux_value_stream) ; s += 2; }
+ else if ('S' == sub_type) {
kputw(*(uint16_t*)s, &aux_value_stream) ; s += 2; }
+ else if ('i' == sub_type) { kputw(*(int32_t*)s,
&aux_value_stream) ; s += 4; }
+ else if ('I' == sub_type) {
kputuw(*(uint32_t*)s, &aux_value_stream) ; s += 4; }
_______________________________________________
checkin-list mailing list
[email protected]
http://mail.monetdb.org/mailman/listinfo/checkin-list
