Changeset: c493be8ab1e2 for MonetDB URL: https://dev.monetdb.org/hg/MonetDB/rev/c493be8ab1e2 Added Files: sql/backends/monet5/vaults/sambam/CMakeLists.txt sql/backends/monet5/vaults/sambam/sambam.c sql/test/file_loader/sambam/Tests/All sql/test/file_loader/sambam/Tests/SingleServer sql/test/file_loader/sambam/Tests/file1.bam sql/test/file_loader/sambam/Tests/file1.sam sql/test/file_loader/sambam/Tests/read_file1.bam.test.in sql/test/file_loader/sambam/Tests/read_file1.sam.test.in Modified Files: sql/backends/monet5/vaults/CMakeLists.txt testing/Mtest.py.in Branch: samloader Log Message:
Add first working version of loader for *.sam and *.bam files using htslib. The module is called sambam and needs to be loaded explicitly using mserver5 arg: --loadmodule sambam There are still some TODO's in sambam.c diffs (truncated from 1564 to 300 lines): diff --git a/sql/backends/monet5/vaults/CMakeLists.txt b/sql/backends/monet5/vaults/CMakeLists.txt --- a/sql/backends/monet5/vaults/CMakeLists.txt +++ b/sql/backends/monet5/vaults/CMakeLists.txt @@ -13,6 +13,7 @@ add_subdirectory(netcdf) add_subdirectory(shp) add_subdirectory(csv) add_subdirectory(parquet) +add_subdirectory(sambam) add_subdirectory(monetdb) add_subdirectory(odbc) diff --git a/sql/backends/monet5/vaults/sambam/CMakeLists.txt b/sql/backends/monet5/vaults/sambam/CMakeLists.txt new file mode 100644 --- /dev/null +++ b/sql/backends/monet5/vaults/sambam/CMakeLists.txt @@ -0,0 +1,60 @@ +#[[ +# SPDX-License-Identifier: MPL-2.0 +# +# This Source Code Form is subject to the terms of the Mozilla Public +# License, v. 2.0. If a copy of the MPL was not distributed with this +# file, You can obtain one at https://mozilla.org/MPL/2.0/. +# +# For copyright information, see the file debian/copyright. +#]] + +if(HAVE_HTSLIB) + add_library(sambam MODULE) + + target_sources(sambam + PRIVATE + sambam.c) + + target_include_directories(sambam + PRIVATE + $<$<BOOL:${HTSLIB_FOUND}>:${HTSLIB_INCLUDE_DIRS}> + $<TARGET_PROPERTY:mal,INTERFACE_INCLUDE_DIRECTORIES> + $<TARGET_PROPERTY:malmodules,INTERFACE_INCLUDE_DIRECTORIES> + $<TARGET_PROPERTY:atoms,INTERFACE_INCLUDE_DIRECTORIES> + $<TARGET_PROPERTY:sql,INTERFACE_INCLUDE_DIRECTORIES> + $<TARGET_PROPERTY:sqlcommon,INTERFACE_INCLUDE_DIRECTORIES> + $<TARGET_PROPERTY:sqlserver,INTERFACE_INCLUDE_DIRECTORIES> + $<TARGET_PROPERTY:sqlstorage,INTERFACE_INCLUDE_DIRECTORIES> + $<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}> + $<INSTALL_INTERFACE:${INCLUDEDIR}/monetdb>) + + set_target_properties(sambam + PROPERTIES + NO_SONAME ON + OUTPUT_NAME + _sambam) + + target_compile_definitions(sambam + PRIVATE + LIBHTSLIB) + + target_link_libraries(sambam + PRIVATE + monetdb_config_header + bat + monetdb5 + sqlinclude + HTSLIB::HTSLIB) + + install(TARGETS sambam + DESTINATION ${CMAKE_INSTALL_LIBDIR}/monetdb5-${MONETDB_VERSION} + COMPONENT sambam) + + if(WIN32) + install(FILES + $<TARGET_PDB_FILE:sambam> + DESTINATION ${CMAKE_INSTALL_LIBDIR}/monetdb5-${MONETDB_VERSION} + OPTIONAL) + endif() +endif() + diff --git a/sql/backends/monet5/vaults/sambam/sambam.c b/sql/backends/monet5/vaults/sambam/sambam.c new file mode 100644 --- /dev/null +++ b/sql/backends/monet5/vaults/sambam/sambam.c @@ -0,0 +1,480 @@ +/* + * SPDX-License-Identifier: MPL-2.0 + * + * This Source Code Form is subject to the terms of the Mozilla Public + * License, v. 2.0. If a copy of the MPL was not distributed with this + * file, You can obtain one at https://mozilla.org/MPL/2.0/. + * + * For copyright information, see the file debian/copyright. + */ + +/* + * This file contains a reader of aligmnent data from *.sam and *.bam files + * and load it into a transient MonetDB table consisting of 11 columns. + * SAM stands for Sequence Alignment Map. Ref: https://en.wikipedia.org/wiki/SAM_(file_format) + * BAM stands for Binary Alignment Map. Ref: https://en.wikipedia.org/wiki/BAM_(file_format) + * BAM is the compressed binary representation of SAM (Sequence Alignment Map). + * BAM is one of the most common file formats for raw genome sequencing data. + * + * A SAM file is a TAB-delimited text format consisting of a header section, which is optional, + * and an alignment section. If present, the header must be prior to the alignments. + * Header lines start with ‘@’, while alignment lines do not. + * Each alignment line has 11 mandatory fields for essential alignment information + * such as mapping position, and variable number of optional fields for flexible + * or aligner specific information. + * Ref: https://samtools.github.io/hts-specs/SAMv1.pdf + * + * We use the htslib library for reading data from *.sam and *.bam files. + * Ref: https://github.com/samtools/htslib + * To install this library and its include files, follow instructions at: + * https://github.com/samtools/htslib/blob/develop/INSTALL + * Note: when you issue last step 'make install' you need root privileges so do a 'sudo make install' instead. + * + * Limitations: + * - Only the alignment data is converted into bats, not the header information as it has a different table format. + * - We currently only read the first 11 mandatory columns. Thus ignoring any additonal columns of info. + * - The topn filter is not yet implemented. + * - *.cram files are not yet supported. + * + * Author: Martin van Dinther + */ + +#include "monetdb_config.h" +#include "gdk.h" // COLnew(), bunfastapp() +#include "mal_builder.h" // newStmtArgs(), pushStr() +#include "rel_exp.h" // exp_column() +#include "rel_file_loader.h" // fl_register(), fl_unregister() +#include "sql_monet_backend.h" // backend, getBackendContext() + +#include <htslib/sam.h> // it includes hts.h +/* note there is no <htslib/bam.h> */ + +/* copied from monetdb5/modules/mal/tablet.c */ +static BAT * +bat_create(int adt, BUN nr) +{ + BAT *b = COLnew(0, adt, nr, TRANSIENT); + + /* check for correct structures */ + if (b == NULL) + return NULL; + if ((b = BATsetaccess(b, BAT_APPEND)) == NULL) { + return NULL; + } + + /* disable all properties here */ + b->tsorted = false; + b->trevsorted = false; + b->tnosorted = 0; + b->tnorevsorted = 0; + b->tseqbase = oid_nil; + b->tkey = false; + b->tnokey[0] = 0; + b->tnokey[1] = 0; + b->tnil = false; + b->tnonil = false; + return b; +} + +static void +sam_list_append(list * nameslist, char * name, list * typelist, sql_subtype * mtype, list *res_exps, mvc *sql, int has_nils) +{ + list_append(nameslist, name); + list_append(typelist, mtype); + if (res_exps) { + sql_exp *ne = exp_column(sql->sa, "", name, mtype, CARD_MULTI, has_nils, 0, 0); + set_basecol(ne); + ne->alias.label = -(sql->nid++); + list_append(res_exps, ne); + } +} + +/* nr of sam/bam file alignment info columns */ +#define SAM_NR_COLS 11 +#define SAM_MAX_STR_SIZE 4096 + +#define SAM_RELATION 1 +#define SAM_LOADER 2 + +/* + * sam_query() contains the logic for both sam_relation() and SAMloader() functions. + * the caller argument is SAM_RELATION when called from sam_relation() or SAM_LOADER when called from SAMloader(). + */ +static str +sam_query(int caller, mvc *sql, sql_subfunc *f, char *filename, list *res_exps, MalBlkPtr mb, MalStkPtr stk, InstrPtr pci) +{ + assert(caller == SAM_RELATION || caller == SAM_LOADER); + + if (sql == NULL) + return "Missing mvc value"; + if (f == NULL) + return "Missing sql_subfunc value"; + if (filename == NULL) + return "Missing filename value"; + + bam1_t *bamdata = NULL; + sam_hdr_t *in_samhdr = NULL; + samFile *infile = NULL; + char * errmsg = NULL; + + if (!(bamdata = bam_init1())) { + errmsg = "baminit1() failed to allocate memory"; + goto end; + } + + // open input file + if (!(infile = sam_open(filename, "r"))) { + /* TODO sam_open() logs error msg to stderr when it fails: + [E::hts_open_format] Failed to open file "file1.sam" : No such file or directory + We need to prevent this logging to stderr. + */ + errmsg = "Could not open sam/bam file"; + /* errmsg = RUNTIME_FILE_NOT_FOUND; */ + goto end; + } + + // Read a header from a SAM, BAM or CRAM file. + if (!(in_samhdr = sam_hdr_read(infile))) { + errmsg = "Failed to read sam/bam file header"; + goto end; + } + + /* Alignment data section in a sam/bam/cram file defines 11 mandatory columns + Ref: https://samtools.github.io/hts-specs/SAMv1.pdf secction 1.4 + Col Field Type Regexp/Range Brief description + 1 QNAME String [!-?A-~]{1,254} Query template NAME + 2 FLAG Int [0, 2^16 − 1] bitwise FLAG + 3 RNAME String \*|[:rname:∧*=][:rname:]* Reference sequence NAME12 + 4 POS Int [0, 2^31 − 1] 1-based leftmost mapping POSition + 5 MAPQ Int [0, 2^8 − 1] MAPping Quality + 6 CIGAR String \*|([0-9]+[MIDNSHP=X])+ CIGAR string + 7 RNEXT String \*|=|[:rname:∧*=][:rname:]* Reference name of the mate/next read + 8 PNEXT Int [0, 2^31 − 1] Position of the mate/next read + 9 TLEN Int [−2^31 + 1, 2^31 − 1] observed Template LENgth + 10 SEQ String \*|[A-Za-z=.]+ segment SEQuence + 11 QUAL String [!-~]+ ASCII of Phred-scaled base QUALity+33 + */ + + /* TODO extend with 2 (or more) optional fields (all varchar4096), + Ref: https://samtools.github.io/hts-specs/SAMv1.pdf secction 1.5 + Col Field Type Regexp/Range Brief description + 12 SAM_AUX String + 13 SAM_RGAUX String + */ + + /* when called from sam_relation() */ + if (caller == SAM_RELATION) { + sql_subtype * varchar255 = sql_bind_subtype(sql->sa, "varchar", (unsigned int) BAM_MAX_QNAME_LEN +1, 0); + sql_subtype * varchar1020 = sql_bind_subtype(sql->sa, "varchar", (unsigned int) 1020, 0); + sql_subtype * varchar4096 = sql_bind_subtype(sql->sa, "varchar", (unsigned int) SAM_MAX_STR_SIZE, 0); + sql_subtype * int31 = sql_bind_subtype(sql->sa, "int", (unsigned int) 0, 0); + sql_subtype * int15 = sql_bind_subtype(sql->sa, "smallint", (unsigned int) 0, 0); + + list * typelist = sa_list(sql->sa); + list * nameslist = sa_list(sql->sa); + + sam_list_append(nameslist, "qname", typelist, varchar255, res_exps, sql, 0); + sam_list_append(nameslist, "flag", typelist, int31, res_exps, sql, 0); + sam_list_append(nameslist, "rname", typelist, varchar1020, res_exps, sql, 0); + sam_list_append(nameslist, "pos", typelist, int31, res_exps, sql, 0); + sam_list_append(nameslist, "mapq", typelist, int15, res_exps, sql, 0); + sam_list_append(nameslist, "cigar", typelist, varchar4096, res_exps, sql, 0); + sam_list_append(nameslist, "rnext", typelist, varchar1020, res_exps, sql, 0); + sam_list_append(nameslist, "pnext", typelist, int31, res_exps, sql, 0); + sam_list_append(nameslist, "tlen", typelist, int31, res_exps, sql, 0); + sam_list_append(nameslist, "seq", typelist, varchar4096, res_exps, sql, 0); + sam_list_append(nameslist, "qual", typelist, varchar4096, res_exps, sql, 0); + /* TODO extend with 2 (or more) optional fields, + sam_list_append(nameslist, "sam_aux", typelist, varchar4096, res_exps, sql, 1); + sam_list_append(nameslist, "sam_rgaux", typelist, varchar4096, res_exps, sql, 1); + */ + + f->tname = ma_strdup(sql->sa, ""); /* TODO maybe use the filename without path and extension */ + f->colnames = nameslist; + f->coltypes = typelist; + f->res = typelist; + goto end; + } + + /* when called from SAMloader() */ + if (caller == SAM_LOADER) { + typedef struct { + int battype; /* MonetDB atom type, used to create the BAT */ + BAT * bat; /* MonetDB BAT */ + } rescol_t; + + rescol_t * colmetadata = (rescol_t *) GDKmalloc(SAM_NR_COLS * sizeof(rescol_t)); + if (colmetadata == NULL) { + errmsg = "GDKzalloc colmetadata[SAM_NR_COLS] failed."; + goto end; + } + + assert(mb); + assert(pci); + int col; + /* make bats with right atom type */ + for (col = 0; col < SAM_NR_COLS; col++) { + int battype = getBatType(getArgType(mb, pci, col)); + TRC_DEBUG(LOADER, "Before create BAT %d type %d\n", col+1, battype); + colmetadata[col].battype = battype; _______________________________________________ checkin-list mailing list -- [email protected] To unsubscribe send an email to [email protected]
