From: Charles Plessy Date: Thu, 13 Jan 2011 09:09:43 +0000 (+0900) Subject: Imported Upstream version 0.2 X-Git-Tag: archive/raspbian/0.22.0+ds-1+rpi1~1^2^2~468 X-Git-Url: https://dgit.raspbian.org/?a=commitdiff_plain;h=209dac4a745ea4365a0738e730ca243a75c0e03f;p=python-pysam.git Imported Upstream version 0.2 --- diff --git a/INSTALL b/INSTALL index 132e50c..65330bc 100644 --- a/INSTALL +++ b/INSTALL @@ -38,3 +38,13 @@ distribution. Type python setup.py install --help for more options. + +Architecture specific options +============================= + +Pysam has been compiled on various linux systems and works +with python 2.6 and python 2.5. + +Python 2.7 and Python 3 have not been tested. + +Windows support does not work yet diff --git a/MANIFEST.in b/MANIFEST.in index c0ca312..11fb9d1 100644 --- a/MANIFEST.in +++ b/MANIFEST.in @@ -17,6 +17,9 @@ include tests/ex1.fa include tests/ex1.sam.gz include tests/ex3.sam include tests/ex4.sam +include tests/ex5.sam +include tests/ex6.sam +include tests/ex7.sam include tests/example.py include tests/pysam_test.py include tests/segfault_tests.py diff --git a/PKG-INFO b/PKG-INFO index 765534d..3e3b745 100644 --- a/PKG-INFO +++ b/PKG-INFO @@ -1,6 +1,6 @@ Metadata-Version: 1.0 Name: pysam -Version: 0.1.2 +Version: 0.2 Summary: pysam Home-page: http://code.google.com/p/pysam/ Author: Andreas Heger diff --git a/pysam/csamtools.pxd b/pysam/csamtools.pxd index f8361b7..7dac38d 100644 --- a/pysam/csamtools.pxd +++ b/pysam/csamtools.pxd @@ -27,10 +27,15 @@ cdef extern from "stdio.h": FILE * stdout int fclose(FILE *) int sscanf(char *str,char *fmt,...) + int printf(char *str,char *fmt,...) int sprintf(char *str,char *fmt,...) int fprintf(FILE *ifile,char *fmt,...) char *fgets(char *str,int size,FILE *ifile) +cdef extern from "ctype.h": + int toupper(int c) + int tolower(int c) + cdef extern from "unistd.h": char *ttyname(int fd) int isatty(int fd) @@ -43,20 +48,24 @@ cdef extern from "string.h": char *strdup(char *) char *strcat(char *,char *) size_t strlen(char *s) + int memcmp( void * s1, void *s2, size_t len ) cdef extern from "razf.h": pass -cdef extern from "bam.h": - - # IF _IOLIB=2, bamFile = BGZF, see bgzf.h - # samtools uses KNETFILE, check how this works +cdef extern from "stdint.h": ctypedef int int64_t ctypedef int int32_t ctypedef int uint32_t ctypedef int uint8_t ctypedef int uint64_t + +cdef extern from "bam.h": + + # IF _IOLIB=2, bamFile = BGZF, see bgzf.h + # samtools uses KNETFILE, check how this works + ctypedef struct tamFile: pass @@ -131,15 +140,6 @@ cdef extern from "bam.h": bam_plbuf_t *bam_plbuf_init(bam_pileup_f func, void *data) - ctypedef struct bam_fetch_iterator_t: - pass - - bam_fetch_iterator_t* bam_init_fetch_iterator(bamFile fp, bam_index_t *idx, int tid, int beg, int end) - - bam1_t * bam_fetch_iterate(bam_fetch_iterator_t *iter) - - void bam_cleanup_fetch_iterator(bam_fetch_iterator_t *iter) - int bam_fetch(bamFile fp, bam_index_t *idx, int tid, int beg, int end, void *data, bam_fetch_f func) int bam_plbuf_push(bam1_t *b, bam_plbuf_t *buf) @@ -163,12 +163,16 @@ cdef extern from "bam.h": bam1_t * bam_copy1(bam1_t *bdst, bam1_t *bsrc) uint8_t *bam_aux_get(bam1_t *b, char tag[2]) + int bam_aux2i(uint8_t *s) float bam_aux2f(uint8_t *s) double bam_aux2d(uint8_t *s) char bam_aux2A( uint8_t *s) char *bam_aux2Z( uint8_t *s) + + int bam_reg2bin(uint32_t beg, uint32_t end) + uint32_t bam_calend(bam1_core_t *c, uint32_t *cigar) cdef extern from "sam.h": @@ -192,11 +196,20 @@ cdef extern from "sam.h": int samwrite(samfile_t *fp, bam1_t *b) -cdef extern from "pysam_util.h": +cdef extern from "faidx.h": + + ctypedef struct faidx_t: + pass + + int fai_build(char *fn) + + void fai_destroy(faidx_t *fai) - ctypedef struct pair64_t: - uint64_t u - uint64_t v + faidx_t *fai_load(char *fn) + + char *fai_fetch(faidx_t *fai, char *reg, int *len) + +cdef extern from "pysam_util.h": int pysam_bam_plbuf_push(bam1_t *b, bam_plbuf_t *buf, int cont) @@ -210,3 +223,29 @@ cdef extern from "pysam_util.h": # stand-in functions for samtools macros void pysam_bam_destroy1( bam1_t * b) + + # add *nbytes* into the variable length data of *src* at *pos* + bam1_t * pysam_bam_update( bam1_t * b, + size_t nbytes_old, + size_t nbytes_new, + uint8_t * pos ) + + # translate char to unsigned char + unsigned char pysam_translate_sequence( char s ) + + # stand-ins for samtools macros + uint32_t * pysam_bam1_cigar( bam1_t * b) + char * pysam_bam1_qname( bam1_t * b) + uint8_t * pysam_bam1_seq( bam1_t * b) + uint8_t * pysam_bam1_qual( bam1_t * b) + uint8_t * pysam_bam1_aux( bam1_t * b) + + # iterator implemenation + ctypedef struct bam_fetch_iterator_t: + pass + + bam_fetch_iterator_t* bam_init_fetch_iterator(bamFile fp, bam_index_t *idx, int tid, int beg, int end) + + bam1_t * bam_fetch_iterate(bam_fetch_iterator_t *iter) + + void bam_cleanup_fetch_iterator(bam_fetch_iterator_t *iter) diff --git a/pysam/csamtools.pyx b/pysam/csamtools.pyx index 826016f..0da8d9e 100644 --- a/pysam/csamtools.pyx +++ b/pysam/csamtools.pyx @@ -1,7 +1,7 @@ # cython: embedsignature=True # adds doc-strings for sphinx -import tempfile, os, sys, types, itertools +import tempfile, os, sys, types, itertools, struct, ctypes # defines imported from samtools DEF SEEK_SET = 0 @@ -52,6 +52,14 @@ cdef makeAlignedRead( bam1_t * src): dest._delegate = bam_dup1(src) return dest +cdef class PileupProxy +cdef makePileupProxy( bam_plbuf_t * buf, int n ): + cdef PileupProxy dest + dest = PileupProxy() + dest.buf = buf + dest.n = n + return dest + cdef class PileupRead cdef makePileupRead( bam_pileup1_t * src ): '''fill a PileupRead object from a bam_pileup1_t * object.''' @@ -79,6 +87,23 @@ cdef int fetch_callback( bam1_t *alignment, void *f): a = makeAlignedRead( alignment ) (f)(a) +class PileupColumn(object): + '''A pileup column. A pileup column contains + all the reads that map to a certain target base. + + tid + chromosome ID as is defined in the header + pos + the target base coordinate (0-based) + n + number of reads mapping to this column + pileups + list of reads (:class:`pysam.PileupRead`) aligned to this column + ''' + def __str__(self): + return "\t".join( map(str, (self.tid, self.pos, self.n))) +\ + "\n" + "\n".join( map(str, self.pileups) ) + cdef int pileup_callback( uint32_t tid, uint32_t pos, int n, bam_pileup1_t *pl, void *f): '''callback for pileup. @@ -118,6 +143,23 @@ cdef int pileup_fetch_callback( bam1_t *b, void *data): bam_plbuf_push(b, buf) return 0 +class StderrStore(): + ''' + stderr is captured. + ''' + def __init__(self): + self.stderr_h, self.stderr_f = tempfile.mkstemp() + self.stderr_save = Outs( sys.stderr.fileno() ) + self.stderr_save.setfd( self.stderr_h ) + + def release(self): + self.stderr_save.restore() + if os.path.exists(self.stderr_f): + os.remove( self.stderr_f ) + + def __del__(self): + self.release() + ###################################################################### ###################################################################### ###################################################################### @@ -125,7 +167,7 @@ cdef int pileup_fetch_callback( bam1_t *b, void *data): VALID_HEADER_TYPES = { "HD" : dict, "SQ" : list, "RG" : list, - "PG" : dict, + "PG" : list, "CO" : list } # order of records within sam headers @@ -152,10 +194,10 @@ cdef class Samfile: '''*(filename, mode='r', template = None, referencenames = None, referencelengths = None, text = NULL, header = None)* A *SAM* file. The file is automatically opened. - - *mode* should be ``r`` for reading or ``w`` for writing. The default is text mode so for binary (:term:`BAM`) I/O you should append - ``b`` for compressed or ``u`` for uncompressed :term:`BAM` output. Use ``h`` to output header information in text (:term:`TAM`) mode. + *mode* should be ``r`` for reading or ``w`` for writing. The default is text mode so for binary + (:term:`BAM`) I/O you should append ``b`` for compressed or ``u`` for uncompressed :term:`BAM` output. + Use ``h`` to output header information in text (:term:`TAM`) mode. If ``b`` is present, it must immediately follow ``r`` or ``w``. Currently valid modes are ``r``, ``w``, ``wh``, ``rb``, ``wb`` and ``wbu``. @@ -176,9 +218,8 @@ cdef class Samfile: 4. The names (*referencenames*) and lengths (*referencelengths*) are supplied directly as lists. - - - + If an index for a BAM file exists (.bai), it will be opened automatically. Without an index random + access to reads via :meth:`fetch` and :meth:`pileup` is disabled. ''' cdef char * filename @@ -189,11 +230,17 @@ cdef class Samfile: # true if file is a bam file cdef int isbam + # current read within iteration + cdef bam1_t * b + def __cinit__(self, *args, **kwargs ): self.samfile = NULL self.isbam = False self._open( *args, **kwargs ) + # allocate memory for iterator + self.b = calloc(1, sizeof(bam1_t)) + def _isOpen( self ): '''return true if samfile has been opened.''' return self.samfile != NULL @@ -203,13 +250,13 @@ cdef class Samfile: return self.index != NULL def _open( self, - char * filename, - mode ='r', - Samfile template = None, - referencenames = None, - referencelengths = None, - char * text = NULL, - header = None, + char * filename, + mode ='r', + Samfile template = None, + referencenames = None, + referencelengths = None, + char * text = NULL, + header = None, ): '''open a sam/bam file. @@ -221,11 +268,11 @@ cdef class Samfile: # close a previously opened file if self.samfile != NULL: self.close() + self.samfile = NULL cdef bam_header_t * header_to_write header_to_write = NULL - if self.samfile != NULL: self.close() self.filename = filename self.isbam = len(mode) > 1 and mode[1] == 'b' @@ -240,9 +287,10 @@ cdef class Samfile: elif header: header_to_write = self._buildHeader( header ) + else: # build header from a target names and lengths - assert referencenames and referencelengths, "supply names and lengths of reference sequences for writing" + assert referencenames and referencelengths, "either supply options `template`, `header` or both `refernencenames` and `referencelengths` for writing" assert len(referencenames) == len(referencelengths), "unequal names and lengths of reference sequences" # allocate and fill header @@ -269,7 +317,9 @@ cdef class Samfile: # open file. Header gets written to file at the same time for bam files # and sam files (in the latter case, the mode needs to be wh) + store = StderrStore() self.samfile = samopen( filename, mode, header_to_write ) + store.release() # bam_header_destroy takes care of cleaning up of all the members if not template and header_to_write != NULL: @@ -277,15 +327,24 @@ cdef class Samfile: elif mode[0] == "r": # open file for reading + if strncmp( filename, "-", 1) != 0 and not os.path.exists( filename ): + raise IOError( "file `%s` not found" % filename) + + store = StderrStore() self.samfile = samopen( filename, mode, NULL ) + store.release() if self.samfile == NULL: raise IOError("could not open file `%s`" % filename ) if mode[0] == "r" and self.isbam: - self.index = bam_index_load(filename) - if self.index == NULL: - raise IOError("could not open index `%s` " % filename ) + if not os.path.exists(filename + ".bai"): + self.index = NULL + else: + # returns NULL if there is no index or index could not be opened + self.index = bam_index_load(filename) + if self.index == NULL: + raise IOError("error while opening index `%s` " % filename ) def getrname( self, tid ): '''(tid ) @@ -294,8 +353,11 @@ cdef class Samfile: raise ValueError( "tid out of range 0<=tid<%i" % self.samfile.header.n_targets ) return self.samfile.header.target_name[tid] - def _parseRegion( self, reference = None, start = None, end = None, - region = None ): + def _parseRegion( self, + reference = None, + start = None, + end = None, + region = None ): '''parse region information. raise Value for for invalid regions. @@ -310,6 +372,8 @@ cdef class Samfile: cdef int rtid cdef int rstart cdef int rend + cdef int max_pos + max_pos = 2 << 29 rtid = rstart = rend = 0 @@ -321,13 +385,14 @@ cdef class Samfile: region = reference if region: + store = StderrStore() bam_parse_region( self.samfile.header, region, &rtid, &rstart, &rend) + store.release() if rtid < 0: raise ValueError( "invalid region `%s`" % region ) - if rstart > rend: raise ValueError( 'invalid region: start (%i) > end (%i)' % (rstart, rend) ) - if rstart < 0: raise ValueError( 'negative start coordinate (%i)' % rstart ) - if rend < 0: raise ValueError( 'negative end coordinate (%i)' % rend ) - + if not 0 <= rstart < max_pos: raise ValueError( 'start out of range (%i)' % rstart ) + if not 0 <= rend < max_pos: raise ValueError( 'end out of range (%i)' % rend ) + return region, rtid, rstart, rend def fetch( self, @@ -361,14 +426,18 @@ cdef class Samfile: cdef int rstart cdef int rend + if not self._isOpen(): + raise ValueError( "I/O operation on closed file" ) + region, rtid, rstart, rend = self._parseRegion( reference, start, end, region ) if self.isbam: - assert self._hasIndex(), "no index available for fetch" if callback: if not region: raise ValueError( "callback functionality requires a region/reference" ) - return bam_fetch(self.samfile.x.bam, self.index, rtid, rstart, rend, callback, fetch_callback ) + if not self._hasIndex(): raise ValueError( "no index available for fetch" ) + return bam_fetch(self.samfile.x.bam, + self.index, rtid, rstart, rend, callback, fetch_callback ) else: if region: return IteratorRow( self, rtid, rstart, rend ) @@ -377,6 +446,7 @@ cdef class Samfile: return IteratorRowAll( self ) else: # return all targets by chaining the individual targets together. + if not self._hasIndex(): raise ValueError( "no index available for fetch" ) i = [] rstart = 0 rend = 1<<29 @@ -406,23 +476,28 @@ cdef class Samfile: an exception is raised. .. Note:: - Note that *all* reads which overlap the region are returned. So the first base returned will be the first base of the first read *not* the first base of the region. + + *all* reads which overlap the region are returned. The first base returned will be the + first base of the first read *not* necessarily the first base of the region used in the query. ''' cdef int rtid cdef int rstart cdef int rend cdef bam_plbuf_t *buf + if not self._isOpen(): + raise ValueError( "I/O operation on closed file" ) + region, rtid, rstart, rend = self._parseRegion( reference, start, end, region ) if self.isbam: - assert self._hasIndex(), "no index available for pileup" + if not self._hasIndex(): raise ValueError( "no index available for pileup" ) if callback: if not region: raise ValueError( "callback functionality requires a region/reference" ) - buf = bam_plbuf_init(pileup_callback, callback ) + buf = bam_plbuf_init( pileup_callback, callback ) bam_fetch(self.samfile.x.bam, self.index, rtid, rstart, rend, buf, pileup_fetch_callback ) @@ -454,8 +529,10 @@ cdef class Samfile: def __dealloc__( self ): '''clean up.''' + # remember: dealloc cannot call other methods # Note that __del__ is not called. self.close() + pysam_bam_destroy1(self.b) def write( self, AlignedRead read ): '''(AlignedRead read ) @@ -523,8 +600,8 @@ cdef class Samfile: # the following is clumsy as generators do not work? x = {} - for field in fields[1:]: - key,value = field.split(":") + for field in fields[1:]: + key, value = field.split(":",1) if key not in VALID_HEADER_FIELDS[record]: raise ValueError( "unknown field code '%s' in record '%s'" % (key, record) ) x[key] = VALID_HEADER_FIELDS[record][key](value) @@ -606,6 +683,114 @@ cdef class Samfile: return dest + def __iter__(self): + return self + + cdef bam1_t * getCurrent( self ): + return self.b + + cdef int cnext(self): + '''cversion of iterator. Used by IteratorColumn''' + cdef int ret + return samread(self.samfile, self.b) + + def __next__(self): + """python version of next(). + + pyrex uses this non-standard name instead of next() + """ + cdef int ret + ret = samread(self.samfile, self.b) + if (ret > 0): + return makeAlignedRead( self.b ) + else: + raise StopIteration + +cdef class Fastafile: + '''*(filename)* + + A *FASTA* file. The file is automatically opened. + + The file expects an indexed fasta file. + + TODO: + add automatic indexing. + add function to get sequence names. + ''' + + cdef char * filename + # pointer to fastafile + cdef faidx_t * fastafile + + def __cinit__(self, *args, **kwargs ): + self.fastafile = NULL + self._open( *args, **kwargs ) + + def _isOpen( self ): + '''return true if samfile has been opened.''' + return self.fastafile != NULL + + def _open( self, + char * filename ): + '''open an indexed fasta file. + + This method expects an indexed fasta file. + ''' + + # close a previously opened file + if self.fastafile != NULL: self.close() + self.filename = filename + self.fastafile = fai_load( filename ) + + if self.fastafile == NULL: + raise IOError("could not open file `%s`" % filename ) + + def close( self ): + if self.fastafile != NULL: + fai_destroy( self.fastafile ) + self.fastafile = NULL + + def fetch( self, + reference = None, + start = None, + end = None, + region = None): + + '''*(reference = None, start = None, end = None, region = None)* + + fetch :meth:`AlignedRead` objects in a :term:`region` using 0-based indexing. The region is specified by + :term:`reference`, *start* and *end*. Alternatively, a samtools :term:`region` string can be supplied. + ''' + + if not self._isOpen(): + raise ValueError( "I/O operation on closed file" ) + + cdef int len, max_pos + cdef char * seq + max_pos = 2 << 29 + + if not region: + if reference == None: raise ValueError( 'no sequence/region supplied.' ) + if start == None and end == None: + region = "%s" % str(reference) + elif start == None or end == None: + raise ValueError( 'only start or only end of region supplied' ) + else: + if start > end: raise ValueError( 'invalid region: start (%i) > end (%i)' % (start, end) ) + # valid ranges are from 0 to 2^29-1 + if not 0 <= start < max_pos: raise ValueError( 'start out of range (%i)' % start ) + if not 0 <= end < max_pos: raise ValueError( 'end out of range (%i)' % end ) + region = "%s:%i-%i" % (reference, start+1, end ) + + # samtools adds a '\0' at the end + seq = fai_fetch( self.fastafile, region, &len ) + # copy to python + result = seq + # clean up + free(seq) + + return result + ## turning callbacks elegantly into iterators is an unsolved problem, see the following threads: ## http://groups.google.com/group/comp.lang.python/browse_frm/thread/0ce55373f128aa4e/1d27a78ca6408134?hl=en&pli=1 ## http://www.velocityreviews.com/forums/t359277-turning-a-callback-function-into-a-generator.html @@ -716,8 +901,26 @@ cdef class IteratorRowAll: cdef class IteratorColumn: '''iterates over columns. - This iterator wraps the pileup functionality of the samtools - function bam_plbuf_push. + This iterator wraps the pileup functionality of samtools. + + For reasons of efficiency, the iterator returns the current + pileup buffer. As this buffer is updated at every iteration, + the contents of this iterator will change accordingly. Hence the conversion to + a list will not produce the expected result:: + + f = Samfile("file.bam", "rb") + result = list( f.pileup() ) + + Here, result will contain ``n`` objects of type :class:`PileupProxy` for ``n`` columns, + but each object will contain the same information. + + If the results of several columns are required at the same time, the results + need to be stored explicitely:: + + result = [ x.pileups() for x in f.pileup() ] + + Here, result will be a list of ``n`` lists of objects of type :class:`PileupRead`. + ''' cdef bam_plbuf_t *buf @@ -786,19 +989,7 @@ cdef class IteratorColumn: cdef bam_pileup1_t * pl if ret > 0 : - p = PileupColumn() - p.tid = pysam_get_tid( self.buf ) - p.pos = pysam_get_pos( self.buf ) - p.n = self.n_pu - pl = pysam_get_pileup( self.buf ) - - pileups = [] - - for x from 0 <= x < p.n: - pileups.append( makePileupRead( &pl[x]) ) - p.pileups = pileups - - return p + return makePileupProxy( self.buf, self.n_pu ) else: raise StopIteration @@ -808,16 +999,38 @@ cdef class IteratorColumn: cdef class AlignedRead: ''' Class representing an aligned read. see SAM format specification for meaning of fields (http://samtools.sourceforge.net/). - + + This class stores a handle to the samtools C-structure representing + an aligned read. Member read access is forwarded to the C-structure + and converted into python objects. This implementation should be fast, + as only the data needed is converted. + + For write access, the C-structure is updated in-place. This is + not the most efficient way to build BAM entries, as the variable + length data is concatenated and thus needs to resized if + a field is updated. Furthermore, the BAM entry might be + in an inconsistent state. The :meth:`~validate` method can + be used to check if an entry is consistent. + + One issue to look out for is that the sequence should always + be set *before* the quality scores. Setting the sequence will + also erase any quality scores that were set previously. ''' cdef: bam1_t * _delegate def __cinit__( self ): - self._delegate = calloc( sizeof( bam1_t), 1 ) + # see bam_init1 + self._delegate = calloc( 1, sizeof( bam1_t) ) + # allocate some memory + # If size is 0, calloc does not return a pointer that can be passed to free() + # so allocate 40 bytes for a new read + self._delegate.m_data = 40 + self._delegate.data = calloc( self._delegate.m_data, 1 ) + self._delegate.data_len = 0 def __dealloc__(self): - """todo is this enough or do we need to free() each string? eg 'qual' etc""" + '''clear up memory.''' pysam_bam_destroy1(self._delegate) def __str__(self): @@ -825,202 +1038,537 @@ cdef class AlignedRead: return "\t".join(map(str, (self.qname, self.rname, self.pos, + self.cigar, self.qual, self.flag, self.seq, - self.mapq ))) + self.mapq, + self.tags))) + + def __cmp__(self, AlignedRead other): + '''return true, if contents in this are binary equal to ``other``.''' + cdef int retval, x + cdef bam1_t *t, *o + t = self._delegate + o = other._delegate + + # uncomment for debugging purposes + # cdef unsigned char * oo, * tt + # tt = (&t.core) + # oo = (&o.core) + # for x from 0 <= x < sizeof( bam1_core_t): print x, tt[x], oo[x] + # tt = (t.data) + # oo = (o.data) + # for x from 0 <= x < max(t.data_len, o.data_len): print x, tt[x], oo[x], chr(tt[x]), chr(oo[x]) + + retval = memcmp( &t.core, + &o.core, + sizeof( bam1_core_t )) + + if retval: return retval + retval = cmp( t.data_len, o.data_len) + if retval: return retval + return memcmp( t.data, + o.data, + sizeof( t.data_len )) + + property qname: + """the query name (None if not present)""" + def __get__(self): + cdef bam1_t * src + src = self._delegate + if src.core.l_qname == 0: return None + return pysam_bam1_qname( src ) + + def __set__(self, qname ): + if qname == None or len(qname) == 0: return + cdef bam1_t * src + cdef int l + cdef char * p + + src = self._delegate + p = pysam_bam1_qname( src ) + + # the qname is \0 terminated + l = len(qname) + 1 + pysam_bam_update( src, + src.core.l_qname, + l, + p ) + + src.core.l_qname = l + + # re-acquire pointer to location in memory + # as it might have moved + p = pysam_bam1_qname(src) + + strncpy( p, qname, l ) + property cigar: - """the :term:`cigar` alignment string (None if not present)""" + """the :term:`cigar` alignment (None if not present). + """ def __get__(self): cdef uint32_t * cigar_p cdef bam1_t * src cdef op, l, cigar src = self._delegate - if src.core.n_cigar > 0: - cigar = [] - cigar_p = < uint32_t *> (src.data + src.core.l_qname) - for k from 0 <= k < src.core.n_cigar: - op = cigar_p[k] & BAM_CIGAR_MASK - l = cigar_p[k] >> BAM_CIGAR_SHIFT - cigar.append((op, l)) - return cigar - return None - property qname: - """the query name (None if not present)""" - def __get__(self): + if src.core.n_cigar == 0: return None + + cigar = [] + cigar_p = pysam_bam1_cigar(src); + for k from 0 <= k < src.core.n_cigar: + op = cigar_p[k] & BAM_CIGAR_MASK + l = cigar_p[k] >> BAM_CIGAR_SHIFT + cigar.append((op, l)) + return cigar + + def __set__(self, values ): + if values == None or len(values) == 0: return + cdef uint32_t * p cdef bam1_t * src + cdef op, l + cdef int k + + k = 0 + src = self._delegate - ## parse qname (bam1_qname) - return < char *> src.data + + # get location of cigar string + p = pysam_bam1_cigar(src) + + # create space for cigar data within src.data + pysam_bam_update( src, + src.core.n_cigar * 4, + len(values) * 4, + p ) + + # length is number of cigar operations, not bytes + src.core.n_cigar = len(values) + + # re-acquire pointer to location in memory + # as it might have moved + p = pysam_bam1_cigar(src) + + # insert cigar operations + for op, l in values: + p[k] = l << BAM_CIGAR_SHIFT | op + k += 1 + + ## setting the cigar string also updates the "bin" attribute + src.core.bin = bam_reg2bin( src.core.pos, bam_calend( &src.core, p)) + property seq: """the query sequence (None if not present)""" def __get__(self): cdef bam1_t * src - cdef uint8_t * seq_p + cdef uint8_t * p cdef char * s src = self._delegate bam_nt16_rev_table = "=ACMGRSVTWYHKDBN" ## parse qseq (bam1_seq) - if src.core.l_qseq: - s = < char *> calloc(src.core.l_qseq + 1 , sizeof(char)) - seq_p = < uint8_t *> (src.data + src.core.n_cigar * 4 + src.core.l_qname) - for k from 0 <= k < src.core.l_qseq: - ## equivalent to bam_nt16_rev_table[bam1_seqi(s, i)] (see bam.c) - s[k] = "=ACMGRSVTWYHKDBN"[((seq_p)[(k) / 2] >> 4 * (1 - (k) % 2) & 0xf)] - retval=s - free(s) - return retval - return None + if src.core.l_qseq == 0: return None + + s = < char *> calloc(src.core.l_qseq + 1 , sizeof(char)) + p = pysam_bam1_seq( src ) + for k from 0 <= k < src.core.l_qseq: + ## equivalent to bam_nt16_rev_table[bam1_seqi(s, i)] (see bam.c) + s[k] = "=ACMGRSVTWYHKDBN"[((p)[(k) / 2] >> 4 * (1 - (k) % 2) & 0xf)] + retval=s + free(s) + return retval + + def __set__(self,seq): + # samtools manages sequence and quality length memory together + # if no quality information is present, the first byte says 0xff. + + if seq == None or len(seq) == 0: return + cdef bam1_t * src + cdef uint8_t * p + cdef char * s + src = self._delegate + cdef int l, k, nbytes_new, nbytes_old + + l = len(seq) + + # as the sequence is stored in half-bytes, the total length (sequence + # plus quality scores) is (l+1)/2 + l + nbytes_new = (l+1)/2 + l + nbytes_old = (src.core.l_qseq+1)/2 + src.core.l_qseq + # acquire pointer to location in memory + p = pysam_bam1_seq( src ) + src.core.l_qseq = l + + pysam_bam_update( src, + nbytes_old, + nbytes_new, + p) + # re-acquire pointer to location in memory + # as it might have moved + p = pysam_bam1_seq( src ) + for k from 0 <= k < nbytes_new: p[k] = 0 + # convert to C string + s = seq + for k from 0 <= k < l: + p[k/2] |= pysam_translate_sequence(s[k]) << 4 * (1 - k % 2) + + # erase qualities + p = pysam_bam1_qual( src ) + p[0] = 0xff + property qual: """the base quality (None if not present)""" def __get__(self): - #todo is this still needed cdef bam1_t * src - cdef uint8_t * qual_p + cdef uint8_t * p cdef char * q src = self._delegate - qual_p = < uint8_t *> (src.data + src.core.n_cigar * 4 + src.core.l_qname + (src.core.l_qseq + 1) / 2) - if qual_p[0] != 0xff: - q = < char *> calloc(src.core.l_qseq + 1 , sizeof(char)) - for k from 0 <= k < src.core.l_qseq: - ## equivalent to t[i] + 33 (see bam.c) - q[k] = qual_p[k] + 33 - retval=q - free(q) - return retval - return None + if src.core.l_qseq == 0: return None + + p = pysam_bam1_qual( src ) + if p[0] == 0xff: return None + + q = < char *>calloc(src.core.l_qseq + 1 , sizeof(char)) + for k from 0 <= k < src.core.l_qseq: + ## equivalent to t[i] + 33 (see bam.c) + q[k] = p[k] + 33 + # convert to python string + retval=q + # clean up + free(q) + return retval + + def __set__(self,qual): + # note that space is already allocated via the sequences + cdef bam1_t * src + cdef uint8_t * p + cdef char * q + src = self._delegate + p = pysam_bam1_qual( src ) + if qual == None or len(qual) == 0: + # if absent - set to 0xff + p[0] = 0xff + return + cdef int l + # convert to C string + q = qual + l = len(qual) + if src.core.l_qseq != l: + raise ValueError("quality and sequence mismatch: %i != %i" % (l, src.core.l_qseq)) + assert src.core.l_qseq == l + for k from 0 <= k < l: + p[k] = q[k] - 33 + + property tags: + """the tags in the AUX field.""" + def __get__(self): + cdef char * ctag + cdef bam1_t * src + cdef uint8_t * s + cdef char tpe + + src = self._delegate + if src.l_aux == 0: return None + + s = pysam_bam1_aux( src ) + result = [] + ctag = calloc( 3, sizeof(char) ) + cdef int x + while s < (src.data + src.data_len): + # get tag + ctag[0] = s[0] + ctag[1] = s[1] + pytag = ctag + + s += 2 + + # convert type - is there a better way? + ctag[0] = s[0] + ctag[1] = 0 + pytype = ctag + # get type and value + # how do I do char literal comparison in cython? + # the code below works (i.e, is C comparison) + tpe = toupper(s[0]) + if tpe == 'S'[0]: + value = bam_aux2i(s) + s += 2 + elif tpe == 'I'[0]: + value = bam_aux2i(s) + s += 4 + elif tpe == 'F'[0]: + value = bam_aux2f(s) + s += 4 + elif tpe == 'D'[0]: + value = bam_aux2d(s) + s += 8 + elif tpe == 'C'[0]: + value = bam_aux2i(s) + s += 1 + elif tpe == 'A'[0]: + # there might a more efficient way + # to convert a char into a string + value = "%c" % bam_aux2A(s) + s += 1 + elif tpe == 'Z'[0]: + value = bam_aux2Z(s) + # +1 for NULL terminated string + s += len(value) + 1 + + # skip over type + s += 1 + + # ignore pytype + result.append( (pytag, value) ) + + free( ctag ) + return result + + def __set__(self, tags): + cdef char * ctag + cdef bam1_t * src + cdef uint8_t * s + cdef uint8_t * new_data + cdef int guessed_size, control_size + src = self._delegate + cdef int max_size, size + max_size = 4000 + + # map samtools code to python.struct code and byte size + buffer = ctypes.create_string_buffer(max_size) + + offset = 0 + for pytag, value in tags: + t = type(value) + if t == types.FloatType: + fmt = "= -127: fmt, pytype = "= -32767: fmt, pytype = " 4294967295: raise ValueError( "integer %i out of range of BAM/SAM specification" % value ) + else: fmt, pytype = " max_size: + raise NotImplementedError("tags field too large") + + struct.pack_into( fmt, + buffer, + offset, + pytag[0], + pytag[1], + pytype, + value ) + offset += size + + # delete the old data and allocate new + pysam_bam_update( src, + src.l_aux, + offset, + pysam_bam1_aux( src ) ) + + src.l_aux = offset + + if offset == 0: return + + # get location of new data + s = pysam_bam1_aux( src ) + + # check if there is direct path from buffer.raw to tmp + cdef char * temp + temp = buffer.raw + memcpy( s, temp, offset ) property flag: """properties flag""" - def __get__(self): - return self._delegate.core.flag + def __get__(self): return self._delegate.core.flag + def __set__(self, flag): self._delegate.core.flag = flag property rname: - """:term:`reference` ID""" - def __get__(self): - return self._delegate.core.tid + """ + :term:`target` ID + + .. note:: + + This field contains the index of the reference sequence + in the sequence dictionary. To obtain the name + of the reference sequence, use :meth:`pysam.Samfile.getrname()` + + """ + def __get__(self): return self._delegate.core.tid + def __set__(self, tid): self._delegate.core.tid = tid property pos: """0-based leftmost coordinate""" - def __get__(self): - return self._delegate.core.pos + def __get__(self): return self._delegate.core.pos + def __set__(self, pos): + ## setting the cigar string also updates the "bin" attribute + cdef bam1_t * src + src = self._delegate + if src.core.n_cigar: + src.core.bin = bam_reg2bin( src.core.pos, bam_calend( &src.core, pysam_bam1_cigar(src)) ) + else: + src.core.bin = bam_reg2bin( src.core.pos, src.core.pos + 1) + self._delegate.core.pos = pos + property bin: + """properties bin""" + def __get__(self): return self._delegate.core.bin + def __set__(self, bin): self._delegate.core.bin = bin property rlen: - '''length of the read. Returns 0 if not given.''' - def __get__(self): - return self._delegate.core.l_qseq + '''length of the read (read only). Returns 0 if not given.''' + def __get__(self): return self._delegate.core.l_qseq property mapq: """mapping quality""" - def __get__(self): - return self._delegate.core.qual + def __get__(self): return self._delegate.core.qual + def __set__(self, qual): self._delegate.core.qual = qual property mrnm: """the :term:`reference` id of the mate """ - def __get__(self): - return self._delegate.core.mtid + def __get__(self): return self._delegate.core.mtid + def __set__(self, mtid): self._delegate.core.mtid = mtid property mpos: """the position of the mate""" - def __get__(self): - return self._delegate.core.mpos + def __get__(self): return self._delegate.core.mpos + def __set__(self, mpos): self._delegate.core.mpos = mpos property isize: """the insert size""" - def __get__(self): - return self._delegate.core.isize + def __get__(self): return self._delegate.core.isize + def __set__(self, isize): self._delegate.core.isize = isize property is_paired: """true if read is paired in sequencing""" def __get__(self): return (self._delegate.core.flag & BAM_FPAIRED) != 0 + def __set__(self,val): + if val: self._delegate.core.flag |= BAM_FPAIRED + else: self._delegate.core.flag &= ~BAM_FPAIRED property is_proper_pair: """true if read is mapped in a proper pair""" def __get__(self): return (self.flag & BAM_FPROPER_PAIR) != 0 + def __set__(self,val): + if val: self._delegate.core.flag |= BAM_FPROPER_PAIR + else: self._delegate.core.flag &= ~BAM_FPROPER_PAIR property is_unmapped: """true if read itself is unmapped""" def __get__(self): return (self.flag & BAM_FUNMAP) != 0 + def __set__(self,val): + if val: self._delegate.core.flag |= BAM_FUNMAP + else: self._delegate.core.flag &= ~BAM_FUNMAP property mate_is_unmapped: """true if the mate is unmapped""" def __get__(self): return (self.flag & BAM_FMUNMAP) != 0 + def __set__(self,val): + if val: self._delegate.core.flag |= BAM_FMUNMAP + else: self._delegate.core.flag &= ~BAM_FMUNMAP property is_reverse: """true if read is mapped to reverse strand""" def __get__(self):return (self.flag & BAM_FREVERSE) != 0 + def __set__(self,val): + if val: self._delegate.core.flag |= BAM_FREVERSE + else: self._delegate.core.flag &= ~BAM_FREVERSE property mate_is_reverse: """true is read is mapped to reverse strand""" def __get__(self): return (self.flag & BAM_FMREVERSE) != 0 + def __set__(self,val): + if val: self._delegate.core.flag |= BAM_FMREVERSE + else: self._delegate.core.flag &= ~BAM_FMREVERSE property is_read1: """true if this is read1""" def __get__(self): return (self.flag & BAM_FREAD1) != 0 + def __set__(self,val): + if val: self._delegate.core.flag |= BAM_FREAD1 + else: self._delegate.core.flag &= ~BAM_FREAD1 property is_read2: """true if this is read2""" def __get__(self): return (self.flag & BAM_FREAD2) != 0 + def __set__(self,val): + if val: self._delegate.core.flag |= BAM_FREAD2 + else: self._delegate.core.flag &= ~BAM_FREAD2 property is_secondary: """true if not primary alignment""" def __get__(self): return (self.flag & BAM_FSECONDARY) != 0 + def __set__(self,val): + if val: self._delegate.core.flag |= BAM_FSECONDARY + else: self._delegate.core.flag &= ~BAM_FSECONDARY property is_qcfail: """true if QC failure""" - def __get__(self): return (self.flag & BAM_FSECONDARY) != 0 + def __get__(self): return (self.flag & BAM_FQCFAIL) != 0 + def __set__(self,val): + if val: self._delegate.core.flag |= BAM_FQCFAIL + else: self._delegate.core.flag &= ~BAM_FQCFAIL property is_duplicate: """ true if optical or PCR duplicate""" def __get__(self): return (self.flag & BAM_FDUP) != 0 + def __set__(self,val): + if val: self._delegate.core.flag |= BAM_FDUP + else: self._delegate.core.flag &= ~BAM_FDUP def opt(self, tag): """retrieves optional data given a two-letter *tag*""" #see bam_aux.c: bam_aux_get() and bam_aux2i() etc cdef uint8_t * v v = bam_aux_get(self._delegate, tag) + if v == NULL: raise KeyError( "tag '%s' not present" % tag ) type = chr(v[0]) if type == 'c' or type == 'C' or type == 's' or type == 'S' or type == 'i': - return bam_aux2i(v) + return bam_aux2i(v) elif type == 'f': - return bam_aux2f(v) + return bam_aux2f(v) elif type == 'd': - return bam_aux2d(v) + return bam_aux2d(v) elif type == 'A': - return bam_aux2A(v) + # there might a more efficient way + # to convert a char into a string + return '%c' % bam_aux2A(v) elif type == 'Z': - return bam_aux2Z(v) + return bam_aux2Z(v) - # def fancy_str (self): - # """ - # returns list of fieldnames/values in pretty format for debugging - # """ - # ret_string = [] - # field_names = { - # "tid": "Contig index", - # "pos": "Mapped position on contig", - # - # "mtid": "Contig index for mate pair", - # "mpos": "Position of mate pair", - # "isize": "Insert size", - # - # "flag": "Binary flag", - # "n_cigar": "Count of cigar entries", - # "cigar": "Cigar entries", - # "qual": "Mapping quality", - # - # "bin": "Bam index bin number", - # - # "l_qname": "Length of query name", - # "qname": "Query name", - # - # "l_qseq": "Length of query sequence", - # "qseq": "Query sequence", - # "bqual": "Quality scores", - # - # - # "l_aux": "Length of auxilary data", - # "m_data": "Maximum data length", - # "data_len": "Current data length", - # } - # fields_names_in_order = ["tid", "pos", "mtid", "mpos", "isize", "flag", - # "n_cigar", "cigar", "qual", "bin", "l_qname", "qname", - # "l_qseq", "qseq", "bqual", "l_aux", "m_data", "data_len"] - # - # for f in fields_names_in_order: - # if not f in self.__dict__: - # continue - # ret_string.append("%-30s %-10s= %s" % (field_names[f], "(" + f + ")", self.__getattribute__(f))) - # - # for f in self.__dict__: - # if not f in field_names: - # ret_string.append("%-30s %-10s= %s" % (f, "", self.__getattribute__(f))) -# return ret_string - -class PileupColumn(object): + def fancy_str (self): + """returns list of fieldnames/values in pretty format for debugging + """ + ret_string = [] + field_names = { + "tid": "Contig index", + "pos": "Mapped position on contig", + "mtid": "Contig index for mate pair", + "mpos": "Position of mate pair", + "isize": "Insert size", + "flag": "Binary flag", + "n_cigar": "Count of cigar entries", + "cigar": "Cigar entries", + "qual": "Mapping quality", + "bin": "Bam index bin number", + "l_qname": "Length of query name", + "qname": "Query name", + "l_qseq": "Length of query sequence", + "qseq": "Query sequence", + "bqual": "Quality scores", + "l_aux": "Length of auxilary data", + "m_data": "Maximum data length", + "data_len": "Current data length", + } + fields_names_in_order = ["tid", "pos", "mtid", "mpos", "isize", "flag", + "n_cigar", "cigar", "qual", "bin", "l_qname", "qname", + "l_qseq", "qseq", "bqual", "l_aux", "m_data", "data_len"] + + for f in fields_names_in_order: + if not f in self.__dict__: + continue + ret_string.append("%-30s %-10s= %s" % (field_names[f], "(" + f + ")", self.__getattribute__(f))) + + for f in self.__dict__: + if not f in field_names: + ret_string.append("%-30s %-10s= %s" % (f, "", self.__getattribute__(f))) + return ret_string + +cdef class PileupProxy: '''A pileup column. A pileup column contains all the reads that map to a certain target base. @@ -1029,15 +1577,49 @@ class PileupColumn(object): pos the target base coordinate (0-based) n - number of reads mapping to this column + number of reads mapping to this column pileups list of reads (:class:`pysam.PileupRead`) aligned to this column + + This class is a proxy for results returned by the samtools pileup engine. + If the underlying engine iterator advances, the results of this column + will change. ''' + cdef bam_plbuf_t * buf + cdef int n_pu + + def __cinit__(self ): + pass + def __str__(self): return "\t".join( map(str, (self.tid, self.pos, self.n))) +\ "\n" +\ "\n".join( map(str, self.pileups) ) + property tid: + '''the chromosome ID as is defined in the header''' + def __get__(self): return pysam_get_tid( self.buf ) + + property n: + '''number of reads mapping to this column.''' + def __get__(self): return self.n_pu + def __set__(self, n): self.n_pu = n + + property pos: + def __get__(self): return pysam_get_pos( self.buf ) + + property pileups: + '''list of reads (:class:`pysam.PileupRead`) aligned to this column''' + def __get__(self): + cdef bam_pileup1_t * pl + pl = pysam_get_pileup( self.buf ) + pileups = [] + # warning: there could be problems if self.n and self.buf are + # out of sync. + for x from 0 <= x < self.n_pu: + pileups.append( makePileupRead( &pl[x]) ) + return pileups + cdef class PileupRead: '''A read aligned to a column. ''' @@ -1083,7 +1665,6 @@ cdef class PileupRead: def __get__(self): return self._level - class Outs: '''http://mail.python.org/pipermail/python-list/2000-June/038406.html''' def __init__(self, id = 1): @@ -1163,7 +1744,6 @@ def _samtools_dispatch( method, args = () ): cdef char ** cargs cdef int i, n, retval - n = len(args) # allocate two more for first (dummy) argument (contains command) cargs = calloc( n+2, sizeof( char *) ) @@ -1188,7 +1768,15 @@ def _samtools_dispatch( method, args = () ): return retval, out_stderr, out_stdout -__all__ = ["Samfile", "IteratorRow", "IteratorRowAll", "IteratorColumn", "AlignedRead", "PileupColumn", "PileupRead" ] +__all__ = ["Samfile", + "Fastafile", + "IteratorRow", + "IteratorRowAll", + "IteratorColumn", + "AlignedRead", + "PileupColumn", + "PileupProxy", + "PileupRead" ] diff --git a/pysam/namedtuple.py b/pysam/namedtuple.py new file mode 100644 index 0000000..a60fb1a --- /dev/null +++ b/pysam/namedtuple.py @@ -0,0 +1,117 @@ +from operator import itemgetter as _itemgetter +from keyword import iskeyword as _iskeyword +import sys as _sys + +def namedtuple(typename, field_names, verbose=False, rename=False): + """Returns a new subclass of tuple with named fields. + + >>> Point = namedtuple('Point', 'x y') + >>> Point.__doc__ # docstring for the new class + 'Point(x, y)' + >>> p = Point(11, y=22) # instantiate with positional args or keywords + >>> p[0] + p[1] # indexable like a plain tuple + 33 + >>> x, y = p # unpack like a regular tuple + >>> x, y + (11, 22) + >>> p.x + p.y # fields also accessable by name + 33 + >>> d = p._asdict() # convert to a dictionary + >>> d['x'] + 11 + >>> Point(**d) # convert from a dictionary + Point(x=11, y=22) + >>> p._replace(x=100) # _replace() is like str.replace() but targets named fields + Point(x=100, y=22) + + """ + + # Parse and validate the field names. Validation serves two purposes, + # generating informative error messages and preventing template injection attacks. + if isinstance(field_names, basestring): + field_names = field_names.replace(',', ' ').split() # names separated by whitespace and/or commas + field_names = tuple(map(str, field_names)) + if rename: + names = list(field_names) + seen = set() + for i, name in enumerate(names): + if (not min(c.isalnum() or c=='_' for c in name) or _iskeyword(name) + or not name or name[0].isdigit() or name.startswith('_') + or name in seen): + names[i] = '_%d' % i + seen.add(name) + field_names = tuple(names) + for name in (typename,) + field_names: + if not min(c.isalnum() or c=='_' for c in name): + raise ValueError('Type names and field names can only contain alphanumeric characters and underscores: %r' % name) + if _iskeyword(name): + raise ValueError('Type names and field names cannot be a keyword: %r' % name) + if name[0].isdigit(): + raise ValueError('Type names and field names cannot start with a number: %r' % name) + seen_names = set() + for name in field_names: + if name.startswith('_') and not rename: + raise ValueError('Field names cannot start with an underscore: %r' % name) + if name in seen_names: + raise ValueError('Encountered duplicate field name: %r' % name) + seen_names.add(name) + + # Create and fill-in the class template + numfields = len(field_names) + argtxt = repr(field_names).replace("'", "")[1:-1] # tuple repr without parens or quotes + reprtxt = ', '.join('%s=%%r' % name for name in field_names) + template = '''class %(typename)s(tuple): + '%(typename)s(%(argtxt)s)' \n + __slots__ = () \n + _fields = %(field_names)r \n + def __new__(_cls, %(argtxt)s): + return _tuple.__new__(_cls, (%(argtxt)s)) \n + @classmethod + def _make(cls, iterable, new=tuple.__new__, len=len): + 'Make a new %(typename)s object from a sequence or iterable' + result = new(cls, iterable) + if len(result) != %(numfields)d: + raise TypeError('Expected %(numfields)d arguments, got %%d' %% len(result)) + return result \n + def __repr__(self): + return '%(typename)s(%(reprtxt)s)' %% self \n + def _asdict(self): + 'Return a new dict which maps field names to their values' + return dict(zip(self._fields, self)) \n + def _replace(_self, **kwds): + 'Return a new %(typename)s object replacing specified fields with new values' + result = _self._make(map(kwds.pop, %(field_names)r, _self)) + if kwds: + raise ValueError('Got unexpected field names: %%r' %% kwds.keys()) + return result \n + def __getnewargs__(self): + return tuple(self) \n\n''' % locals() + for i, name in enumerate(field_names): + template += ' %s = _property(_itemgetter(%d))\n' % (name, i) + if verbose: + print template + + # Execute the template string in a temporary namespace + namespace = dict(_itemgetter=_itemgetter, __name__='namedtuple_%s' % typename, + _property=property, _tuple=tuple) + try: + exec template in namespace + except SyntaxError, e: + raise SyntaxError(e.message + ':\n' + template) + result = namespace[typename] + + # For pickling to work, the __module__ variable needs to be set to the frame + # where the named tuple is created. Bypass this step in enviroments where + # sys._getframe is not defined (Jython for example) or sys._getframe is not + # defined for arguments greater than 0 (IronPython). + try: + result.__module__ = _sys._getframe(1).f_globals.get('__name__', '__main__') + except (AttributeError, ValueError): + pass + + return result + + + + + diff --git a/pysam/pysam_util.c b/pysam/pysam_util.c index 02d9db3..5360626 100644 --- a/pysam/pysam_util.c +++ b/pysam/pysam_util.c @@ -13,6 +13,11 @@ // The order of the following declarations is important. // ####################################################### +typedef struct +{ + uint64_t u, v; +} pair64_t; + #define pair64_lt(a,b) ((a).u < (b).u) typedef struct { @@ -138,6 +143,9 @@ static inline int resolve_cigar(bam_pileup1_t *p, uint32_t pos) return ret; } + + + // the following code has been taken from bam_plbuf_push // and modified such that instead of a function call // the function returns and will continue (if cont is true). @@ -244,7 +252,20 @@ bam_pileup1_t * pysam_get_pileup( const bam_plbuf_t *buf) // pysam dispatch function to emulate the samtools // command line within python. // taken from the main function in bamtk.c -// add code to reset getopt +// added code to reset getopt +extern int main_samview(int argc, char *argv[]); +extern int main_import(int argc, char *argv[]); +extern int bam_pileup(int argc, char *argv[]); +extern int bam_merge(int argc, char *argv[]); +extern int bam_sort(int argc, char *argv[]); +extern int bam_index(int argc, char *argv[]); +extern int faidx_main(int argc, char *argv[]); +extern int bam_mating(int argc, char *argv[]); +extern int bam_rmdup(int argc, char *argv[]); +extern int glf3_view_main(int argc, char *argv[]); +extern int bam_flagstat(int argc, char *argv[]); +extern int bam_fillmd(int argc, char *argv[]); + int pysam_dispatch(int argc, char *argv[] ) { @@ -256,6 +277,8 @@ int pysam_dispatch(int argc, char *argv[] ) #endif #endif + extern int optind; + // reset getop optind = 1; @@ -270,10 +293,8 @@ int pysam_dispatch(int argc, char *argv[] ) else if (strcmp(argv[1], "faidx") == 0) return faidx_main(argc-1, argv+1); else if (strcmp(argv[1], "fixmate") == 0) return bam_mating(argc-1, argv+1); else if (strcmp(argv[1], "rmdup") == 0) return bam_rmdup(argc-1, argv+1); - else if (strcmp(argv[1], "rmdupse") == 0) return bam_rmdupse(argc-1, argv+1); else if (strcmp(argv[1], "glfview") == 0) return glf3_view_main(argc-1, argv+1); else if (strcmp(argv[1], "flagstat") == 0) return bam_flagstat(argc-1, argv+1); - // else if (strcmp(argv[1], "tagview") == 0) return bam_tagview(argc-1, argv+1); else if (strcmp(argv[1], "calmd") == 0) return bam_fillmd(argc-1, argv+1); else if (strcmp(argv[1], "fillmd") == 0) return bam_fillmd(argc-1, argv+1); @@ -289,9 +310,190 @@ int pysam_dispatch(int argc, char *argv[] ) } // standin for bam_destroy1 in bam.h +// deletes all variable length data void pysam_bam_destroy1( bam1_t * b ) { - free((b)->data); + if (b == NULL) return; + if (b->data != NULL) free(b->data); free(b); } +// taken from samtools/bam_import.c +static inline uint8_t *alloc_data(bam1_t *b, size_t size) +{ + if (b->m_data < size) + { + b->m_data = size; + kroundup32(b->m_data); + b->data = (uint8_t*)realloc(b->data, b->m_data); + } + return b->data; +} + +// update the variable length data within a bam1_t entry. +// Adds *nbytes_new* - *nbytes_old* into the variable length data of *src* at *pos*. +// Data within the bam1_t entry is moved so that it is +// consistent with the data field lengths. +bam1_t * pysam_bam_update( bam1_t * b, + const size_t nbytes_old, + const size_t nbytes_new, + uint8_t * pos ) +{ + int d = nbytes_new-nbytes_old; + + // no change + if (d == 0) return b; + + int new_size = d + b->data_len; + size_t offset = pos - b->data; + + //printf("d=%i, old=%i, new=%i, old_size=%i, new_size=%i\n", + // d, nbytes_old, nbytes_new, b->data_len, new_size); + + // increase memory if required + if (d > 0) + { + alloc_data( b, new_size ); + pos = b->data + offset; + } + + if (b->data_len != 0) + { + if (offset < 0 || offset > b->data_len) + fprintf(stderr, "[pysam_bam_insert] illegal offset: '%i'\n", (int)offset); + } + + // printf("dest=%p, src=%p, n=%i\n", pos+nbytes_new, pos + nbytes_old, b->data_len - (offset+nbytes_old)); + memmove( pos + nbytes_new, + pos + nbytes_old, + b->data_len - (offset + nbytes_old)); + + b->data_len = new_size; + + return b; +} + +// translate a nucleotide character to binary code +unsigned char pysam_translate_sequence( const unsigned char s ) +{ + return bam_nt16_table[s]; +} + +// stand-ins for samtools macros in bam.h +char * pysam_bam1_qname( const bam1_t * b) +{ + return (char*)b->data; +} + +uint32_t * pysam_bam1_cigar( const bam1_t * b) +{ + return (uint32_t*)(b->data + b->core.l_qname); +} + +uint8_t * pysam_bam1_seq( const bam1_t * b) +{ + return (uint8_t*)(b->data + b->core.n_cigar*4 + b->core.l_qname); +} + +uint8_t * pysam_bam1_qual( const bam1_t * b) +{ + return (uint8_t*)(b->data + b->core.n_cigar*4 + b->core.l_qname + (b->core.l_qseq + 1)/2); +} + +uint8_t * pysam_bam1_aux( const bam1_t * b) +{ + return (uint8_t*)(b->data + b->core.n_cigar*4 + b->core.l_qname + b->core.l_qseq + (b->core.l_qseq + 1)/2); +} + +// ####################################################### +// Iterator implementation +// ####################################################### + +// functions defined in bam_index.c +extern pair64_t * get_chunk_coordinates(const bam_index_t *idx, int tid, int beg, int end, int* cnt_off); + +static inline int is_overlap(uint32_t beg, uint32_t end, const bam1_t *b) +{ + uint32_t rbeg = b->core.pos; + uint32_t rend = b->core.n_cigar? bam_calend(&b->core, bam1_cigar(b)) : b->core.pos + 1; + return (rend > beg && rbeg < end); +} + +struct __bam_fetch_iterator_t +{ + bam1_t * b; + pair64_t * off; + int n_off; + uint64_t curr_off; + int curr_chunk; + bamFile fp; + int tid; + int beg; + int end; + int n_seeks; +}; + +bam_fetch_iterator_t* bam_init_fetch_iterator(bamFile fp, const bam_index_t *idx, int tid, int beg, int end) +{ + // iterator contains current alignment position + // and will contain actual alignment during iterations + bam_fetch_iterator_t* iter = (bam_fetch_iterator_t*)calloc(1, sizeof(bam_fetch_iterator_t)); + iter->b = (bam1_t*)calloc(1, sizeof(bam1_t)); + + // list of chunks containing our alignments + iter->off = get_chunk_coordinates(idx, tid, beg, end, &iter->n_off); + + // initialise other state variables in iterator + iter->fp = fp; + iter->curr_chunk = -1; + iter->curr_off = 0; + iter->n_seeks = 0; + iter->tid = tid; + iter->beg = beg; + iter->end = end; + return iter; +} + +bam1_t * bam_fetch_iterate(bam_fetch_iterator_t *iter) +{ + if (!iter->off) { + return 0; + } + + int ret; + // iterate through all alignments in chunks + for (;;) { + if (iter->curr_off == 0 || iter->curr_off >= iter->off[iter->curr_chunk].v) { // then jump to the next chunk + if (iter->curr_chunk == iter->n_off - 1) break; // no more chunks + if (iter->curr_chunk >= 0) assert(iter->curr_off == iter->off[iter->curr_chunk].v); // otherwise bug + if (iter->curr_chunk < 0 || iter->off[iter->curr_chunk].v != iter->off[iter->curr_chunk+1].u) { // not adjacent chunks; then seek + bam_seek(iter->fp, iter->off[iter->curr_chunk+1].u, SEEK_SET); + iter->curr_off = bam_tell(iter->fp); + ++iter->n_seeks; + } + ++iter->curr_chunk; + } + if ((ret = bam_read1(iter->fp, iter->b)) > 0) { + iter->curr_off = bam_tell(iter->fp); + if (iter->b->core.tid != iter->tid || iter->b->core.pos >= iter->end) break; // no need to proceed + else if (is_overlap(iter->beg, iter->end, iter->b)) + // + //func(iter->b, data); + // + return iter->b; + } else + return 0; // end of file + } + return 0; +} + +void bam_cleanup_fetch_iterator(bam_fetch_iterator_t *iter) +{ + // fprintf(stderr, "[bam_fetch] # seek calls: %d\n", iter->n_seeks); + bam_destroy1(iter->b); + free(iter->off); +} + + + + diff --git a/pysam/pysam_util.h b/pysam/pysam_util.h index 5334d70..ff5d569 100644 --- a/pysam/pysam_util.h +++ b/pysam/pysam_util.h @@ -1,33 +1,54 @@ #ifndef PYSAM_UTIL_H #define PYSAM_UTIL_H -// This file contains some of the definitions from bam_index.c - -#define BAM_MIN_CHUNK_GAP 32768 -// 1<<14 is the size of minimum bin. -#define BAM_LIDX_SHIFT 14 -// =(8^6-1)/7+1 -#define MAX_BIN 37450 - -typedef struct -{ - uint64_t u, v; -} pair64_t; - -struct bam_fetch_iterator_t { - bam1_t * b; - pair64_t * off; - int n_off; - uint64_t curr_off; - int curr_chunk; - bamFile fp; - int tid; - int beg; - int end; - int n_seeks; -}; - -int is_overlap(uint32_t beg, uint32_t end, const bam1_t *b); +////////////////////////////////////////////////////////////////// +////////////////////////////////////////////////////////////////// +////////////////////////////////////////////////////////////////// +// code for iterator + +/*! @typedef + @Structure for holding current state (current alignment etc.) for iterating through + alignments overlapping a specified region. + @field b pointer to the current alignment + @field off pointer to an array of chunk loci (each with beg/end positions) + @field n_off The number of chunks + @field curr_off The current file positon + @field curr_chunk The item in a list of chunk + @discussion See also bam_fetch_iterate +*/ +struct __bam_fetch_iterator_t; +typedef struct __bam_fetch_iterator_t bam_fetch_iterator_t; + +/*! + @abstract Retrieve the alignments that are overlapped with the + specified region. + + @discussion Returns iterator object to retrieve successive alignments ordered by + start position. + @param fp BAM file handler + @param idx pointer to the alignment index + @param tid chromosome ID as is defined in the header + @param beg start coordinate, 0-based + @param end end coordinate, 0-based +*/ +bam_fetch_iterator_t * bam_init_fetch_iterator(bamFile fp, const bam_index_t *idx, int tid, int beg, int end); + + +/*! + @abstract Iterates through alignments overlapped the specified region. + @discussion Returns pointer to successive alignments ordered by start position. + Returns null pointer to signal the end of the iteration. + The alignment data is nested within the iterator to avoid unnecessary allocations. +*/ +bam1_t * bam_fetch_iterate(bam_fetch_iterator_t *iter); + +bam_fetch_iterator_t* bam_init_fetchall_iterator(bamFile fp, const bam_index_t *idx); +bam1_t * bam_fetchall_iterate(bam_fetch_iterator_t *iter); + +////////////////////////////////////////////////////////////////// +////////////////////////////////////////////////////////////////// +////////////////////////////////////////////////////////////////// +// various helper functions int pysam_bam_plbuf_push(const bam1_t *b, bam_plbuf_t *buf, int cont); @@ -42,4 +63,33 @@ int pysam_dispatch(int argc, char *argv[] ); // stand-in for macro - not wrappable in pyrex void pysam_bam_destroy1( bam1_t * b ); +// stand-in for other samtools macros +uint32_t * pysam_bam1_cigar( const bam1_t * b); +char * pysam_bam1_qname( const bam1_t * b); +uint8_t * pysam_bam1_seq( const bam1_t * b); +uint8_t * pysam_bam1_qual( const bam1_t * b); +uint8_t * pysam_bam1_aux( const bam1_t * b); + +/*! + @abstract Update the variable length data within a bam1_t entry + + Old data is deleted and the data within b are re-arranged to + make place for new data. + + @discussion Returns b + + @param b bam1_t data + @param nbytes_old size of old data + @param nbytes_new size of new data + @param pos position of data +*/ +bam1_t * pysam_bam_update( bam1_t * b, + const size_t nbytes_old, + const size_t nbytes_new, + uint8_t * pos ); + +// translate a nucleotide character to binary code +unsigned char pysam_translate_sequence( const unsigned char s ); + + #endif diff --git a/samtools/bam.c b/samtools/bam.c index 619b46a..ee7642b 100644 --- a/samtools/bam.c +++ b/samtools/bam.c @@ -1,9 +1,11 @@ #include #include +#include #include #include "bam.h" #include "bam_endian.h" #include "kstring.h" +#include "sam_header.h" int bam_is_be = 0; char *bam_flag2char_table = "pPuUrR12sfd\0\0\0\0\0"; @@ -12,28 +14,6 @@ char *bam_flag2char_table = "pPuUrR12sfd\0\0\0\0\0"; * CIGAR related routines * **************************/ -int bam_segreg(int32_t pos, const bam1_core_t *c, const uint32_t *cigar, bam_segreg_t *reg) -{ - unsigned k; - int32_t x = c->pos, y = 0; - int state = 0; - for (k = 0; k < c->n_cigar; ++k) { - int op = cigar[k] & BAM_CIGAR_MASK; // operation - int l = cigar[k] >> BAM_CIGAR_SHIFT; // length - if (state == 0 && (op == BAM_CMATCH || op == BAM_CDEL || op == BAM_CINS) && x + l > pos) { - reg->tbeg = x; reg->qbeg = y; reg->cbeg = k; - state = 1; - } - if (op == BAM_CMATCH) { x += l; y += l; } - else if (op == BAM_CDEL || op == BAM_CREF_SKIP) x += l; - else if (op == BAM_CINS || op == BAM_CSOFT_CLIP) y += l; - if (state == 1 && (op == BAM_CSOFT_CLIP || op == BAM_CHARD_CLIP || op == BAM_CREF_SKIP || k == c->n_cigar - 1)) { - reg->tend = x; reg->qend = y; reg->cend = k; - } - } - return state? 0 : -1; -} - uint32_t bam_calend(const bam1_core_t *c, const uint32_t *cigar) { uint32_t k, end; @@ -80,10 +60,9 @@ void bam_header_destroy(bam_header_t *header) free(header->target_len); } free(header->text); -#ifndef BAM_NO_HASH - if (header->rg2lib) bam_strmap_destroy(header->rg2lib); + if (header->dict) sam_header_free(header->dict); + if (header->rg2lib) sam_tbl_destroy(header->rg2lib); bam_destroy_header_hash(header); -#endif free(header); } @@ -94,7 +73,11 @@ bam_header_t *bam_header_read(bamFile fp) int32_t i = 1, name_len; // check EOF i = bgzf_check_EOF(fp); - if (i < 0) fprintf(stderr, "[bam_header_read] read from pipe; skip EOF checking.\n"); + if (i < 0) { + // If the file is a pipe, checking the EOF marker will *always* fail + // with ESPIPE. Suppress the error message in this case. + if (errno != ESPIPE) perror("[bam_header_read] bgzf_check_EOF"); + } else if (i == 0) fprintf(stderr, "[bam_header_read] EOF marker is absent.\n"); // read "BAM1" if (bam_read(fp, buf, 4) != 4) return 0; @@ -308,3 +291,13 @@ void bam_view1(const bam_header_t *header, const bam1_t *b) printf("%s\n", s); free(s); } + +// FIXME: we should also check the LB tag associated with each alignment +const char *bam_get_library(bam_header_t *h, const bam1_t *b) +{ + const uint8_t *rg; + if (h->dict == 0) h->dict = sam_header_parse2(h->text); + if (h->rg2lib == 0) h->rg2lib = sam_header2tbl(h->dict, "RG", "ID", "LB"); + rg = bam_aux_get(b, "RG"); + return (rg == 0)? 0 : sam_tbl_get(h->rg2lib, (const char*)(rg + 1)); +} diff --git a/samtools/bam.h b/samtools/bam.h index c739b62..291b303 100644 --- a/samtools/bam.h +++ b/samtools/bam.h @@ -73,6 +73,7 @@ typedef gzFile bamFile; @field n_targets number of reference sequences @field target_name names of the reference sequences @field target_len lengths of the referene sequences + @field dict header dictionary @field hash hash table for fast name lookup @field rg2lib hash table for @RG-ID -> LB lookup @field l_text length of the plain text in the header @@ -85,7 +86,7 @@ typedef struct { int32_t n_targets; char **target_name; uint32_t *target_len; - void *hash, *rg2lib; + void *dict, *hash, *rg2lib; int l_text; char *text; } bam_header_t; @@ -423,8 +424,8 @@ extern "C" { @abstract Free the memory allocated for an alignment. @param b pointer to an alignment */ -#define bam_destroy1(b) do { \ - free((b)->data); free(b); \ +#define bam_destroy1(b) do { \ + if (b) { free((b)->data); free(b); } \ } while (0) /*! @@ -437,6 +438,8 @@ extern "C" { char *bam_format1_core(const bam_header_t *header, const bam1_t *b, int of); + const char *bam_get_library(bam_header_t *header, const bam1_t *b); + /*! @typedef @abstract Structure for one alignment covering the pileup position. @field b pointer to the alignment @@ -579,51 +582,6 @@ extern "C" { */ int bam_fetch(bamFile fp, const bam_index_t *idx, int tid, int beg, int end, void *data, bam_fetch_f func); - /*! @typedef - @Structure for holding current state (current alignment etc.) for iterating through - alignments overlapping a specified region. - @field b pointer to the current alignment - @field off pointer to an array of chunk loci (each with beg/end positions) - @field n_off The number of chunks - @field curr_off The current file positon - @field curr_chunk The item in a list of chunk - - @discussion See also bam_fetch_iterate - */ - struct __bam_fetch_iterator_t; - typedef struct __bam_fetch_iterator_t bam_fetch_iterator_t; - - - - /*! - @abstract Retrieve the alignments that are overlapped with the - specified region. - - @discussion Returns iterator object to retrieve successive alignments ordered by - start position. - @param fp BAM file handler - @param idx pointer to the alignment index - @param tid chromosome ID as is defined in the header - @param beg start coordinate, 0-based - @param end end coordinate, 0-based - */ - bam_fetch_iterator_t * bam_init_fetch_iterator(bamFile fp, const bam_index_t *idx, int tid, int beg, int end); - - - /*! - @abstract Iterates through alignments overlapped the specified region. - @discussion Returns pointer to successive alignments ordered by start position. - Returns null pointer to signal the end of the iteration. - The alignment data is nested within the iterator to avoid unnecessary allocations. - */ - bam1_t * bam_fetch_iterate(bam_fetch_iterator_t *iter); - - /*! - @abstract Cleans up the iterator data - @discussion Destroys data allocated on the heap - */ - void bam_cleanup_fetch_iterator(bam_fetch_iterator_t *iter); - /*! @abstract Parse a region in the format: "chr2:100,000-200,000". @discussion bam_header_t::hash will be initialized if empty. @@ -677,14 +635,6 @@ extern "C" { */ int32_t bam_cigar2qlen(const bam1_core_t *c, const uint32_t *cigar); - typedef struct { - int32_t qbeg, qend; - int32_t tbeg, tend; - int32_t cbeg, cend; - } bam_segreg_t; - - int bam_segreg(int32_t pos, const bam1_core_t *c, const uint32_t *cigar, bam_segreg_t *reg); - #ifdef __cplusplus } #endif @@ -744,7 +694,4 @@ static inline bam1_t *bam_dup1(const bam1_t *src) return b; } -bam_fetch_iterator_t* bam_init_fetchall_iterator(bamFile fp, const bam_index_t *idx); -bam1_t * bam_fetchall_iterate(bam_fetch_iterator_t *iter); - #endif diff --git a/samtools/bam_aux.c b/samtools/bam_aux.c index d0d733f..89e99f2 100644 --- a/samtools/bam_aux.c +++ b/samtools/bam_aux.c @@ -180,70 +180,3 @@ char *bam_aux2Z(const uint8_t *s) if (type == 'Z' || type == 'H') return (char*)s; else return 0; } - -/****************** - * rg2lib related * - ******************/ - -int bam_strmap_put(void *rg2lib, const char *rg, const char *lib) -{ - int ret; - khint_t k; - khash_t(r2l) *h = (khash_t(r2l)*)rg2lib; - char *key; - if (h == 0) return 1; - key = strdup(rg); - k = kh_put(r2l, h, key, &ret); - if (ret) kh_val(h, k) = strdup(lib); - else { - fprintf(stderr, "[bam_rg2lib_put] duplicated @RG ID: %s\n", rg); - free(key); - } - return 0; -} - -const char *bam_strmap_get(const void *rg2lib, const char *rg) -{ - const khash_t(r2l) *h = (const khash_t(r2l)*)rg2lib; - khint_t k; - if (h == 0) return 0; - k = kh_get(r2l, h, rg); - if (k != kh_end(h)) return (const char*)kh_val(h, k); - else return 0; -} - -void *bam_strmap_dup(const void *rg2lib) -{ - const khash_t(r2l) *h = (const khash_t(r2l)*)rg2lib; - khash_t(r2l) *g; - khint_t k, l; - int ret; - if (h == 0) return 0; - g = kh_init(r2l); - for (k = kh_begin(h); k < kh_end(h); ++k) { - if (kh_exist(h, k)) { - char *key = strdup(kh_key(h, k)); - l = kh_put(r2l, g, key, &ret); - kh_val(g, l) = strdup(kh_val(h, k)); - } - } - return g; -} - -void *bam_strmap_init() -{ - return (void*)kh_init(r2l); -} - -void bam_strmap_destroy(void *rg2lib) -{ - khash_t(r2l) *h = (khash_t(r2l)*)rg2lib; - khint_t k; - if (h == 0) return; - for (k = kh_begin(h); k < kh_end(h); ++k) { - if (kh_exist(h, k)) { - free((char*)kh_key(h, k)); free(kh_val(h, k)); - } - } - kh_destroy(r2l, h); -} diff --git a/samtools/bam_import.c b/samtools/bam_import.c index 1dc906e..9d463d1 100644 --- a/samtools/bam_import.c +++ b/samtools/bam_import.c @@ -10,6 +10,7 @@ #endif #include "kstring.h" #include "bam.h" +#include "sam_header.h" #include "kseq.h" #include "khash.h" @@ -170,107 +171,26 @@ static inline void append_text(bam_header_t *header, kstring_t *str) header->text[header->l_text] = 0; } -int sam_header_parse_rg(bam_header_t *h) -{ - kstring_t *rgid, *rglib; - char *p, *q, *s, *r; - int n = 0; - - // free - if (h == 0) return 0; - bam_strmap_destroy(h->rg2lib); h->rg2lib = 0; - if (h->l_text < 3) return 0; - // parse @RG lines - h->rg2lib = bam_strmap_init(); - rgid = calloc(1, sizeof(kstring_t)); - rglib = calloc(1, sizeof(kstring_t)); - s = h->text; - while ((s = strstr(s, "@RG")) != 0) { - if (rgid->l && rglib->l) { - bam_strmap_put(h->rg2lib, rgid->s, rglib->s); - ++n; - } - rgid->l = rglib->l = 0; - s += 3; - r = s; - if ((p = strstr(s, "ID:")) != 0) { - q = p + 3; - for (p = q; *p && *p != '\t' && *p != '\r' && *p != '\n'; ++p); - kputsn(q, p - q, rgid); - } else { - fprintf(stderr, "[bam_header_parse] missing ID tag in @RG lines.\n"); - break; - } - if (r < p) r = p; - if ((p = strstr(s, "LB:")) != 0) { - q = p + 3; - for (p = q; *p && *p != '\t' && *p != '\r' && *p != '\n'; ++p); - kputsn(q, p - q, rglib); - } else { - fprintf(stderr, "[bam_header_parse] missing LB tag in @RG lines.\n"); - break; - } - if (r < p) r = p; - s = r + 3; - } - if (rgid->l && rglib->l) { - bam_strmap_put(h->rg2lib, rgid->s, rglib->s); - ++n; - } - free(rgid->s); free(rgid); - free(rglib->s); free(rglib); - if (n == 0) { - bam_strmap_destroy(h->rg2lib); - h->rg2lib = 0; - } - return n; -} - int sam_header_parse(bam_header_t *h) { + char **tmp; int i; - char *s, *p, *q, *r; - - // free free(h->target_len); free(h->target_name); h->n_targets = 0; h->target_len = 0; h->target_name = 0; if (h->l_text < 3) return 0; - // count number of @SQ - s = h->text; - while ((s = strstr(s, "@SQ")) != 0) { - ++h->n_targets; - s += 3; - } + if (h->dict == 0) h->dict = sam_header_parse2(h->text); + tmp = sam_header2list(h->dict, "SQ", "SN", &h->n_targets); if (h->n_targets == 0) return 0; - h->target_len = (uint32_t*)calloc(h->n_targets, 4); - h->target_name = (char**)calloc(h->n_targets, sizeof(void*)); - // parse @SQ lines - i = 0; - s = h->text; - while ((s = strstr(s, "@SQ")) != 0) { - s += 3; - r = s; - if ((p = strstr(s, "SN:")) != 0) { - q = p + 3; - for (p = q; *p && *p != '\t' && *p != '\r' && *p != '\n'; ++p); - h->target_name[i] = (char*)calloc(p - q + 1, 1); - strncpy(h->target_name[i], q, p - q); - } else goto header_err_ret; - if (r < p) r = p; - if ((p = strstr(s, "LN:")) != 0) h->target_len[i] = strtol(p + 3, 0, 10); - else goto header_err_ret; - if (r < p) r = p; - s = r + 3; - ++i; - } - sam_header_parse_rg(h); + h->target_name = calloc(h->n_targets, sizeof(void*)); + for (i = 0; i < h->n_targets; ++i) + h->target_name[i] = strdup(tmp[i]); + free(tmp); + tmp = sam_header2list(h->dict, "SQ", "LN", &h->n_targets); + h->target_len = calloc(h->n_targets, 4); + for (i = 0; i < h->n_targets; ++i) + h->target_len[i] = atoi(tmp[i]); + free(tmp); return h->n_targets; - -header_err_ret: - fprintf(stderr, "[bam_header_parse] missing SN or LN tag in @SQ lines.\n"); - free(h->target_len); free(h->target_name); - h->n_targets = 0; h->target_len = 0; h->target_name = 0; - return 0; } bam_header_t *sam_header_read(tamFile fp) diff --git a/samtools/bam_index.c b/samtools/bam_index.c index 26ea53c..a627884 100644 --- a/samtools/bam_index.c +++ b/samtools/bam_index.c @@ -538,125 +538,37 @@ pair64_t * get_chunk_coordinates(const bam_index_t *idx, int tid, int beg, int e return off; } - -//int bam_fetch(bamFile fp, const bam_index_t *idx, int tid, int beg, int end, void *data, bam_fetch_f func) -//{ -// int n_off; -// pair64_t *off = get_chunk_coordinates(idx, tid, beg, end, &n_off); -// { -// // retrive alignments -// uint64_t curr_off; -// int i, ret, n_seeks; -// n_seeks = 0; i = -1; curr_off = 0; -// bam1_t *b = (bam1_t*)calloc(1, sizeof(bam1_t)); -// for (;;) { -// if (curr_off == 0 || curr_off >= off[i].v) { // then jump to the next chunk -// if (i == n_off - 1) break; // no more chunks -// if (i >= 0) assert(curr_off == off[i].v); // otherwise bug -// if (i < 0 || off[i].v != off[i+1].u) { // not adjacent chunks; then seek -// bam_seek(fp, off[i+1].u, SEEK_SET); -// curr_off = bam_tell(fp); -// ++n_seeks; -// } -// ++i; -// } -// if ((ret = bam_read1(fp, b)) > 0) { -// curr_off = bam_tell(fp); -// if (b->core.tid != tid || b->core.pos >= end) break; // no need to proceed -// else if (is_overlap(beg, end, b)) func(b, data); -// } else break; // end of file -// } -//// fprintf(stderr, "[bam_fetch] # seek calls: %d\n", n_seeks); -// bam_destroy1(b); -// } -// free(off); -// return 0; -//} - -struct __bam_fetch_iterator_t{ - bam1_t * b; - pair64_t * off; - int n_off; - uint64_t curr_off; - int curr_chunk; - bamFile fp; - int tid; - int beg; - int end; - int n_seeks; -}; - - - -bam_fetch_iterator_t* bam_init_fetch_iterator(bamFile fp, const bam_index_t *idx, int tid, int beg, int end) -{ - // iterator contains current alignment position - // and will contain actual alignment during iterations - bam_fetch_iterator_t* iter = (bam_fetch_iterator_t*)calloc(1, sizeof(bam_fetch_iterator_t)); - iter->b = (bam1_t*)calloc(1, sizeof(bam1_t)); - - // list of chunks containing our alignments - iter->off = get_chunk_coordinates(idx, tid, beg, end, &iter->n_off); - - // initialise other state variables in iterator - iter->fp = fp; - iter->curr_chunk = -1; - iter->curr_off = 0; - iter->n_seeks = 0; - iter->tid = tid; - iter->beg = beg; - iter->end = end; - return iter; -} - -bam1_t * bam_fetch_iterate(bam_fetch_iterator_t *iter) +int bam_fetch(bamFile fp, const bam_index_t *idx, int tid, int beg, int end, void *data, bam_fetch_f func) { - if (!iter->off) { - return 0; - } - - int ret; - // iterate through all alignments in chunks - for (;;) { - if (iter->curr_off == 0 || iter->curr_off >= iter->off[iter->curr_chunk].v) { // then jump to the next chunk - if (iter->curr_chunk == iter->n_off - 1) break; // no more chunks - if (iter->curr_chunk >= 0) assert(iter->curr_off == iter->off[iter->curr_chunk].v); // otherwise bug - if (iter->curr_chunk < 0 || iter->off[iter->curr_chunk].v != iter->off[iter->curr_chunk+1].u) { // not adjacent chunks; then seek - bam_seek(iter->fp, iter->off[iter->curr_chunk+1].u, SEEK_SET); - iter->curr_off = bam_tell(iter->fp); - ++iter->n_seeks; + int n_off; + pair64_t *off = get_chunk_coordinates(idx, tid, beg, end, &n_off); + if (off == 0) return 0; + { + // retrive alignments + uint64_t curr_off; + int i, ret, n_seeks; + n_seeks = 0; i = -1; curr_off = 0; + bam1_t *b = (bam1_t*)calloc(1, sizeof(bam1_t)); + for (;;) { + if (curr_off == 0 || curr_off >= off[i].v) { // then jump to the next chunk + if (i == n_off - 1) break; // no more chunks + if (i >= 0) assert(curr_off == off[i].v); // otherwise bug + if (i < 0 || off[i].v != off[i+1].u) { // not adjacent chunks; then seek + bam_seek(fp, off[i+1].u, SEEK_SET); + curr_off = bam_tell(fp); + ++n_seeks; + } + ++i; } - ++iter->curr_chunk; + if ((ret = bam_read1(fp, b)) > 0) { + curr_off = bam_tell(fp); + if (b->core.tid != tid || b->core.pos >= end) break; // no need to proceed + else if (is_overlap(beg, end, b)) func(b, data); + } else break; // end of file } - if ((ret = bam_read1(iter->fp, iter->b)) > 0) { - iter->curr_off = bam_tell(iter->fp); - if (iter->b->core.tid != iter->tid || iter->b->core.pos >= iter->end) break; // no need to proceed - else if (is_overlap(iter->beg, iter->end, iter->b)) - // - //func(iter->b, data); - // - return iter->b; - } else - return 0; // end of file +// fprintf(stderr, "[bam_fetch] # seek calls: %d\n", n_seeks); + bam_destroy1(b); } + free(off); return 0; } - -void bam_cleanup_fetch_iterator(bam_fetch_iterator_t *iter) -{ -// fprintf(stderr, "[bam_fetch] # seek calls: %d\n", iter->n_seeks); - bam_destroy1(iter->b); - free(iter->off); -} - - -int bam_fetch(bamFile fp, const bam_index_t *idx, int tid, int beg, int end, void *data, bam_fetch_f func) -{ - bam_fetch_iterator_t* iter = bam_init_fetch_iterator(fp, idx, tid, beg, end); - bam1_t * b; - while (b = bam_fetch_iterate(iter)) - func(b, data); - bam_cleanup_fetch_iterator(iter); - return 0; -} - diff --git a/samtools/bam_maqcns.c b/samtools/bam_maqcns.c index f36b0ee..71c2185 100644 --- a/samtools/bam_maqcns.c +++ b/samtools/bam_maqcns.c @@ -1,10 +1,13 @@ #include +#include #include "bam.h" #include "bam_maqcns.h" #include "ksort.h" +#include "kaln.h" KSORT_INIT_GENERIC(uint32_t) -#define MAX_WINDOW 33 +#define INDEL_WINDOW_SIZE 50 +#define INDEL_EXT_DEP 0.9 typedef struct __bmc_aux_t { int max; @@ -22,14 +25,13 @@ char bam_nt16_nt4_table[] = { 4, 0, 1, 4, 2, 4, 4, 4, 3, 4, 4, 4, 4, 4, 4, 4 }; /* P() = \theta \sum_{i=1}^{N-1} 1/i P(D|) = \sum_{k=1}^{N-1} p_k 1/2 [(k/N)^n_2(1-k/N)^n_1 + (k/N)^n1(1-k/N)^n_2] - p_k = i/k / \sum_{i=1}^{N-1} 1/i + p_k = 1/k / \sum_{i=1}^{N-1} 1/i */ static void cal_het(bam_maqcns_t *aa) { int k, n1, n2; double sum_harmo; // harmonic sum double poly_rate; - double p1 = 0.0, p3 = 0.0; // just for testing free(aa->lhet); aa->lhet = (double*)calloc(256 * 256, sizeof(double)); @@ -39,7 +41,7 @@ static void cal_het(bam_maqcns_t *aa) for (n1 = 0; n1 < 256; ++n1) { for (n2 = 0; n2 < 256; ++n2) { long double sum = 0.0; - double lC = lgamma(n1+n2+1) - lgamma(n1+1) - lgamma(n2+1); // \binom{n1+n2}{n1} + double lC = aa->is_soap? 0 : lgamma(n1+n2+1) - lgamma(n1+1) - lgamma(n2+1); // \binom{n1+n2}{n1} for (k = 1; k <= aa->n_hap - 1; ++k) { double pk = 1.0 / k / sum_harmo; double log1 = log((double)k/aa->n_hap); @@ -47,8 +49,6 @@ static void cal_het(bam_maqcns_t *aa) sum += pk * 0.5 * (expl(log1*n2) * expl(log2*n1) + expl(log1*n1) * expl(log2*n2)); } aa->lhet[n1<<8|n2] = lC + logl(sum); - if (n1 == 17 && n2 == 3) p3 = lC + logl(expl(logl(0.5) * 20)); - if (n1 == 19 && n2 == 1) p1 = lC + logl(expl(logl(0.5) * 20)); } } poly_rate = aa->het_rate * sum_harmo; @@ -62,16 +62,18 @@ static void cal_coef(bam_maqcns_t *aa) long double sum_a[257], b[256], q_c[256], tmp[256], fk2[256]; double *lC; - lC = (double*)calloc(256 * 256, sizeof(double)); // aa->lhet will be allocated and initialized free(aa->fk); free(aa->coef); + aa->coef = 0; aa->fk = (double*)calloc(256, sizeof(double)); - aa->coef = (double*)calloc(256*256*64, sizeof(double)); aa->fk[0] = fk2[0] = 1.0; for (n = 1; n != 256; ++n) { aa->fk[n] = pow(aa->theta, n) * (1.0 - aa->eta) + aa->eta; fk2[n] = aa->fk[n>>1]; // this is an approximation, assuming reads equally likely come from both strands } + if (aa->is_soap) return; + aa->coef = (double*)calloc(256*256*64, sizeof(double)); + lC = (double*)calloc(256 * 256, sizeof(double)); for (n = 1; n != 256; ++n) for (k = 1; k <= n; ++k) lC[n<<8|k] = lgamma(n+1) - lgamma(k+1) - lgamma(n-k+1); @@ -170,7 +172,7 @@ glf1_t *bam_maqcns_glfgen(int _n, const bam_pileup1_t *pl, uint8_t ref_base, bam if (w[k] < 0xff) ++w[k]; ++b->c[k&3]; } - tmp = (int)(info&0x7f) < bm->cap_mapQ? (int)(info&0x7f) : bm->cap_mapQ; + tmp = (int)(info&0xff) < bm->cap_mapQ? (int)(info&0xff) : bm->cap_mapQ; rms += tmp * tmp; } b->rms_mapQ = (uint8_t)(sqrt((double)rms / n) + .499); @@ -180,56 +182,73 @@ glf1_t *bam_maqcns_glfgen(int _n, const bam_pileup1_t *pl, uint8_t ref_base, bam for (j = 0; j != 4; ++j) b->c[j] = (int)(254.0 * b->c[j] / c + 0.5); for (j = c = 0; j != 4; ++j) c += b->c[j]; } - // generate likelihood - for (j = 0; j != 4; ++j) { - // homozygous - float tmp1, tmp3; - int tmp2, bar_e; - for (k = 0, tmp1 = tmp3 = 0.0, tmp2 = 0; k != 4; ++k) { - if (j == k) continue; - tmp1 += b->esum[k]; tmp2 += b->c[k]; tmp3 += b->fsum[k]; - } - if (tmp2) { - bar_e = (int)(tmp1 / tmp3 + 0.5); - if (bar_e < 4) bar_e = 4; // should not happen - if (bar_e > 63) bar_e = 63; - p[j<<2|j] = tmp1 + bm->coef[bar_e<<16|c<<8|tmp2]; - } else p[j<<2|j] = 0.0; // all the bases are j - // heterozygous - for (k = j + 1; k < 4; ++k) { - for (i = 0, tmp2 = 0, tmp1 = tmp3 = 0.0; i != 4; ++i) { - if (i == j || i == k) continue; - tmp1 += b->esum[i]; tmp2 += b->c[i]; tmp3 += b->fsum[i]; + if (!bm->is_soap) { + // generate likelihood + for (j = 0; j != 4; ++j) { + // homozygous + float tmp1, tmp3; + int tmp2, bar_e; + for (k = 0, tmp1 = tmp3 = 0.0, tmp2 = 0; k != 4; ++k) { + if (j == k) continue; + tmp1 += b->esum[k]; tmp2 += b->c[k]; tmp3 += b->fsum[k]; } if (tmp2) { bar_e = (int)(tmp1 / tmp3 + 0.5); - if (bar_e < 4) bar_e = 4; + if (bar_e < 4) bar_e = 4; // should not happen if (bar_e > 63) bar_e = 63; - p[j<<2|k] = p[k<<2|j] = -4.343 * bm->lhet[b->c[j]<<8|b->c[k]] + tmp1 + bm->coef[bar_e<<16|c<<8|tmp2]; - } else p[j<<2|k] = p[k<<2|j] = -4.343 * bm->lhet[b->c[j]<<8|b->c[k]]; // all the bases are either j or k + p[j<<2|j] = tmp1 + bm->coef[bar_e<<16|c<<8|tmp2]; + } else p[j<<2|j] = 0.0; // all the bases are j + // heterozygous + for (k = j + 1; k < 4; ++k) { + for (i = 0, tmp2 = 0, tmp1 = tmp3 = 0.0; i != 4; ++i) { + if (i == j || i == k) continue; + tmp1 += b->esum[i]; tmp2 += b->c[i]; tmp3 += b->fsum[i]; + } + if (tmp2) { + bar_e = (int)(tmp1 / tmp3 + 0.5); + if (bar_e < 4) bar_e = 4; + if (bar_e > 63) bar_e = 63; + p[j<<2|k] = p[k<<2|j] = -4.343 * bm->lhet[b->c[j]<<8|b->c[k]] + tmp1 + bm->coef[bar_e<<16|c<<8|tmp2]; + } else p[j<<2|k] = p[k<<2|j] = -4.343 * bm->lhet[b->c[j]<<8|b->c[k]]; // all the bases are either j or k + } + // + for (k = 0; k != 4; ++k) + if (p[j<<2|k] < 0.0) p[j<<2|k] = 0.0; } - // - for (k = 0; k != 4; ++k) - if (p[j<<2|k] < 0.0) p[j<<2|k] = 0.0; - } - { // fix p[k<<2|k] - float max1, max2, min1, min2; - int max_k, min_k; - max_k = min_k = -1; - max1 = max2 = -1.0; min1 = min2 = 1e30; - for (k = 0; k < 4; ++k) { - if (b->esum[k] > max1) { - max2 = max1; max1 = b->esum[k]; max_k = k; - } else if (b->esum[k] > max2) max2 = b->esum[k]; + { // fix p[k<<2|k] + float max1, max2, min1, min2; + int max_k, min_k; + max_k = min_k = -1; + max1 = max2 = -1.0; min1 = min2 = 1e30; + for (k = 0; k < 4; ++k) { + if (b->esum[k] > max1) { + max2 = max1; max1 = b->esum[k]; max_k = k; + } else if (b->esum[k] > max2) max2 = b->esum[k]; + } + for (k = 0; k < 4; ++k) { + if (p[k<<2|k] < min1) { + min2 = min1; min1 = p[k<<2|k]; min_k = k; + } else if (p[k<<2|k] < min2) min2 = p[k<<2|k]; + } + if (max1 > max2 && (min_k != max_k || min1 + 1.0 > min2)) + p[max_k<<2|max_k] = min1 > 1.0? min1 - 1.0 : 0.0; } - for (k = 0; k < 4; ++k) { - if (p[k<<2|k] < min1) { - min2 = min1; min1 = p[k<<2|k]; min_k = k; - } else if (p[k<<2|k] < min2) min2 = p[k<<2|k]; + } else { // apply the SOAP model + // generate likelihood + for (j = 0; j != 4; ++j) { + float tmp; + // homozygous + for (k = 0, tmp = 0.0; k != 4; ++k) + if (j != k) tmp += b->esum[k]; + p[j<<2|j] = tmp; + // heterozygous + for (k = j + 1; k < 4; ++k) { + for (i = 0, tmp = 0.0; i != 4; ++i) + if (i != j && i != k) tmp += b->esum[i]; + p[j<<2|k] = p[k<<2|j] = -4.343 * bm->lhet[b->c[j]<<8|b->c[k]] + tmp; + } } - if (max1 > max2 && (min_k != max_k || min1 + 1.0 > min2)) - p[max_k<<2|max_k] = min1 > 1.0? min1 - 1.0 : 0.0; } // convert necessary information to glf1_t @@ -304,12 +323,40 @@ void bam_maqindel_ret_destroy(bam_maqindel_ret_t *mir) free(mir->s[0]); free(mir->s[1]); free(mir); } +int bam_tpos2qpos(const bam1_core_t *c, const uint32_t *cigar, int32_t tpos, int is_left, int32_t *_tpos) +{ + int k, x = c->pos, y = 0, last_y = 0; + *_tpos = c->pos; + for (k = 0; k < c->n_cigar; ++k) { + int op = cigar[k] & BAM_CIGAR_MASK; + int l = cigar[k] >> BAM_CIGAR_SHIFT; + if (op == BAM_CMATCH) { + if (c->pos > tpos) return y; + if (x + l > tpos) { + *_tpos = tpos; + return y + (tpos - x); + } + x += l; y += l; + last_y = y; + } else if (op == BAM_CINS || op == BAM_CSOFT_CLIP) y += l; + else if (op == BAM_CDEL || op == BAM_CREF_SKIP) { + if (x + l > tpos) { + *_tpos = is_left? x : x + l; + return y; + } + x += l; + } + } + *_tpos = x; + return last_y; +} + #define MINUS_CONST 0x10000000 bam_maqindel_ret_t *bam_maqindel(int n, int pos, const bam_maqindel_opt_t *mi, const bam_pileup1_t *pl, const char *ref, int _n_types, int *_types) { - int i, j, n_types, *types, left, right; + int i, j, n_types, *types, left, right, max_rd_len = 0; bam_maqindel_ret_t *ret = 0; // if there is no proposed indel, check if there is an indel from the alignment if (_n_types == 0) { @@ -329,6 +376,8 @@ bam_maqindel_ret_t *bam_maqindel(int n, int pos, const bam_maqindel_opt_t *mi, c const bam_pileup1_t *p = pl + i; if (!(p->b->core.flag&BAM_FUNMAP) && p->indel != 0) aux[m++] = MINUS_CONST + p->indel; + j = bam_cigar2qlen(&p->b->core, bam1_cigar(p->b)); + if (j > max_rd_len) max_rd_len = j; } if (_n_types) // then also add this to aux[] for (i = 0; i < _n_types; ++i) @@ -347,23 +396,17 @@ bam_maqindel_ret_t *bam_maqindel(int n, int pos, const bam_maqindel_opt_t *mi, c free(aux); } { // calculate left and right boundary - bam_segreg_t seg; - left = 0x7fffffff; right = 0; - for (i = 0; i < n; ++i) { - const bam_pileup1_t *p = pl + i; - if (!(p->b->core.flag&BAM_FUNMAP)) { - bam_segreg(pos, &p->b->core, bam1_cigar(p->b), &seg); - if (seg.tbeg < left) left = seg.tbeg; - if (seg.tend > right) right = seg.tend; - } - } - if (pos - left > MAX_WINDOW) left = pos - MAX_WINDOW; - if (right - pos> MAX_WINDOW) right = pos + MAX_WINDOW; + left = pos > INDEL_WINDOW_SIZE? pos - INDEL_WINDOW_SIZE : 0; + right = pos + INDEL_WINDOW_SIZE; + if (types[0] < 0) right -= types[0]; + // in case the alignments stand out the reference + for (i = pos; i < right; ++i) + if (ref[i] == 0) break; + right = i; } { // the core part - char *ref2, *inscns = 0; + char *ref2, *rs, *inscns = 0; int k, l, *score, *pscore, max_ins = types[n_types-1]; - ref2 = (char*)calloc(right - left + types[n_types-1] + 2, 1); if (max_ins > 0) { // get the consensus of inserted sequences int *inscns_aux = (int*)calloc(4 * n_types * max_ins, sizeof(int)); // count occurrences @@ -396,52 +439,68 @@ bam_maqindel_ret_t *bam_maqindel(int n, int pos, const bam_maqindel_opt_t *mi, c free(inscns_aux); } // calculate score + ref2 = (char*)calloc(right - left + types[n_types-1] + 2, 1); + rs = (char*)calloc(right - left + max_rd_len + types[n_types-1] + 2, 1); score = (int*)calloc(n_types * n, sizeof(int)); pscore = (int*)calloc(n_types * n, sizeof(int)); for (i = 0; i < n_types; ++i) { + ka_param_t ap = ka_param_blast; + ap.band_width = 2 * types[n_types - 1] + 2; // write ref2 for (k = 0, j = left; j <= pos; ++j) - ref2[k++] = bam_nt16_table[(int)ref[j]]; + ref2[k++] = bam_nt16_nt4_table[bam_nt16_table[(int)ref[j]]]; if (types[i] <= 0) j += -types[i]; else for (l = 0; l < types[i]; ++l) - ref2[k++] = inscns[i*max_ins + l]; + ref2[k++] = bam_nt16_nt4_table[(int)inscns[i*max_ins + l]]; for (; j < right && ref[j]; ++j) - ref2[k++] = bam_nt16_table[(int)ref[j]]; + ref2[k++] = bam_nt16_nt4_table[bam_nt16_table[(int)ref[j]]]; + if (j < right) right = j; // calculate score for each read for (j = 0; j < n; ++j) { const bam_pileup1_t *p = pl + j; - uint32_t *cigar; - bam1_core_t *c = &p->b->core; - int s, ps; - bam_segreg_t seg; - if (c->flag&BAM_FUNMAP) continue; - cigar = bam1_cigar(p->b); - bam_segreg(pos, c, cigar, &seg); - for (ps = s = 0, l = seg.qbeg; c->pos + l < right && l < seg.qend; ++l) { - int cq = bam1_seqi(bam1_seq(p->b), l), ct; - // in the following line, "<" will happen if reads are too long - ct = c->pos + l - seg.qbeg >= left? ref2[c->pos + l - seg.qbeg - left] : 15; - if (cq < 15 && ct < 15) { - s += cq == ct? 1 : -mi->mm_penalty; - if (cq != ct) ps += bam1_qual(p->b)[l]; + int qbeg, qend, tbeg, tend; + if (p->b->core.flag & BAM_FUNMAP) continue; + qbeg = bam_tpos2qpos(&p->b->core, bam1_cigar(p->b), left, 0, &tbeg); + qend = bam_tpos2qpos(&p->b->core, bam1_cigar(p->b), right, 1, &tend); + assert(tbeg >= left); + for (l = qbeg; l < qend; ++l) + rs[l - qbeg] = bam_nt16_nt4_table[bam1_seqi(bam1_seq(p->b), l)]; + { + int x, y, n_acigar, ps; + uint32_t *acigar; + ps = 0; + if (tend - tbeg + types[i] <= 0) { + score[i*n+j] = -(1<<20); + pscore[i*n+j] = 1<<20; + continue; } - } - score[i*n + j] = s; pscore[i*n + j] = ps; - if (types[i] != 0) { // then try the other way to calculate the score - for (ps = s = 0, l = seg.qbeg; c->pos + l + types[i] < right && l < seg.qend; ++l) { - int cq = bam1_seqi(bam1_seq(p->b), l), ct; - ct = c->pos + l - seg.qbeg + types[i] >= left? ref2[c->pos + l - seg.qbeg + types[i] - left] : 15; - if (cq < 15 && ct < 15) { - s += cq == ct? 1 : -mi->mm_penalty; - if (cq != ct) ps += bam1_qual(p->b)[l]; + acigar = ka_global_core((uint8_t*)ref2 + tbeg - left, tend - tbeg + types[i], (uint8_t*)rs, qend - qbeg, &ap, &score[i*n+j], &n_acigar); + x = tbeg - left; y = 0; + for (l = 0; l < n_acigar; ++l) { + int op = acigar[l]&0xf; + int len = acigar[l]>>4; + if (op == BAM_CMATCH) { + int k; + for (k = 0; k < len; ++k) + if (ref2[x+k] != rs[y+k]) ps += bam1_qual(p->b)[y+k]; + x += len; y += len; + } else if (op == BAM_CINS || op == BAM_CSOFT_CLIP) { + if (op == BAM_CINS) ps += mi->q_indel * len; + y += len; + } else if (op == BAM_CDEL) { + ps += mi->q_indel * len; + x += len; } } + pscore[i*n+j] = ps; + /*if (pos == 2618517) { // for debugging only + fprintf(stderr, "pos=%d, type=%d, j=%d, score=%d, psore=%d, %d, %d, %d, %d, ", pos+1, types[i], j, score[i*n+j], pscore[i*n+j], tbeg, tend, qbeg, qend); + for (l = 0; l < n_acigar; ++l) fprintf(stderr, "%d%c", acigar[l]>>4, "MIDS"[acigar[l]&0xf]); fprintf(stderr, "\n"); + for (l = 0; l < tend - tbeg + types[i]; ++l) fputc("ACGTN"[ref2[l]], stderr); fputc('\n', stderr); + for (l = 0; l < qend - qbeg; ++l) fputc("ACGTN"[rs[l]], stderr); fputc('\n', stderr); + }*/ + free(acigar); } - if (score[i*n+j] < s) score[i*n+j] = s; // choose the higher of the two scores - if (pscore[i*n+j] > ps) pscore[i*n+j] = ps; - //if (types[i] != 0) score[i*n+j] -= mi->indel_err; - //printf("%d, %d, %d, %d, %d, %d, %d\n", p->b->core.pos + 1, seg.qbeg, i, types[i], j, - // score[i*n+j], pscore[i*n+j]); } } { // get final result @@ -491,13 +550,20 @@ bam_maqindel_ret_t *bam_maqindel(int n, int pos, const bam_maqindel_opt_t *mi, c else if (p->indel == ret->indel2) ++ret->cnt2; else ++ret->cnt_anti; } - // write gl[] - ret->gl[0] = ret->gl[1] = 0; - for (j = 0; j < n; ++j) { - int s1 = pscore[max1_i*n + j], s2 = pscore[max2_i*n + j]; - //printf("%d, %d, %d, %d, %d\n", pl[j].b->core.pos+1, max1_i, max2_i, s1, s2); - if (s1 > s2) ret->gl[0] += s1 - s2 < mi->q_indel? s1 - s2 : mi->q_indel; - else ret->gl[1] += s2 - s1 < mi->q_indel? s2 - s1 : mi->q_indel; + { // write gl[] + int tmp, seq_err = 0; + double x = 1.0; + tmp = max1_i - max2_i; + if (tmp < 0) tmp = -tmp; + for (j = 0; j < tmp + 1; ++j) x *= INDEL_EXT_DEP; + seq_err = mi->q_indel * (1.0 - x) / (1.0 - INDEL_EXT_DEP); + ret->gl[0] = ret->gl[1] = 0; + for (j = 0; j < n; ++j) { + int s1 = pscore[max1_i*n + j], s2 = pscore[max2_i*n + j]; + //printf("%d, %d, %d, %d, %d\n", pl[j].b->core.pos+1, max1_i, max2_i, s1, s2); + if (s1 > s2) ret->gl[0] += s1 - s2 < seq_err? s1 - s2 : seq_err; + else ret->gl[1] += s2 - s1 < seq_err? s2 - s1 : seq_err; + } } // write cnt_ref and cnt_ambi if (max1_i != 0 && max2_i != 0) { @@ -509,7 +575,7 @@ bam_maqindel_ret_t *bam_maqindel(int n, int pos, const bam_maqindel_opt_t *mi, c } } } - free(score); free(pscore); free(ref2); free(inscns); + free(score); free(pscore); free(ref2); free(rs); free(inscns); } { // call genotype int q[3], qr_indel = (int)(-4.343 * log(mi->r_indel) + 0.5); diff --git a/samtools/bam_maqcns.h b/samtools/bam_maqcns.h index 2a82aee..fa5489d 100644 --- a/samtools/bam_maqcns.h +++ b/samtools/bam_maqcns.h @@ -7,7 +7,7 @@ struct __bmc_aux_t; typedef struct { float het_rate, theta; - int n_hap, cap_mapQ; + int n_hap, cap_mapQ, is_soap; float eta, q_r; double *fk, *coef; diff --git a/samtools/bam_md.c b/samtools/bam_md.c index 8d07487..3ca7309 100644 --- a/samtools/bam_md.c +++ b/samtools/bam_md.c @@ -132,8 +132,11 @@ int bam_fillmd(int argc, char *argv[]) free(ref); ref = fai_fetch(fai, fp->header->target_name[b->core.tid], &len); tid = b->core.tid; + if (ref == 0) + fprintf(stderr, "[bam_fillmd] fail to find sequence '%s' in the reference.\n", + fp->header->target_name[tid]); } - bam_fillmd1(b, ref, is_equal); + if (ref) bam_fillmd1(b, ref, is_equal); } samwrite(fpout, b); } diff --git a/samtools/bam_plcmd.c b/samtools/bam_plcmd.c index 5bf1ed0..ba787a9 100644 --- a/samtools/bam_plcmd.c +++ b/samtools/bam_plcmd.c @@ -121,9 +121,11 @@ static int glt3_func(uint32_t tid, uint32_t pos, int n, const bam_pileup1_t *pu, g3->offset = pos - d->last_pos; d->last_pos = pos; glf3_write1(d->fp_glf, g3); - if (proposed_indels) - r = bam_maqindel(n, pos, d->ido, pu, d->ref, proposed_indels[0], proposed_indels+1); - else r = bam_maqindel(n, pos, d->ido, pu, d->ref, 0, 0); + if (pos < d->len) { + if (proposed_indels) + r = bam_maqindel(n, pos, d->ido, pu, d->ref, proposed_indels[0], proposed_indels+1); + else r = bam_maqindel(n, pos, d->ido, pu, d->ref, 0, 0); + } if (r) { // then write indel line int het = 3 * n, min; min = het; @@ -182,7 +184,7 @@ static int pileup_func(uint32_t tid, uint32_t pos, int n, const bam_pileup1_t *p // call the consensus and indel if (d->format & BAM_PLF_CNS) // call consensus cns = bam_maqcns_call(n, pu, d->c); - if ((d->format & (BAM_PLF_CNS|BAM_PLF_INDEL_ONLY)) && d->ref) { // call indels + if ((d->format & (BAM_PLF_CNS|BAM_PLF_INDEL_ONLY)) && d->ref && pos < d->len) { // call indels if (proposed_indels) // the first element gives the size of the array r = bam_maqindel(n, pos, d->ido, pu, d->ref, proposed_indels[0], proposed_indels+1); else r = bam_maqindel(n, pos, d->ido, pu, d->ref, 0, 0); @@ -299,8 +301,9 @@ int bam_pileup(int argc, char *argv[]) d->tid = -1; d->mask = BAM_DEF_MASK; d->c = bam_maqcns_init(); d->ido = bam_maqindel_opt_init(); - while ((c = getopt(argc, argv, "st:f:cT:N:r:l:im:gI:G:vM:S2")) >= 0) { + while ((c = getopt(argc, argv, "st:f:cT:N:r:l:im:gI:G:vM:S2a")) >= 0) { switch (c) { + case 'a': d->c->is_soap = 1; break; case 's': d->format |= BAM_PLF_SIMPLE; break; case 't': fn_list = strdup(optarg); break; case 'l': fn_pos = strdup(optarg); break; @@ -327,6 +330,7 @@ int bam_pileup(int argc, char *argv[]) fprintf(stderr, "Usage: samtools pileup [options] |\n\n"); fprintf(stderr, "Option: -s simple (yet incomplete) pileup format\n"); fprintf(stderr, " -S the input is in SAM\n"); + fprintf(stderr, " -a use the SOAPsnp model for SNP calling\n"); fprintf(stderr, " -2 output the 2nd best call and quality\n"); fprintf(stderr, " -i only show lines/consensus with indels\n"); fprintf(stderr, " -m INT filtering reads with bits in INT [%d]\n", d->mask); diff --git a/samtools/bam_rmdup.c b/samtools/bam_rmdup.c index 5da9460..f0d2b5d 100644 --- a/samtools/bam_rmdup.c +++ b/samtools/bam_rmdup.c @@ -2,15 +2,23 @@ #include #include #include -#include "bam.h" +#include +#include "sam.h" typedef bam1_t *bam1_p; + #include "khash.h" KHASH_SET_INIT_STR(name) KHASH_MAP_INIT_INT64(pos, bam1_p) #define BUFFER_SIZE 0x40000 +typedef struct { + uint64_t n_checked, n_removed; + khash_t(pos) *best_hash; +} lib_aux_t; +KHASH_MAP_INIT_STR(lib, lib_aux_t) + typedef struct { int n, max; bam1_t **a; @@ -25,15 +33,14 @@ static inline void stack_insert(tmp_stack_t *stack, bam1_t *b) stack->a[stack->n++] = b; } -static inline void dump_best(tmp_stack_t *stack, khash_t(pos) *best_hash, bamFile out) +static inline void dump_best(tmp_stack_t *stack, samfile_t *out) { int i; for (i = 0; i != stack->n; ++i) { - bam_write1(out, stack->a[i]); + samwrite(out, stack->a[i]); bam_destroy1(stack->a[i]); } stack->n = 0; - if (kh_size(best_hash) > BUFFER_SIZE) kh_clear(pos, best_hash); } static void clear_del_set(khash_t(name) *del_set) @@ -45,101 +52,155 @@ static void clear_del_set(khash_t(name) *del_set) kh_clear(name, del_set); } -void bam_rmdup_core(bamFile in, bamFile out) +static lib_aux_t *get_aux(khash_t(lib) *aux, const char *lib) +{ + khint_t k = kh_get(lib, aux, lib); + if (k == kh_end(aux)) { + int ret; + char *p = strdup(lib); + lib_aux_t *q; + k = kh_put(lib, aux, p, &ret); + q = &kh_val(aux, k); + q->n_checked = q->n_removed = 0; + q->best_hash = kh_init(pos); + return q; + } else return &kh_val(aux, k); +} + +static void clear_best(khash_t(lib) *aux, int max) +{ + khint_t k; + for (k = kh_begin(aux); k != kh_end(aux); ++k) { + if (kh_exist(aux, k)) { + lib_aux_t *q = &kh_val(aux, k); + if (kh_size(q->best_hash) >= max) + kh_clear(pos, q->best_hash); + } + } +} + +static inline int sum_qual(const bam1_t *b) +{ + int i, q; + uint8_t *qual = bam1_qual(b); + for (i = q = 0; i < b->core.l_qseq; ++i) q += qual[i]; + return q; +} + +void bam_rmdup_core(samfile_t *in, samfile_t *out) { - bam_header_t *header; bam1_t *b; int last_tid = -1, last_pos = -1; - uint64_t n_checked = 0, n_removed = 0; tmp_stack_t stack; khint_t k; - khash_t(pos) *best_hash; + khash_t(lib) *aux; khash_t(name) *del_set; - - best_hash = kh_init(pos); + + aux = kh_init(lib); del_set = kh_init(name); b = bam_init1(); memset(&stack, 0, sizeof(tmp_stack_t)); - header = bam_header_read(in); - bam_header_write(out, header); kh_resize(name, del_set, 4 * BUFFER_SIZE); - kh_resize(pos, best_hash, 3 * BUFFER_SIZE); - while (bam_read1(in, b) >= 0) { + while (samread(in, b) >= 0) { bam1_core_t *c = &b->core; if (c->tid != last_tid || last_pos != c->pos) { - dump_best(&stack, best_hash, out); // write the result + dump_best(&stack, out); // write the result + clear_best(aux, BUFFER_SIZE); if (c->tid != last_tid) { - kh_clear(pos, best_hash); + clear_best(aux, 0); if (kh_size(del_set)) { // check fprintf(stderr, "[bam_rmdup_core] %llu unmatched pairs\n", (long long)kh_size(del_set)); clear_del_set(del_set); } if ((int)c->tid == -1) { // append unmapped reads - bam_write1(out, b); - while (bam_read1(in, b) >= 0) bam_write1(out, b); + samwrite(out, b); + while (samread(in, b) >= 0) samwrite(out, b); break; } last_tid = c->tid; - fprintf(stderr, "[bam_rmdup_core] processing reference %s...\n", header->target_name[c->tid]); + fprintf(stderr, "[bam_rmdup_core] processing reference %s...\n", in->header->target_name[c->tid]); } } if (!(c->flag&BAM_FPAIRED) || (c->flag&(BAM_FUNMAP|BAM_FMUNMAP)) || (c->mtid >= 0 && c->tid != c->mtid)) { - bam_write1(out, b); + samwrite(out, b); } else if (c->isize > 0) { // paired, head uint64_t key = (uint64_t)c->pos<<32 | c->isize; + const char *lib; + lib_aux_t *q; int ret; - ++n_checked; - k = kh_put(pos, best_hash, key, &ret); + lib = bam_get_library(in->header, b); + q = lib? get_aux(aux, lib) : get_aux(aux, "\t"); + ++q->n_checked; + k = kh_put(pos, q->best_hash, key, &ret); if (ret == 0) { // found in best_hash - bam1_t *p = kh_val(best_hash, k); - ++n_removed; - if (p->core.qual < c->qual) { // the current alignment is better + bam1_t *p = kh_val(q->best_hash, k); + ++q->n_removed; + if (sum_qual(p) < sum_qual(b)) { // the current alignment is better; this can be accelerated in principle kh_put(name, del_set, strdup(bam1_qname(p)), &ret); // p will be removed bam_copy1(p, b); // replaced as b } else kh_put(name, del_set, strdup(bam1_qname(b)), &ret); // b will be removed if (ret == 0) fprintf(stderr, "[bam_rmdup_core] inconsistent BAM file for pair '%s'. Continue anyway.\n", bam1_qname(b)); } else { // not found in best_hash - kh_val(best_hash, k) = bam_dup1(b); - stack_insert(&stack, kh_val(best_hash, k)); + kh_val(q->best_hash, k) = bam_dup1(b); + stack_insert(&stack, kh_val(q->best_hash, k)); } } else { // paired, tail k = kh_get(name, del_set, bam1_qname(b)); if (k != kh_end(del_set)) { free((char*)kh_key(del_set, k)); kh_del(name, del_set, k); - } else bam_write1(out, b); + } else samwrite(out, b); } last_pos = c->pos; } - dump_best(&stack, best_hash, out); - bam_header_destroy(header); + for (k = kh_begin(aux); k != kh_end(aux); ++k) { + if (kh_exist(aux, k)) { + lib_aux_t *q = &kh_val(aux, k); + dump_best(&stack, out); + fprintf(stderr, "[bam_rmdup_core] %lld / %lld = %.4lf in library '%s'\n", (long long)q->n_removed, + (long long)q->n_checked, (double)q->n_removed/q->n_checked, kh_key(aux, k)); + kh_destroy(pos, q->best_hash); + free((char*)kh_key(aux, k)); + } + } + kh_destroy(lib, aux); + clear_del_set(del_set); kh_destroy(name, del_set); - kh_destroy(pos, best_hash); free(stack.a); bam_destroy1(b); - fprintf(stderr, "[bam_rmdup_core] %lld / %lld = %.4lf\n", (long long)n_removed, (long long)n_checked, - (double)n_removed/n_checked); } + +void bam_rmdupse_core(samfile_t *in, samfile_t *out, int force_se); + int bam_rmdup(int argc, char *argv[]) { - bamFile in, out; - if (argc < 3) { - fprintf(stderr, "Usage: samtools rmdup \n\n"); - fprintf(stderr, "Note: Picard is recommended for this task.\n"); + int c, is_se = 0, force_se = 0; + samfile_t *in, *out; + while ((c = getopt(argc, argv, "sS")) >= 0) { + switch (c) { + case 's': is_se = 1; break; + case 'S': force_se = is_se = 1; break; + } + } + if (optind + 2 > argc) { + fprintf(stderr, "\n"); + fprintf(stderr, "Usage: samtools rmdup [-sS] \n\n"); + fprintf(stderr, "Option: -s rmdup for SE reads\n"); + fprintf(stderr, " -S treat PE reads as SE in rmdup (force -s)\n\n"); return 1; } - in = (strcmp(argv[1], "-") == 0)? bam_dopen(fileno(stdin), "r") : bam_open(argv[1], "r"); - out = (strcmp(argv[2], "-") == 0)? bam_dopen(fileno(stdout), "w") : bam_open(argv[2], "w"); + in = samopen(argv[optind], "rb", 0); + out = samopen(argv[optind+1], "wb", in->header); if (in == 0 || out == 0) { fprintf(stderr, "[bam_rmdup] fail to read/write input files\n"); return 1; } - bam_rmdup_core(in, out); - bam_close(in); - bam_close(out); + if (is_se) bam_rmdupse_core(in, out, force_se); + else bam_rmdup_core(in, out); + samclose(in); samclose(out); return 0; } diff --git a/samtools/bam_rmdupse.c b/samtools/bam_rmdupse.c index cf1b7bd..e7dbdc7 100644 --- a/samtools/bam_rmdupse.c +++ b/samtools/bam_rmdupse.c @@ -1,178 +1,159 @@ #include #include "sam.h" #include "khash.h" +#include "klist.h" -typedef struct { - int n, m; - int *a; -} listelem_t; - -KHASH_MAP_INIT_INT(32, listelem_t) - -#define BLOCK_SIZE 65536 +#define QUEUE_CLEAR_SIZE 0x100000 +#define MAX_POS 0x7fffffff typedef struct { + int endpos; + uint32_t score:31, discarded:1; bam1_t *b; - int rpos, score; -} elem_t; +} elem_t, *elem_p; +#define __free_elem(p) bam_destroy1((p)->data.b) +KLIST_INIT(q, elem_t, __free_elem) +typedef klist_t(q) queue_t; + +KHASH_MAP_INIT_INT(best, elem_p) +typedef khash_t(best) besthash_t; typedef struct { - int n, max, x; - elem_t *buf; -} buffer_t; + uint64_t n_checked, n_removed; + besthash_t *left, *rght; +} lib_aux_t; +KHASH_MAP_INIT_STR(lib, lib_aux_t) -static int fill_buf(samfile_t *in, buffer_t *buf) +static lib_aux_t *get_aux(khash_t(lib) *aux, const char *lib) { - int i, ret, last_tid, min_rpos = 0x7fffffff, capacity; - bam1_t *b = bam_init1(); - bam1_core_t *c = &b->core; - // squeeze out the empty cells at the beginning - for (i = 0; i < buf->n; ++i) - if (buf->buf[i].b) break; - if (i < buf->n) { // squeeze - if (i > 0) { - memmove(buf->buf, buf->buf + i, sizeof(elem_t) * (buf->n - i)); - buf->n = buf->n - i; - } - } else buf->n = 0; - // calculate min_rpos - for (i = 0; i < buf->n; ++i) { - elem_t *e = buf->buf + i; - if (e->b && e->rpos >= 0 && e->rpos < min_rpos) - min_rpos = buf->buf[i].rpos; - } - // fill the buffer - buf->x = -1; - last_tid = buf->n? buf->buf[0].b->core.tid : -1; - capacity = buf->n + BLOCK_SIZE; - while ((ret = samread(in, b)) >= 0) { - elem_t *e; - uint8_t *qual = bam1_qual(b); - int is_mapped; - if (last_tid < 0) last_tid = c->tid; - if (c->tid != last_tid) { - if (buf->x < 0) buf->x = buf->n; - } - if (buf->n >= buf->max) { // enlarge - buf->max = buf->max? buf->max<<1 : 8; - buf->buf = (elem_t*)realloc(buf->buf, sizeof(elem_t) * buf->max); - } - e = &buf->buf[buf->n++]; - e->b = bam_dup1(b); - e->rpos = -1; e->score = 0; - for (i = 0; i < c->l_qseq; ++i) e->score += qual[i] + 1; - e->score = (double)e->score / sqrt(c->l_qseq + 1); - is_mapped = (c->tid < 0 || c->tid >= in->header->n_targets || (c->flag&BAM_FUNMAP))? 0 : 1; - if (!is_mapped) e->score = -1; - if (is_mapped && (c->flag & BAM_FREVERSE)) { - e->rpos = b->core.pos + bam_calend(&b->core, bam1_cigar(b)); - if (min_rpos > e->rpos) min_rpos = e->rpos; - } - if (buf->n >= capacity) { - if (is_mapped && c->pos <= min_rpos) capacity += BLOCK_SIZE; - else break; - } - } - if (ret >= 0 && buf->x < 0) buf->x = buf->n; - bam_destroy1(b); - return buf->n; + khint_t k = kh_get(lib, aux, lib); + if (k == kh_end(aux)) { + int ret; + char *p = strdup(lib); + lib_aux_t *q; + k = kh_put(lib, aux, p, &ret); + q = &kh_val(aux, k); + q->left = kh_init(best); + q->rght = kh_init(best); + q->n_checked = q->n_removed = 0; + return q; + } else return &kh_val(aux, k); +} + +static inline int sum_qual(const bam1_t *b) +{ + int i, q; + uint8_t *qual = bam1_qual(b); + for (i = q = 0; i < b->core.l_qseq; ++i) q += qual[i]; + return q; +} + +static inline elem_t *push_queue(queue_t *queue, const bam1_t *b, int endpos, int score) +{ + elem_t *p = kl_pushp(q, queue); + p->discarded = 0; + p->endpos = endpos; p->score = score; + if (p->b == 0) p->b = bam_init1(); + bam_copy1(p->b, b); + return p; } -static void rmdupse_buf(buffer_t *buf) +static void clear_besthash(besthash_t *h, int32_t pos) { - khash_t(32) *h; - uint32_t key; khint_t k; - int mpos, i, upper; - listelem_t *p; - mpos = 0x7fffffff; - mpos = (buf->x == buf->n)? buf->buf[buf->x-1].b->core.pos : 0x7fffffff; - upper = (buf->x < 0)? buf->n : buf->x; - // fill the hash table - h = kh_init(32); - for (i = 0; i < upper; ++i) { - elem_t *e = buf->buf + i; - int ret; - if (e->score < 0) continue; - if (e->rpos >= 0) { - if (e->rpos <= mpos) key = (uint32_t)e->rpos<<1 | 1; - else continue; - } else { - if (e->b->core.pos < mpos) key = (uint32_t)e->b->core.pos<<1; - else continue; - } - k = kh_put(32, h, key, &ret); - p = &kh_val(h, k); - if (ret == 0) { // present in the hash table - if (p->n == p->m) { - p->m <<= 1; - p->a = (int*)realloc(p->a, p->m * sizeof(int)); + for (k = kh_begin(h); k != kh_end(h); ++k) + if (kh_exist(h, k) && kh_val(h, k)->endpos <= pos) + kh_del(best, h, k); +} + +static void dump_alignment(samfile_t *out, queue_t *queue, int32_t pos, khash_t(lib) *h) +{ + if (queue->size > QUEUE_CLEAR_SIZE || pos == MAX_POS) { + khint_t k; + while (1) { + elem_t *q; + if (queue->head == queue->tail) break; + q = &kl_val(queue->head); + if (q->discarded) { + q->b->data_len = 0; + kl_shift(q, queue, 0); + continue; } - p->a[p->n++] = i; - } else { - p->m = p->n = 1; - p->a = (int*)calloc(p->m, sizeof(int)); - p->a[0] = i; + if ((q->b->core.flag&BAM_FREVERSE) && q->endpos > pos) break; + samwrite(out, q->b); + q->b->data_len = 0; + kl_shift(q, queue, 0); } - } - // rmdup - for (k = kh_begin(h); k < kh_end(h); ++k) { - if (kh_exist(h, k)) { - int max, maxi; - p = &kh_val(h, k); - // get the max - for (i = max = 0, maxi = -1; i < p->n; ++i) { - if (buf->buf[p->a[i]].score > max) { - max = buf->buf[p->a[i]].score; - maxi = i; - } - } - // mark the elements - for (i = 0; i < p->n; ++i) { - buf->buf[p->a[i]].score = -1; - if (i != maxi) { - bam_destroy1(buf->buf[p->a[i]].b); - buf->buf[p->a[i]].b = 0; - } + for (k = kh_begin(h); k != kh_end(h); ++k) { + if (kh_exist(h, k)) { + clear_besthash(kh_val(h, k).left, pos); + clear_besthash(kh_val(h, k).rght, pos); } - // free - free(p->a); } } - kh_destroy(32, h); } -static void dump_buf(buffer_t *buf, samfile_t *out) +void bam_rmdupse_core(samfile_t *in, samfile_t *out, int force_se) { - int i; - for (i = 0; i < buf->n; ++i) { - elem_t *e = buf->buf + i; - if (e->score != -1) break; - if (e->b) { - samwrite(out, e->b); - bam_destroy1(e->b); - e->b = 0; + bam1_t *b; + queue_t *queue; + khint_t k; + int last_tid = -2; + khash_t(lib) *aux; + + aux = kh_init(lib); + b = bam_init1(); + queue = kl_init(q); + while (samread(in, b) >= 0) { + bam1_core_t *c = &b->core; + int endpos = bam_calend(c, bam1_cigar(b)); + int score = sum_qual(b); + + if (last_tid != c->tid) { + if (last_tid >= 0) dump_alignment(out, queue, MAX_POS, aux); + last_tid = c->tid; + } else dump_alignment(out, queue, c->pos, aux); + if ((c->flag&BAM_FUNMAP) || ((c->flag&BAM_FPAIRED) && !force_se)) { + push_queue(queue, b, endpos, score); + } else { + const char *lib; + lib_aux_t *q; + besthash_t *h; + uint32_t key; + int ret; + lib = bam_get_library(in->header, b); + q = lib? get_aux(aux, lib) : get_aux(aux, "\t"); + ++q->n_checked; + h = (c->flag&BAM_FREVERSE)? q->rght : q->left; + key = (c->flag&BAM_FREVERSE)? endpos : c->pos; + k = kh_put(best, h, key, &ret); + if (ret == 0) { // in the hash table + elem_t *p = kh_val(h, k); + ++q->n_removed; + if (p->score < score) { + if (c->flag&BAM_FREVERSE) { // mark "discarded" and push the queue + p->discarded = 1; + kh_val(h, k) = push_queue(queue, b, endpos, score); + } else { // replace + p->score = score; p->endpos = endpos; + bam_copy1(p->b, b); + } + } // otherwise, discard the alignment + } else kh_val(h, k) = push_queue(queue, b, endpos, score); } } -} + dump_alignment(out, queue, MAX_POS, aux); -int bam_rmdupse(int argc, char *argv[]) -{ - samfile_t *in, *out; - buffer_t *buf; - if (argc < 3) { - fprintf(stderr, "Usage: samtools rmdupse \n\n"); - fprintf(stderr, "Note: Picard is recommended for this task.\n"); - return 1; - } - buf = calloc(1, sizeof(buffer_t)); - in = samopen(argv[1], "rb", 0); - out = samopen(argv[2], "wb", in->header); - while (fill_buf(in, buf)) { - rmdupse_buf(buf); - dump_buf(buf, out); + for (k = kh_begin(aux); k != kh_end(aux); ++k) { + if (kh_exist(aux, k)) { + lib_aux_t *q = &kh_val(aux, k); + fprintf(stderr, "[bam_rmdupse_core] %lld / %lld = %.4lf in library '%s'\n", (long long)q->n_removed, + (long long)q->n_checked, (double)q->n_removed/q->n_checked, kh_key(aux, k)); + kh_destroy(best, q->left); kh_destroy(best, q->rght); + free((char*)kh_key(aux, k)); + } } - samclose(in); samclose(out); - free(buf->buf); free(buf); - return 0; + kh_destroy(lib, aux); + bam_destroy1(b); + kl_destroy(q, queue); } diff --git a/samtools/bam_sort.c b/samtools/bam_sort.c index 651f7f9..9884f3d 100644 --- a/samtools/bam_sort.c +++ b/samtools/bam_sort.c @@ -33,18 +33,20 @@ static inline int strnum_cmp(const char *a, const char *b) typedef struct { int i; - uint64_t pos; + uint64_t pos, idx; bam1_t *b; } heap1_t; +#define __pos_cmp(a, b) ((a).pos > (b).pos || ((a).pos == (b).pos && ((a).i > (b).i || ((a).i == (b).i && (a).idx > (b).idx)))) + static inline int heap_lt(const heap1_t a, const heap1_t b) { if (g_is_by_qname) { int t; if (a.b == 0 || b.b == 0) return a.b == 0? 1 : 0; t = strnum_cmp(bam1_qname(a.b), bam1_qname(b.b)); - return (t > 0 || (t == 0 && a.pos > b.pos)); - } else return (a.pos > b.pos); + return (t > 0 || (t == 0 && __pos_cmp(a, b))); + } else return __pos_cmp(a, b); } KSORT_INIT(heap, heap1_t, heap_lt) @@ -69,13 +71,15 @@ static void swap_header_text(bam_header_t *h1, bam_header_t *h2) @discussion Padding information may NOT correctly maintained. This function is NOT thread safe. */ -void bam_merge_core(int by_qname, const char *out, const char *headers, int n, char * const *fn) +void bam_merge_core(int by_qname, const char *out, const char *headers, int n, char * const *fn, int add_RG) { bamFile fpout, *fp; heap1_t *heap; bam_header_t *hout = 0; bam_header_t *hheaders = NULL; - int i, j; + int i, j, *RG_len = 0; + uint64_t idx = 0; + char **RG = 0; if (headers) { tamFile fpheaders = sam_open(headers); @@ -90,11 +94,34 @@ void bam_merge_core(int by_qname, const char *out, const char *headers, int n, c g_is_by_qname = by_qname; fp = (bamFile*)calloc(n, sizeof(bamFile)); heap = (heap1_t*)calloc(n, sizeof(heap1_t)); + // prepare RG tag + if (add_RG) { + RG = (char**)calloc(n, sizeof(void*)); + RG_len = (int*)calloc(n, sizeof(int)); + for (i = 0; i != n; ++i) { + int l = strlen(fn[i]); + const char *s = fn[i]; + if (l > 4 && strcmp(s + l - 4, ".bam") == 0) l -= 4; + for (j = l - 1; j >= 0; --j) if (s[j] == '/') break; + ++j; l -= j; + RG[i] = calloc(l + 1, 1); + RG_len[i] = l; + strncpy(RG[i], s + j, l); + } + } + // read the first for (i = 0; i != n; ++i) { heap1_t *h; bam_header_t *hin; fp[i] = bam_open(fn[i], "r"); - assert(fp[i]); + if (fp[i] == 0) { + int j; + fprintf(stderr, "[bam_merge_core] fail to open file %s\n", fn[i]); + for (j = 0; j < i; ++j) bam_close(fp[j]); + free(fp); free(heap); + // FIXME: possible memory leak + return; + } hin = bam_header_read(fp[i]); if (i == 0) { // the first SAM hout = hin; @@ -130,8 +157,10 @@ void bam_merge_core(int by_qname, const char *out, const char *headers, int n, c h = heap + i; h->i = i; h->b = (bam1_t*)calloc(1, sizeof(bam1_t)); - if (bam_read1(fp[i], h->b) >= 0) + if (bam_read1(fp[i], h->b) >= 0) { h->pos = ((uint64_t)h->b->core.tid<<32) | (uint32_t)h->b->core.pos<<1 | bam1_strand(h->b); + h->idx = idx++; + } else h->pos = HEAP_EMPTY; } fpout = strcmp(out, "-")? bam_open(out, "w") : bam_dopen(fileno(stdout), "w"); @@ -142,9 +171,12 @@ void bam_merge_core(int by_qname, const char *out, const char *headers, int n, c ks_heapmake(heap, n, heap); while (heap->pos != HEAP_EMPTY) { bam1_t *b = heap->b; + if (add_RG && bam_aux_get(b, "RG") == 0) + bam_aux_append(b, "RG", 'Z', RG_len[heap->i] + 1, (uint8_t*)RG[heap->i]); bam_write1_core(fpout, &b->core, b->data_len, b->data); if ((j = bam_read1(fp[heap->i], b)) >= 0) { heap->pos = ((uint64_t)b->core.tid<<32) | (uint32_t)b->core.pos<<1 | bam1_strand(b); + heap->idx = idx++; } else if (j == -1) { heap->pos = HEAP_EMPTY; free(heap->b->data); free(heap->b); @@ -153,32 +185,38 @@ void bam_merge_core(int by_qname, const char *out, const char *headers, int n, c ks_heapadjust(heap, 0, n, heap); } + if (add_RG) { + for (i = 0; i != n; ++i) free(RG[i]); + free(RG); free(RG_len); + } for (i = 0; i != n; ++i) bam_close(fp[i]); bam_close(fpout); free(fp); free(heap); } int bam_merge(int argc, char *argv[]) { - int c, is_by_qname = 0; + int c, is_by_qname = 0, add_RG = 0; char *fn_headers = NULL; - while ((c = getopt(argc, argv, "h:n")) >= 0) { + while ((c = getopt(argc, argv, "h:nr")) >= 0) { switch (c) { + case 'r': add_RG = 1; break; case 'h': fn_headers = strdup(optarg); break; case 'n': is_by_qname = 1; break; } } if (optind + 2 >= argc) { fprintf(stderr, "\n"); - fprintf(stderr, "Usage: samtools merge [-n] [-h inh.sam] [...]\n\n"); + fprintf(stderr, "Usage: samtools merge [-nr] [-h inh.sam] [...]\n\n"); fprintf(stderr, "Options: -n sort by read names\n"); + fprintf(stderr, " -r attach RG tag (inferred from file names)\n"); fprintf(stderr, " -h FILE copy the header in FILE to [in1.bam]\n\n"); fprintf(stderr, "Note: Samtools' merge does not reconstruct the @RG dictionary in the header. Users\n"); fprintf(stderr, " must provide the correct header with -h, or uses Picard which properly maintains\n"); fprintf(stderr, " the header dictionary in merging.\n\n"); return 1; } - bam_merge_core(is_by_qname, argv[optind], fn_headers, argc - optind - 1, argv + optind + 1); + bam_merge_core(is_by_qname, argv[optind], fn_headers, argc - optind - 1, argv + optind + 1, add_RG); free(fn_headers); return 0; } @@ -194,7 +232,7 @@ static inline int bam1_lt(const bam1_p a, const bam1_p b) } KSORT_INIT(sort, bam1_p, bam1_lt) -static void sort_blocks(int n, int k, bam1_p *buf, const char *prefix, const bam_header_t *h) +static void sort_blocks(int n, int k, bam1_p *buf, const char *prefix, const bam_header_t *h, int is_stdout) { char *name; int i; @@ -203,8 +241,13 @@ static void sort_blocks(int n, int k, bam1_p *buf, const char *prefix, const bam name = (char*)calloc(strlen(prefix) + 20, 1); if (n >= 0) sprintf(name, "%s.%.4d.bam", prefix, n); else sprintf(name, "%s.bam", prefix); - fp = bam_open(name, "w"); - assert(fp); + fp = is_stdout? bam_dopen(fileno(stdout), "w") : bam_open(name, "w"); + if (fp == 0) { + fprintf(stderr, "[sort_blocks] fail to create file %s.\n", name); + free(name); + // FIXME: possible memory leak + return; + } free(name); bam_header_write(fp, h); for (i = 0; i < k; ++i) @@ -226,7 +269,7 @@ static void sort_blocks(int n, int k, bam1_p *buf, const char *prefix, const bam and then merge them by calling bam_merge_core(). This function is NOT thread safe. */ -void bam_sort_core(int is_by_qname, const char *fn, const char *prefix, size_t max_mem) +void bam_sort_core_ext(int is_by_qname, const char *fn, const char *prefix, size_t max_mem, int is_stdout) { int n, ret, k, i; size_t mem; @@ -237,7 +280,10 @@ void bam_sort_core(int is_by_qname, const char *fn, const char *prefix, size_t m g_is_by_qname = is_by_qname; n = k = 0; mem = 0; fp = strcmp(fn, "-")? bam_open(fn, "r") : bam_dopen(fileno(stdin), "r"); - assert(fp); + if (fp == 0) { + fprintf(stderr, "[bam_sort_core] fail to open file %s\n", fn); + return; + } header = bam_header_read(fp); buf = (bam1_t**)calloc(max_mem / BAM_CORE_SIZE, sizeof(bam1_t*)); // write sub files @@ -248,25 +294,26 @@ void bam_sort_core(int is_by_qname, const char *fn, const char *prefix, size_t m mem += ret; ++k; if (mem >= max_mem) { - sort_blocks(n++, k, buf, prefix, header); + sort_blocks(n++, k, buf, prefix, header, is_stdout); mem = 0; k = 0; } } if (ret != -1) fprintf(stderr, "[bam_sort_core] truncated file. Continue anyway.\n"); - if (n == 0) sort_blocks(-1, k, buf, prefix, header); + if (n == 0) sort_blocks(-1, k, buf, prefix, header, is_stdout); else { // then merge char **fns, *fnout; fprintf(stderr, "[bam_sort_core] merging from %d files...\n", n+1); - sort_blocks(n++, k, buf, prefix, header); + sort_blocks(n++, k, buf, prefix, header, is_stdout); fnout = (char*)calloc(strlen(prefix) + 20, 1); - sprintf(fnout, "%s.bam", prefix); + if (is_stdout) sprintf(fnout, "-"); + else sprintf(fnout, "%s.bam", prefix); fns = (char**)calloc(n, sizeof(char*)); for (i = 0; i < n; ++i) { fns[i] = (char*)calloc(strlen(prefix) + 20, 1); sprintf(fns[i], "%s.%.4d.bam", prefix, i); } - bam_merge_core(is_by_qname, fnout, 0, n, fns); + bam_merge_core(is_by_qname, fnout, 0, n, fns, 0); free(fnout); for (i = 0; i < n; ++i) { unlink(fns[i]); @@ -285,22 +332,26 @@ void bam_sort_core(int is_by_qname, const char *fn, const char *prefix, size_t m bam_close(fp); } +void bam_sort_core(int is_by_qname, const char *fn, const char *prefix, size_t max_mem) +{ + bam_sort_core_ext(is_by_qname, fn, prefix, max_mem, 0); +} + int bam_sort(int argc, char *argv[]) { size_t max_mem = 500000000; - int c, is_by_qname = 0; - - while ((c = getopt(argc, argv, "nm:")) >= 0) { - switch (c) { - case 'n': is_by_qname = 1; break; - case 'm': max_mem = atol(optarg); break; - } + int c, is_by_qname = 0, is_stdout = 0; + while ((c = getopt(argc, argv, "nom:")) >= 0) { + switch (c) { + case 'o': is_stdout = 1; break; + case 'n': is_by_qname = 1; break; + case 'm': max_mem = atol(optarg); break; + } } if (optind + 2 > argc) { - fprintf(stderr, "Usage: samtools sort [-n] [-m ] \n"); + fprintf(stderr, "Usage: samtools sort [-on] [-m ] \n"); return 1; } - - bam_sort_core(is_by_qname, argv[optind], argv[optind+1], max_mem); + bam_sort_core_ext(is_by_qname, argv[optind], argv[optind+1], max_mem, is_stdout); return 0; } diff --git a/samtools/bgzf.c b/samtools/bgzf.c index 646b2b4..59f902f 100644 --- a/samtools/bgzf.c +++ b/samtools/bgzf.c @@ -199,7 +199,7 @@ bgzf_open(const char* __restrict path, const char* __restrict mode) #ifdef _WIN32 oflag |= O_BINARY; #endif - fd = open(path, oflag, 0644); + fd = open(path, oflag, 0666); if (fd == -1) return 0; fp = open_write(fd, strstr(mode, "u")? 1 : 0); } diff --git a/samtools/faidx.c b/samtools/faidx.c index 055445f..811bdf8 100644 --- a/samtools/faidx.c +++ b/samtools/faidx.c @@ -28,6 +28,9 @@ extern int fseeko(FILE *stream, off_t offset, int whence); #define razf_seek(fp, offset, whence) fseeko(fp, offset, whence) #define razf_tell(fp) ftello(fp) #endif +#ifdef _USE_KNETFILE +#include "knetfile.h" +#endif struct __faidx_t { RAZF *rz; @@ -194,7 +197,7 @@ int fai_build(const char *fn) sprintf(str, "%s.fai", fn); rz = razf_open(fn, "r"); if (rz == 0) { - fprintf(stderr, "[fai_build] fail to open the FASTA file.\n"); + fprintf(stderr, "[fai_build] fail to open the FASTA file %s\n",str); free(str); return -1; } @@ -202,7 +205,7 @@ int fai_build(const char *fn) razf_close(rz); fp = fopen(str, "wb"); if (fp == 0) { - fprintf(stderr, "[fai_build] fail to write FASTA index.\n"); + fprintf(stderr, "[fai_build] fail to write FASTA index %s\n",str); fai_destroy(fai); free(str); return -1; } @@ -213,6 +216,47 @@ int fai_build(const char *fn) return 0; } +#ifdef _USE_KNETFILE +FILE *download_and_open(const char *fn) +{ + const int buf_size = 1 * 1024 * 1024; + uint8_t *buf; + FILE *fp; + knetFile *fp_remote; + const char *url = fn; + const char *p; + int l = strlen(fn); + for (p = fn + l - 1; p >= fn; --p) + if (*p == '/') break; + fn = p + 1; + + // First try to open a local copy + fp = fopen(fn, "r"); + if (fp) + return fp; + + // If failed, download from remote and open + fp_remote = knet_open(url, "rb"); + if (fp_remote == 0) { + fprintf(stderr, "[download_from_remote] fail to open remote file %s\n",url); + return NULL; + } + if ((fp = fopen(fn, "wb")) == 0) { + fprintf(stderr, "[download_from_remote] fail to create file in the working directory %s\n",fn); + knet_close(fp_remote); + return NULL; + } + buf = (uint8_t*)calloc(buf_size, 1); + while ((l = knet_read(fp_remote, buf, buf_size)) != 0) + fwrite(buf, 1, l, fp); + free(buf); + fclose(fp); + knet_close(fp_remote); + + return fopen(fn, "r"); +} +#endif + faidx_t *fai_load(const char *fn) { char *str; @@ -220,19 +264,35 @@ faidx_t *fai_load(const char *fn) faidx_t *fai; str = (char*)calloc(strlen(fn) + 5, 1); sprintf(str, "%s.fai", fn); - fp = fopen(str, "rb"); + +#ifdef _USE_KNETFILE + if (strstr(fn, "ftp://") == fn || strstr(fn, "http://") == fn) + { + fp = download_and_open(str); + if ( !fp ) + { + fprintf(stderr, "[fai_load] failed to open remote FASTA index %s\n", str); + free(str); + return 0; + } + } + else +#endif + fp = fopen(str, "rb"); if (fp == 0) { fprintf(stderr, "[fai_load] build FASTA index.\n"); fai_build(fn); - fp = fopen(str, "r"); + fp = fopen(str, "rb"); if (fp == 0) { fprintf(stderr, "[fai_load] fail to open FASTA index.\n"); free(str); return 0; } } + fai = fai_read(fp); fclose(fp); + fai->rz = razf_open(fn, "rb"); free(str); if (fai->rz == 0) { @@ -287,7 +347,7 @@ char *fai_fetch(const faidx_t *fai, const char *str, int *len) l = 0; s = (char*)malloc(end - beg + 2); razf_seek(fai->rz, val.offset + beg / val.line_blen * val.line_len + beg % val.line_blen, SEEK_SET); - while (razf_read(fai->rz, &c, 1) == 1 && l < end - beg) + while (razf_read(fai->rz, &c, 1) == 1 && l < end - beg && !fai->rz->z_err) if (isgraph(c)) s[l++] = c; s[l] = '\0'; *len = l; @@ -323,6 +383,40 @@ int faidx_main(int argc, char *argv[]) return 0; } +int faidx_fetch_nseq(const faidx_t *fai) +{ + return fai->n; +} + +char *faidx_fetch_seq(const faidx_t *fai, char *c_name, int p_beg_i, int p_end_i, int *len) +{ + int l; + char c; + khiter_t iter; + faidx1_t val; + char *seq=NULL; + + // Adjust position + iter = kh_get(s, fai->hash, c_name); + if(iter == kh_end(fai->hash)) return 0; + val = kh_value(fai->hash, iter); + if(p_end_i < p_beg_i) p_beg_i = p_end_i; + if(p_beg_i < 0) p_beg_i = 0; + else if(val.len <= p_beg_i) p_beg_i = val.len - 1; + if(p_end_i < 0) p_end_i = 0; + else if(val.len <= p_end_i) p_end_i = val.len - 1; + + // Now retrieve the sequence + l = 0; + seq = (char*)malloc(p_end_i - p_beg_i + 2); + razf_seek(fai->rz, val.offset + p_beg_i / val.line_blen * val.line_len + p_beg_i % val.line_blen, SEEK_SET); + while (razf_read(fai->rz, &c, 1) == 1 && l < p_end_i - p_beg_i + 1) + if (isgraph(c)) seq[l++] = c; + seq[l] = '\0'; + *len = l; + return seq; +} + #ifdef FAIDX_MAIN int main(int argc, char *argv[]) { return faidx_main(argc, argv); } #endif diff --git a/samtools/faidx.h b/samtools/faidx.h index 1a52fb7..1fb1b1f 100644 --- a/samtools/faidx.h +++ b/samtools/faidx.h @@ -75,6 +75,27 @@ extern "C" { */ char *fai_fetch(const faidx_t *fai, const char *reg, int *len); + /*! + @abstract Fetch the number of sequences. + @param fai Pointer to the faidx_t struct + @return The number of sequences + */ + int faidx_fetch_nseq(const faidx_t *fai); + + /*! + @abstract Fetch the sequence in a region. + @param fai Pointer to the faidx_t struct + @param c_name Region name + @param p_beg_i Beginning position number (zero-based) + @param p_end_i End position number (zero-based) + @param len Length of the region + @return Pointer to the sequence; null on failure + + @discussion The returned sequence is allocated by malloc family + and should be destroyed by end users by calling free() on it. + */ + char *faidx_fetch_seq(const faidx_t *fai, char *c_name, int p_beg_i, int p_end_i, int *len); + #ifdef __cplusplus } #endif diff --git a/samtools/kaln.c b/samtools/kaln.c new file mode 100644 index 0000000..9fa40d0 --- /dev/null +++ b/samtools/kaln.c @@ -0,0 +1,370 @@ +/* The MIT License + + Copyright (c) 2003-2006, 2008, 2009, by Heng Li + + Permission is hereby granted, free of charge, to any person obtaining + a copy of this software and associated documentation files (the + "Software"), to deal in the Software without restriction, including + without limitation the rights to use, copy, modify, merge, publish, + distribute, sublicense, and/or sell copies of the Software, and to + permit persons to whom the Software is furnished to do so, subject to + the following conditions: + + The above copyright notice and this permission notice shall be + included in all copies or substantial portions of the Software. + + THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, + EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF + MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND + NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS + BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN + ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN + CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE + SOFTWARE. +*/ + +#include +#include +#include +#include +#include "kaln.h" + +#define FROM_M 0 +#define FROM_I 1 +#define FROM_D 2 + +typedef struct { + int i, j; + unsigned char ctype; +} path_t; + +int aln_sm_blosum62[] = { +/* A R N D C Q E G H I L K M F P S T W Y V * X */ + 4,-1,-2,-2, 0,-1,-1, 0,-2,-1,-1,-1,-1,-2,-1, 1, 0,-3,-2, 0,-4, 0, + -1, 5, 0,-2,-3, 1, 0,-2, 0,-3,-2, 2,-1,-3,-2,-1,-1,-3,-2,-3,-4,-1, + -2, 0, 6, 1,-3, 0, 0, 0, 1,-3,-3, 0,-2,-3,-2, 1, 0,-4,-2,-3,-4,-1, + -2,-2, 1, 6,-3, 0, 2,-1,-1,-3,-4,-1,-3,-3,-1, 0,-1,-4,-3,-3,-4,-1, + 0,-3,-3,-3, 9,-3,-4,-3,-3,-1,-1,-3,-1,-2,-3,-1,-1,-2,-2,-1,-4,-2, + -1, 1, 0, 0,-3, 5, 2,-2, 0,-3,-2, 1, 0,-3,-1, 0,-1,-2,-1,-2,-4,-1, + -1, 0, 0, 2,-4, 2, 5,-2, 0,-3,-3, 1,-2,-3,-1, 0,-1,-3,-2,-2,-4,-1, + 0,-2, 0,-1,-3,-2,-2, 6,-2,-4,-4,-2,-3,-3,-2, 0,-2,-2,-3,-3,-4,-1, + -2, 0, 1,-1,-3, 0, 0,-2, 8,-3,-3,-1,-2,-1,-2,-1,-2,-2, 2,-3,-4,-1, + -1,-3,-3,-3,-1,-3,-3,-4,-3, 4, 2,-3, 1, 0,-3,-2,-1,-3,-1, 3,-4,-1, + -1,-2,-3,-4,-1,-2,-3,-4,-3, 2, 4,-2, 2, 0,-3,-2,-1,-2,-1, 1,-4,-1, + -1, 2, 0,-1,-3, 1, 1,-2,-1,-3,-2, 5,-1,-3,-1, 0,-1,-3,-2,-2,-4,-1, + -1,-1,-2,-3,-1, 0,-2,-3,-2, 1, 2,-1, 5, 0,-2,-1,-1,-1,-1, 1,-4,-1, + -2,-3,-3,-3,-2,-3,-3,-3,-1, 0, 0,-3, 0, 6,-4,-2,-2, 1, 3,-1,-4,-1, + -1,-2,-2,-1,-3,-1,-1,-2,-2,-3,-3,-1,-2,-4, 7,-1,-1,-4,-3,-2,-4,-2, + 1,-1, 1, 0,-1, 0, 0, 0,-1,-2,-2, 0,-1,-2,-1, 4, 1,-3,-2,-2,-4, 0, + 0,-1, 0,-1,-1,-1,-1,-2,-2,-1,-1,-1,-1,-2,-1, 1, 5,-2,-2, 0,-4, 0, + -3,-3,-4,-4,-2,-2,-3,-2,-2,-3,-2,-3,-1, 1,-4,-3,-2,11, 2,-3,-4,-2, + -2,-2,-2,-3,-2,-1,-2,-3, 2,-1,-1,-2,-1, 3,-3,-2,-2, 2, 7,-1,-4,-1, + 0,-3,-3,-3,-1,-2,-2,-3,-3, 3, 1,-2, 1,-1,-2,-2, 0,-3,-1, 4,-4,-1, + -4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4, 1,-4, + 0,-1,-1,-1,-2,-1,-1,-1,-1,-1,-1,-1,-1,-1,-2, 0, 0,-2,-1,-1,-4,-1 +}; + +int aln_sm_blast[] = { + 1, -3, -3, -3, -2, + -3, 1, -3, -3, -2, + -3, -3, 1, -3, -2, + -3, -3, -3, 1, -2, + -2, -2, -2, -2, -2 +}; + +ka_param_t ka_param_blast = { 5, 2, 2, aln_sm_blast, 5, 50 }; +ka_param_t ka_param_aa2aa = { 10, 2, 2, aln_sm_blosum62, 22, 50 }; + +static uint32_t *ka_path2cigar32(const path_t *path, int path_len, int *n_cigar) +{ + int i, n; + uint32_t *cigar; + unsigned char last_type; + + if (path_len == 0 || path == 0) { + *n_cigar = 0; + return 0; + } + + last_type = path->ctype; + for (i = n = 1; i < path_len; ++i) { + if (last_type != path[i].ctype) ++n; + last_type = path[i].ctype; + } + *n_cigar = n; + cigar = (uint32_t*)calloc(*n_cigar, 4); + + cigar[0] = 1u << 4 | path[path_len-1].ctype; + last_type = path[path_len-1].ctype; + for (i = path_len - 2, n = 0; i >= 0; --i) { + if (path[i].ctype == last_type) cigar[n] += 1u << 4; + else { + cigar[++n] = 1u << 4 | path[i].ctype; + last_type = path[i].ctype; + } + } + + return cigar; +} + +/***************************/ +/* START OF common_align.c */ +/***************************/ + +#define SET_INF(s) (s).M = (s).I = (s).D = MINOR_INF; + +#define set_M(MM, cur, p, sc) \ +{ \ + if ((p)->M >= (p)->I) { \ + if ((p)->M >= (p)->D) { \ + (MM) = (p)->M + (sc); (cur)->Mt = FROM_M; \ + } else { \ + (MM) = (p)->D + (sc); (cur)->Mt = FROM_D; \ + } \ + } else { \ + if ((p)->I > (p)->D) { \ + (MM) = (p)->I + (sc); (cur)->Mt = FROM_I; \ + } else { \ + (MM) = (p)->D + (sc); (cur)->Mt = FROM_D; \ + } \ + } \ +} +#define set_I(II, cur, p) \ +{ \ + if ((p)->M - gap_open > (p)->I) { \ + (cur)->It = FROM_M; \ + (II) = (p)->M - gap_open - gap_ext; \ + } else { \ + (cur)->It = FROM_I; \ + (II) = (p)->I - gap_ext; \ + } \ +} +#define set_end_I(II, cur, p) \ +{ \ + if (gap_end >= 0) { \ + if ((p)->M - gap_open > (p)->I) { \ + (cur)->It = FROM_M; \ + (II) = (p)->M - gap_open - gap_end; \ + } else { \ + (cur)->It = FROM_I; \ + (II) = (p)->I - gap_end; \ + } \ + } else set_I(II, cur, p); \ +} +#define set_D(DD, cur, p) \ +{ \ + if ((p)->M - gap_open > (p)->D) { \ + (cur)->Dt = FROM_M; \ + (DD) = (p)->M - gap_open - gap_ext; \ + } else { \ + (cur)->Dt = FROM_D; \ + (DD) = (p)->D - gap_ext; \ + } \ +} +#define set_end_D(DD, cur, p) \ +{ \ + if (gap_end >= 0) { \ + if ((p)->M - gap_open > (p)->D) { \ + (cur)->Dt = FROM_M; \ + (DD) = (p)->M - gap_open - gap_end; \ + } else { \ + (cur)->Dt = FROM_D; \ + (DD) = (p)->D - gap_end; \ + } \ + } else set_D(DD, cur, p); \ +} + +typedef struct { + uint8_t Mt:3, It:2, Dt:2; +} dpcell_t; + +typedef struct { + int M, I, D; +} dpscore_t; + +/*************************** + * banded global alignment * + ***************************/ +uint32_t *ka_global_core(uint8_t *seq1, int len1, uint8_t *seq2, int len2, const ka_param_t *ap, int *_score, int *n_cigar) +{ + int i, j; + dpcell_t **dpcell, *q; + dpscore_t *curr, *last, *s; + int b1, b2, tmp_end; + int *mat, end, max = 0; + uint8_t type, ctype; + uint32_t *cigar = 0; + + int gap_open, gap_ext, gap_end, b; + int *score_matrix, N_MATRIX_ROW; + + /* initialize some align-related parameters. just for compatibility */ + gap_open = ap->gap_open; + gap_ext = ap->gap_ext; + gap_end = ap->gap_end; + b = ap->band_width; + score_matrix = ap->matrix; + N_MATRIX_ROW = ap->row; + + *n_cigar = 0; + if (len1 == 0 || len2 == 0) return 0; + + /* calculate b1 and b2 */ + if (len1 > len2) { + b1 = len1 - len2 + b; + b2 = b; + } else { + b1 = b; + b2 = len2 - len1 + b; + } + if (b1 > len1) b1 = len1; + if (b2 > len2) b2 = len2; + --seq1; --seq2; + + /* allocate memory */ + end = (b1 + b2 <= len1)? (b1 + b2 + 1) : (len1 + 1); + dpcell = (dpcell_t**)malloc(sizeof(dpcell_t*) * (len2 + 1)); + for (j = 0; j <= len2; ++j) + dpcell[j] = (dpcell_t*)malloc(sizeof(dpcell_t) * end); + for (j = b2 + 1; j <= len2; ++j) + dpcell[j] -= j - b2; + curr = (dpscore_t*)malloc(sizeof(dpscore_t) * (len1 + 1)); + last = (dpscore_t*)malloc(sizeof(dpscore_t) * (len1 + 1)); + + /* set first row */ + SET_INF(*curr); curr->M = 0; + for (i = 1, s = curr + 1; i < b1; ++i, ++s) { + SET_INF(*s); + set_end_D(s->D, dpcell[0] + i, s - 1); + } + s = curr; curr = last; last = s; + + /* core dynamic programming, part 1 */ + tmp_end = (b2 < len2)? b2 : len2 - 1; + for (j = 1; j <= tmp_end; ++j) { + q = dpcell[j]; s = curr; SET_INF(*s); + set_end_I(s->I, q, last); + end = (j + b1 <= len1 + 1)? (j + b1 - 1) : len1; + mat = score_matrix + seq2[j] * N_MATRIX_ROW; + ++s; ++q; + for (i = 1; i != end; ++i, ++s, ++q) { + set_M(s->M, q, last + i - 1, mat[seq1[i]]); /* this will change s->M ! */ + set_I(s->I, q, last + i); + set_D(s->D, q, s - 1); + } + set_M(s->M, q, last + i - 1, mat[seq1[i]]); + set_D(s->D, q, s - 1); + if (j + b1 - 1 > len1) { /* bug fixed, 040227 */ + set_end_I(s->I, q, last + i); + } else s->I = MINOR_INF; + s = curr; curr = last; last = s; + } + /* last row for part 1, use set_end_D() instead of set_D() */ + if (j == len2 && b2 != len2 - 1) { + q = dpcell[j]; s = curr; SET_INF(*s); + set_end_I(s->I, q, last); + end = (j + b1 <= len1 + 1)? (j + b1 - 1) : len1; + mat = score_matrix + seq2[j] * N_MATRIX_ROW; + ++s; ++q; + for (i = 1; i != end; ++i, ++s, ++q) { + set_M(s->M, q, last + i - 1, mat[seq1[i]]); /* this will change s->M ! */ + set_I(s->I, q, last + i); + set_end_D(s->D, q, s - 1); + } + set_M(s->M, q, last + i - 1, mat[seq1[i]]); + set_end_D(s->D, q, s - 1); + if (j + b1 - 1 > len1) { /* bug fixed, 040227 */ + set_end_I(s->I, q, last + i); + } else s->I = MINOR_INF; + s = curr; curr = last; last = s; + ++j; + } + + /* core dynamic programming, part 2 */ + for (; j <= len2 - b2 + 1; ++j) { + SET_INF(curr[j - b2]); + mat = score_matrix + seq2[j] * N_MATRIX_ROW; + end = j + b1 - 1; + for (i = j - b2 + 1, q = dpcell[j] + i, s = curr + i; i != end; ++i, ++s, ++q) { + set_M(s->M, q, last + i - 1, mat[seq1[i]]); + set_I(s->I, q, last + i); + set_D(s->D, q, s - 1); + } + set_M(s->M, q, last + i - 1, mat[seq1[i]]); + set_D(s->D, q, s - 1); + s->I = MINOR_INF; + s = curr; curr = last; last = s; + } + + /* core dynamic programming, part 3 */ + for (; j < len2; ++j) { + SET_INF(curr[j - b2]); + mat = score_matrix + seq2[j] * N_MATRIX_ROW; + for (i = j - b2 + 1, q = dpcell[j] + i, s = curr + i; i < len1; ++i, ++s, ++q) { + set_M(s->M, q, last + i - 1, mat[seq1[i]]); + set_I(s->I, q, last + i); + set_D(s->D, q, s - 1); + } + set_M(s->M, q, last + len1 - 1, mat[seq1[i]]); + set_end_I(s->I, q, last + i); + set_D(s->D, q, s - 1); + s = curr; curr = last; last = s; + } + /* last row */ + if (j == len2) { + SET_INF(curr[j - b2]); + mat = score_matrix + seq2[j] * N_MATRIX_ROW; + for (i = j - b2 + 1, q = dpcell[j] + i, s = curr + i; i < len1; ++i, ++s, ++q) { + set_M(s->M, q, last + i - 1, mat[seq1[i]]); + set_I(s->I, q, last + i); + set_end_D(s->D, q, s - 1); + } + set_M(s->M, q, last + len1 - 1, mat[seq1[i]]); + set_end_I(s->I, q, last + i); + set_end_D(s->D, q, s - 1); + s = curr; curr = last; last = s; + } + + *_score = last[len1].M; + if (n_cigar) { /* backtrace */ + path_t *p, *path = (path_t*)malloc(sizeof(path_t) * (len1 + len2 + 2)); + i = len1; j = len2; + q = dpcell[j] + i; + s = last + len1; + max = s->M; type = q->Mt; ctype = FROM_M; + if (s->I > max) { max = s->I; type = q->It; ctype = FROM_I; } + if (s->D > max) { max = s->D; type = q->Dt; ctype = FROM_D; } + + p = path; + p->ctype = ctype; p->i = i; p->j = j; /* bug fixed 040408 */ + ++p; + do { + switch (ctype) { + case FROM_M: --i; --j; break; + case FROM_I: --j; break; + case FROM_D: --i; break; + } + q = dpcell[j] + i; + ctype = type; + switch (type) { + case FROM_M: type = q->Mt; break; + case FROM_I: type = q->It; break; + case FROM_D: type = q->Dt; break; + } + p->ctype = ctype; p->i = i; p->j = j; + ++p; + } while (i || j); + cigar = ka_path2cigar32(path, p - path - 1, n_cigar); + free(path); + } + + /* free memory */ + for (j = b2 + 1; j <= len2; ++j) + dpcell[j] += j - b2; + for (j = 0; j <= len2; ++j) + free(dpcell[j]); + free(dpcell); + free(curr); free(last); + + return cigar; +} diff --git a/samtools/kaln.h b/samtools/kaln.h new file mode 100644 index 0000000..b04d8cc --- /dev/null +++ b/samtools/kaln.h @@ -0,0 +1,55 @@ +/* The MIT License + + Copyright (c) 2003-2006, 2008, 2009 by Heng Li + + Permission is hereby granted, free of charge, to any person obtaining + a copy of this software and associated documentation files (the + "Software"), to deal in the Software without restriction, including + without limitation the rights to use, copy, modify, merge, publish, + distribute, sublicense, and/or sell copies of the Software, and to + permit persons to whom the Software is furnished to do so, subject to + the following conditions: + + The above copyright notice and this permission notice shall be + included in all copies or substantial portions of the Software. + + THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, + EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF + MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND + NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS + BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN + ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN + CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE + SOFTWARE. +*/ + +#ifndef LH3_KALN_H_ +#define LH3_KALN_H_ + +#include + +#define MINOR_INF -1073741823 + +typedef struct { + int gap_open; + int gap_ext; + int gap_end; + + int *matrix; + int row; + int band_width; +} ka_param_t; + +#ifdef __cplusplus +extern "C" { +#endif + + uint32_t *ka_global_core(uint8_t *seq1, int len1, uint8_t *seq2, int len2, const ka_param_t *ap, int *_score, int *n_cigar); + +#ifdef __cplusplus +} +#endif + +extern ka_param_t ka_param_blast; /* = { 5, 2, 2, aln_sm_blast, 5, 50 }; */ + +#endif diff --git a/samtools/klist.h b/samtools/klist.h new file mode 100644 index 0000000..2f17016 --- /dev/null +++ b/samtools/klist.h @@ -0,0 +1,96 @@ +#ifndef _LH3_KLIST_H +#define _LH3_KLIST_H + +#include + +#define KMEMPOOL_INIT(name, kmptype_t, kmpfree_f) \ + typedef struct { \ + size_t cnt, n, max; \ + kmptype_t **buf; \ + } kmp_##name##_t; \ + static inline kmp_##name##_t *kmp_init_##name() { \ + return calloc(1, sizeof(kmp_##name##_t)); \ + } \ + static inline void kmp_destroy_##name(kmp_##name##_t *mp) { \ + size_t k; \ + for (k = 0; k < mp->n; ++k) { \ + kmpfree_f(mp->buf[k]); free(mp->buf[k]); \ + } \ + free(mp->buf); free(mp); \ + } \ + static inline kmptype_t *kmp_alloc_##name(kmp_##name##_t *mp) { \ + ++mp->cnt; \ + if (mp->n == 0) return calloc(1, sizeof(kmptype_t)); \ + return mp->buf[--mp->n]; \ + } \ + static inline void kmp_free_##name(kmp_##name##_t *mp, kmptype_t *p) { \ + --mp->cnt; \ + if (mp->n == mp->max) { \ + mp->max = mp->max? mp->max<<1 : 16; \ + mp->buf = realloc(mp->buf, sizeof(void*) * mp->max); \ + } \ + mp->buf[mp->n++] = p; \ + } + +#define kmempool_t(name) kmp_##name##_t +#define kmp_init(name) kmp_init_##name() +#define kmp_destroy(name, mp) kmp_destroy_##name(mp) +#define kmp_alloc(name, mp) kmp_alloc_##name(mp) +#define kmp_free(name, mp, p) kmp_free_##name(mp, p) + +#define KLIST_INIT(name, kltype_t, kmpfree_t) \ + struct __kl1_##name { \ + kltype_t data; \ + struct __kl1_##name *next; \ + }; \ + typedef struct __kl1_##name kl1_##name; \ + KMEMPOOL_INIT(name, kl1_##name, kmpfree_t) \ + typedef struct { \ + kl1_##name *head, *tail; \ + kmp_##name##_t *mp; \ + size_t size; \ + } kl_##name##_t; \ + static inline kl_##name##_t *kl_init_##name() { \ + kl_##name##_t *kl = calloc(1, sizeof(kl_##name##_t)); \ + kl->mp = kmp_init(name); \ + kl->head = kl->tail = kmp_alloc(name, kl->mp); \ + kl->head->next = 0; \ + return kl; \ + } \ + static inline void kl_destroy_##name(kl_##name##_t *kl) { \ + kl1_##name *p; \ + for (p = kl->head; p != kl->tail; p = p->next) \ + kmp_free(name, kl->mp, p); \ + kmp_free(name, kl->mp, p); \ + kmp_destroy(name, kl->mp); \ + free(kl); \ + } \ + static inline kltype_t *kl_pushp_##name(kl_##name##_t *kl) { \ + kl1_##name *q, *p = kmp_alloc(name, kl->mp); \ + q = kl->tail; p->next = 0; kl->tail->next = p; kl->tail = p; \ + ++kl->size; \ + return &q->data; \ + } \ + static inline int kl_shift_##name(kl_##name##_t *kl, kltype_t *d) { \ + kl1_##name *p; \ + if (kl->head->next == 0) return -1; \ + --kl->size; \ + p = kl->head; kl->head = kl->head->next; \ + if (d) *d = p->data; \ + kmp_free(name, kl->mp, p); \ + return 0; \ + } + +#define kliter_t(name) kl1_##name +#define klist_t(name) kl_##name##_t +#define kl_val(iter) ((iter)->data) +#define kl_next(iter) ((iter)->next) +#define kl_begin(kl) ((kl)->head) +#define kl_end(kl) ((kl)->tail) + +#define kl_init(name) kl_init_##name() +#define kl_destroy(name, kl) kl_destroy_##name(kl) +#define kl_pushp(name, kl) kl_pushp_##name(kl) +#define kl_shift(name, kl, d) kl_shift_##name(kl, d) + +#endif diff --git a/samtools/knetfile.c b/samtools/knetfile.c index e110aa7..994babb 100644 --- a/samtools/knetfile.c +++ b/samtools/knetfile.c @@ -34,6 +34,7 @@ #include #include #include +#include #include #include @@ -70,7 +71,14 @@ static int socket_wait(int fd, int is_read) if (is_read) fdr = &fds; else fdw = &fds; ret = select(fd+1, fdr, fdw, 0, &tv); +#ifndef _WIN32 if (ret == -1) perror("select"); +#else + if (ret == 0) + fprintf(stderr, "select time-out\n"); + else if (ret == SOCKET_ERROR) + fprintf(stderr, "select: %d\n", WSAGetLastError()); +#endif return ret; } @@ -103,16 +111,28 @@ static int socket_connect(const char *host, const char *port) } #else /* MinGW's printf has problem with "%lld" */ -char *uint64tostr(char *buf, uint64_t x) +char *int64tostr(char *buf, int64_t x) { - int i, cnt; - for (i = 0; x; x /= 10) buf[i++] = '0' + x%10; + int cnt; + int i = 0; + do { + buf[i++] = '0' + x % 10; + x /= 10; + } while (x); buf[i] = 0; for (cnt = i, i = 0; i < cnt/2; ++i) { int c = buf[i]; buf[i] = buf[cnt-i-1]; buf[cnt-i-1] = c; } return buf; } + +int64_t strtoint64(const char *buf) +{ + int64_t x; + for (x = 0; *buf != '\0'; ++buf) + x = x * 10 + ((int64_t) *buf - 48); + return x; +} /* In windows, the first thing is to establish the TCP connection. */ int knet_win32_init() { @@ -129,7 +149,11 @@ void knet_win32_destroy() * non-Windows OS, I do not use this one. */ static SOCKET socket_connect(const char *host, const char *port) { -#define __err_connect(func) do { perror(func); return -1; } while (0) +#define __err_connect(func) \ + do { \ + fprintf(stderr, "%s: %d\n", func, WSAGetLastError()); \ + return -1; \ + } while (0) int on = 1; SOCKET fd; @@ -182,7 +206,11 @@ static off_t my_netread(int fd, void *buf, off_t len) static int kftp_get_response(knetFile *ftp) { +#ifndef _WIN32 unsigned char c; +#else + char c; +#endif int n = 0; char *p; if (socket_wait(ftp->ctrl_fd, 1) <= 0) return 0; @@ -259,10 +287,11 @@ int kftp_reconnect(knetFile *ftp) ftp->ctrl_fd = -1; } netclose(ftp->fd); + ftp->fd = -1; return kftp_connect(ftp); } -// initialize ->type, ->host and ->retr +// initialize ->type, ->host, ->retr and ->size knetFile *kftp_parse_url(const char *fn, const char *mode) { knetFile *fp; @@ -283,25 +312,42 @@ knetFile *kftp_parse_url(const char *fn, const char *mode) strncpy(fp->host, fn + 6, l); fp->retr = calloc(strlen(p) + 8, 1); sprintf(fp->retr, "RETR %s\r\n", p); - fp->seek_offset = -1; + fp->size_cmd = calloc(strlen(p) + 8, 1); + sprintf(fp->size_cmd, "SIZE %s\r\n", p); + fp->seek_offset = 0; return fp; } // place ->fd at offset off int kftp_connect_file(knetFile *fp) { int ret; + long long file_size; if (fp->fd != -1) { netclose(fp->fd); if (fp->no_reconnect) kftp_get_response(fp); } kftp_pasv_prep(fp); - if (fp->offset) { + kftp_send_cmd(fp, fp->size_cmd, 1); +#ifndef _WIN32 + if ( sscanf(fp->response,"%*d %lld", &file_size) != 1 ) + { + fprintf(stderr,"[kftp_connect_file] %s\n", fp->response); + return -1; + } +#else + const char *p = fp->response; + while (*p != ' ') ++p; + while (*p < '0' || *p > '9') ++p; + file_size = strtoint64(p); +#endif + fp->file_size = file_size; + if (fp->offset>=0) { char tmp[32]; #ifndef _WIN32 sprintf(tmp, "REST %lld\r\n", (long long)fp->offset); #else strcpy(tmp, "REST "); - uint64tostr(tmp + 5, fp->offset); + int64tostr(tmp + 5, fp->offset); strcat(tmp, "\r\n"); #endif kftp_send_cmd(fp, tmp, 1); @@ -319,6 +365,7 @@ int kftp_connect_file(knetFile *fp) return 0; } + /************************** * HTTP specific routines * **************************/ @@ -354,7 +401,7 @@ knetFile *khttp_parse_url(const char *fn, const char *mode) } fp->type = KNF_TYPE_HTTP; fp->ctrl_fd = fp->fd = -1; - fp->seek_offset = -1; + fp->seek_offset = 0; return fp; } @@ -366,8 +413,7 @@ int khttp_connect_file(knetFile *fp) fp->fd = socket_connect(fp->host, fp->port); buf = calloc(0x10000, 1); // FIXME: I am lazy... But in principle, 64KB should be large enough. l += sprintf(buf + l, "GET %s HTTP/1.0\r\nHost: %s\r\n", fp->path, fp->http_host); - if (fp->offset) - l += sprintf(buf + l, "Range: bytes=%lld-\r\n", (long long)fp->offset); + l += sprintf(buf + l, "Range: bytes=%lld-\r\n", (long long)fp->offset); l += sprintf(buf + l, "\r\n"); netwrite(fp->fd, buf, l); l = 0; @@ -383,7 +429,7 @@ int khttp_connect_file(knetFile *fp) return -1; } ret = strtol(buf + 8, &p, 0); // HTTP return code - if (ret == 200 && fp->offset) { // 200 (complete result); then skip beginning of the file + if (ret == 200 && fp->offset>0) { // 200 (complete result); then skip beginning of the file off_t rest = fp->offset; while (rest) { off_t l = rest < 0x10000? rest : 0x10000; @@ -482,7 +528,7 @@ off_t knet_read(knetFile *fp, void *buf, off_t len) return l; } -int knet_seek(knetFile *fp, off_t off, int whence) +off_t knet_seek(knetFile *fp, int64_t off, int whence) { if (whence == SEEK_SET && off == fp->offset) return 0; if (fp->type == KNF_TYPE_LOCAL) { @@ -490,20 +536,40 @@ int knet_seek(knetFile *fp, off_t off, int whence) * while fseek() returns zero on success. */ off_t offset = lseek(fp->fd, off, whence); if (offset == -1) { - perror("lseek"); + // Be silent, it is OK for knet_seek to fail when the file is streamed + // fprintf(stderr,"[knet_seek] %s\n", strerror(errno)); return -1; } fp->offset = offset; return 0; - } else if (fp->type == KNF_TYPE_FTP || fp->type == KNF_TYPE_HTTP) { - if (whence != SEEK_SET) { // FIXME: we can surely allow SEEK_CUR and SEEK_END in future - fprintf(stderr, "[knet_seek] only SEEK_SET is supported for FTP/HTTP. Offset is unchanged.\n"); + } + else if (fp->type == KNF_TYPE_FTP) + { + if (whence==SEEK_CUR) + fp->offset += off; + else if (whence==SEEK_SET) + fp->offset = off; + else if ( whence==SEEK_END) + fp->offset = fp->file_size+off; + fp->is_ready = 0; + return 0; + } + else if (fp->type == KNF_TYPE_HTTP) + { + if (whence == SEEK_END) { // FIXME: can we allow SEEK_END in future? + fprintf(stderr, "[knet_seek] SEEK_END is not supported for HTTP. Offset is unchanged.\n"); + errno = ESPIPE; return -1; } - fp->offset = off; + if (whence==SEEK_CUR) + fp->offset += off; + else if (whence==SEEK_SET) + fp->offset = off; fp->is_ready = 0; - return 0; + return fp->offset; } + errno = EINVAL; + fprintf(stderr,"[knet_seek] %s\n", strerror(errno)); return -1; } diff --git a/samtools/knetfile.h b/samtools/knetfile.h index 9021b93..0a0e66f 100644 --- a/samtools/knetfile.h +++ b/samtools/knetfile.h @@ -9,7 +9,7 @@ #define netwrite(fd, ptr, len) write(fd, ptr, len) #define netclose(fd) close(fd) #else -#include +#include #define netread(fd, ptr, len) recv(fd, ptr, len, 0) #define netwrite(fd, ptr, len) send(fd, ptr, len, 0) #define netclose(fd) closesocket(fd) @@ -28,8 +28,9 @@ typedef struct knetFile_s { // the following are for FTP only int ctrl_fd, pasv_ip[4], pasv_port, max_response, no_reconnect, is_ready; - char *response, *retr; + char *response, *retr, *size_cmd; int64_t seek_offset; // for lazy seek + int64_t file_size; // the following are for HTTP only char *path, *http_host; @@ -64,7 +65,7 @@ extern "C" { This routine only sets ->offset and ->is_ready=0. It does not communicate with the FTP server. */ - int knet_seek(knetFile *fp, off_t off, int whence); + off_t knet_seek(knetFile *fp, int64_t off, int whence); int knet_close(knetFile *fp); #ifdef __cplusplus diff --git a/samtools/razf.c b/samtools/razf.c index a5e8f51..e7499f9 100644 --- a/samtools/razf.c +++ b/samtools/razf.c @@ -38,6 +38,7 @@ #include #include "razf.h" + #if ZLIB_VERNUM < 0x1221 struct _gz_header_s { int text; @@ -107,20 +108,36 @@ static void save_zindex(RAZF *rz, int fd){ } #endif +#ifdef _USE_KNETFILE +static void load_zindex(RAZF *rz, knetFile *fp){ +#else static void load_zindex(RAZF *rz, int fd){ +#endif int32_t i, v32; int is_be; if(!rz->load_index) return; if(rz->index == NULL) rz->index = malloc(sizeof(ZBlockIndex)); is_be = is_big_endian(); +#ifdef _USE_KNETFILE + knet_read(fp, &rz->index->size, sizeof(int)); +#else read(fd, &rz->index->size, sizeof(int)); +#endif if(!is_be) rz->index->size = byte_swap_4((uint32_t)rz->index->size); rz->index->cap = rz->index->size; v32 = rz->index->size / RZ_BIN_SIZE + 1; rz->index->bin_offsets = malloc(sizeof(int64_t) * v32); +#ifdef _USE_KNETFILE + knet_read(fp, rz->index->bin_offsets, sizeof(int64_t) * v32); +#else read(fd, rz->index->bin_offsets, sizeof(int64_t) * v32); +#endif rz->index->cell_offsets = malloc(sizeof(int) * rz->index->size); +#ifdef _USE_KNETFILE + knet_read(fp, rz->index->cell_offsets, sizeof(int) * rz->index->size); +#else read(fd, rz->index->cell_offsets, sizeof(int) * rz->index->size); +#endif if(!is_be){ for(i=0;iindex->bin_offsets[i] = byte_swap_8((uint64_t)rz->index->bin_offsets[i]); for(i=0;iindex->size;i++) rz->index->cell_offsets[i] = byte_swap_4((uint32_t)rz->index->cell_offsets[i]); @@ -141,7 +158,11 @@ static RAZF* razf_open_w(int fd){ #endif rz = calloc(1, sizeof(RAZF)); rz->mode = 'w'; +#ifdef _USE_KNETFILE + rz->x.fpw = fd; +#else rz->filedes = fd; +#endif rz->stream = calloc(sizeof(z_stream), 1); rz->inbuf = malloc(RZ_BUFFER_SIZE); rz->outbuf = malloc(RZ_BUFFER_SIZE); @@ -176,7 +197,11 @@ static void _razf_write(RAZF* rz, const void *data, int size){ deflate(rz->stream, Z_NO_FLUSH); rz->out += tout - rz->stream->avail_out; if(rz->stream->avail_out) break; +#ifdef _USE_KNETFILE + write(rz->x.fpw, rz->outbuf, RZ_BUFFER_SIZE - rz->stream->avail_out); +#else write(rz->filedes, rz->outbuf, RZ_BUFFER_SIZE - rz->stream->avail_out); +#endif rz->stream->avail_out = RZ_BUFFER_SIZE; rz->stream->next_out = rz->outbuf; if(rz->stream->avail_in == 0) break; @@ -192,7 +217,11 @@ static void razf_flush(RAZF *rz){ rz->buf_off = rz->buf_len = 0; } if(rz->stream->avail_out){ +#ifdef _USE_KNETFILE + write(rz->x.fpw, rz->outbuf, RZ_BUFFER_SIZE - rz->stream->avail_out); +#else write(rz->filedes, rz->outbuf, RZ_BUFFER_SIZE - rz->stream->avail_out); +#endif rz->stream->avail_out = RZ_BUFFER_SIZE; rz->stream->next_out = rz->outbuf; } @@ -201,7 +230,11 @@ static void razf_flush(RAZF *rz){ deflate(rz->stream, Z_FULL_FLUSH); rz->out += tout - rz->stream->avail_out; if(rz->stream->avail_out == 0){ +#ifdef _USE_KNETFILE + write(rz->x.fpw, rz->outbuf, RZ_BUFFER_SIZE - rz->stream->avail_out); +#else write(rz->filedes, rz->outbuf, RZ_BUFFER_SIZE - rz->stream->avail_out); +#endif rz->stream->avail_out = RZ_BUFFER_SIZE; rz->stream->next_out = rz->outbuf; } else break; @@ -221,7 +254,11 @@ static void razf_end_flush(RAZF *rz){ deflate(rz->stream, Z_FINISH); rz->out += tout - rz->stream->avail_out; if(rz->stream->avail_out < RZ_BUFFER_SIZE){ +#ifdef _USE_KNETFILE + write(rz->x.fpw, rz->outbuf, RZ_BUFFER_SIZE - rz->stream->avail_out); +#else write(rz->filedes, rz->outbuf, RZ_BUFFER_SIZE - rz->stream->avail_out); +#endif rz->stream->avail_out = RZ_BUFFER_SIZE; rz->stream->next_out = rz->outbuf; } else break; @@ -308,23 +345,35 @@ static int _read_gz_header(unsigned char *data, int size, int *extra_off, int *e return n; } +#ifdef _USE_KNETFILE +static RAZF* razf_open_r(knetFile *fp, int _load_index){ +#else static RAZF* razf_open_r(int fd, int _load_index){ +#endif RAZF *rz; int ext_off, ext_len; int n, is_be, ret; int64_t end; unsigned char c[] = "RAZF"; + rz = calloc(1, sizeof(RAZF)); + rz->mode = 'r'; +#ifdef _USE_KNETFILE + rz->x.fpr = fp; +#else #ifdef _WIN32 setmode(fd, O_BINARY); #endif - rz = calloc(1, sizeof(RAZF)); - rz->mode = 'r'; rz->filedes = fd; +#endif rz->stream = calloc(sizeof(z_stream), 1); rz->inbuf = malloc(RZ_BUFFER_SIZE); rz->outbuf = malloc(RZ_BUFFER_SIZE); rz->end = rz->src_end = 0x7FFFFFFFFFFFFFFFLL; +#ifdef _USE_KNETFILE + n = knet_read(rz->x.fpr, rz->inbuf, RZ_BUFFER_SIZE); +#else n = read(rz->filedes, rz->inbuf, RZ_BUFFER_SIZE); +#endif ret = _read_gz_header(rz->inbuf, n, &ext_off, &ext_len); if(ret == 0){ PLAIN_FILE: @@ -355,7 +404,11 @@ static RAZF* razf_open_r(int fd, int _load_index){ } rz->load_index = _load_index; rz->file_type = FILE_TYPE_RZ; +#ifdef _USE_KNETFILE + if(knet_seek(fp, -16, SEEK_END) == -1){ +#else if(lseek(fd, -16, SEEK_END) == -1){ +#endif UNSEEKABLE: rz->seekable = 0; rz->index = NULL; @@ -363,10 +416,19 @@ static RAZF* razf_open_r(int fd, int _load_index){ } else { is_be = is_big_endian(); rz->seekable = 1; +#ifdef _USE_KNETFILE + knet_read(fp, &end, sizeof(int64_t)); +#else read(fd, &end, sizeof(int64_t)); +#endif if(!is_be) rz->src_end = (int64_t)byte_swap_8((uint64_t)end); else rz->src_end = end; + +#ifdef _USE_KNETFILE + knet_read(fp, &end, sizeof(int64_t)); +#else read(fd, &end, sizeof(int64_t)); +#endif if(!is_be) rz->end = (int64_t)byte_swap_8((uint64_t)end); else rz->end = end; if(n > rz->end){ @@ -374,19 +436,47 @@ static RAZF* razf_open_r(int fd, int _load_index){ n = rz->end; } if(rz->end > rz->src_end){ +#ifdef _USE_KNETFILE + knet_seek(fp, rz->in, SEEK_SET); +#else lseek(fd, rz->in, SEEK_SET); +#endif goto UNSEEKABLE; } +#ifdef _USE_KNETFILE + knet_seek(fp, rz->end, SEEK_SET); + if(knet_tell(fp) != rz->end){ + knet_seek(fp, rz->in, SEEK_SET); +#else if(lseek(fd, rz->end, SEEK_SET) != rz->end){ lseek(fd, rz->in, SEEK_SET); +#endif goto UNSEEKABLE; } +#ifdef _USE_KNETFILE + load_zindex(rz, fp); + knet_seek(fp, n, SEEK_SET); +#else load_zindex(rz, fd); lseek(fd, n, SEEK_SET); +#endif } return rz; } +#ifdef _USE_KNETFILE +RAZF* razf_dopen(int fd, const char *mode){ + if (strstr(mode, "r")) fprintf(stderr,"[razf_dopen] implement me\n"); + else if(strstr(mode, "w")) return razf_open_w(fd); + return NULL; +} + +RAZF* razf_dopen2(int fd, const char *mode) +{ + fprintf(stderr,"[razf_dopen2] implement me\n"); + return NULL; +} +#else RAZF* razf_dopen(int fd, const char *mode){ if(strstr(mode, "r")) return razf_open_r(fd, 1); else if(strstr(mode, "w")) return razf_open_w(fd); @@ -399,23 +489,34 @@ RAZF* razf_dopen2(int fd, const char *mode) else if(strstr(mode, "w")) return razf_open_w(fd); else return NULL; } +#endif static inline RAZF* _razf_open(const char *filename, const char *mode, int _load_index){ int fd; RAZF *rz; if(strstr(mode, "r")){ +#ifdef _USE_KNETFILE + knetFile *fd = knet_open(filename, "r"); + if (fd == 0) { + fprintf(stderr, "[_razf_open] fail to open %s\n", filename); + return NULL; + } +#else #ifdef _WIN32 fd = open(filename, O_RDONLY | O_BINARY); #else fd = open(filename, O_RDONLY); #endif +#endif + if(fd < 0) return NULL; rz = razf_open_r(fd, _load_index); } else if(strstr(mode, "w")){ #ifdef _WIN32 - fd = open(filename, O_WRONLY | O_CREAT | O_TRUNC | O_BINARY, 0644); + fd = open(filename, O_WRONLY | O_CREAT | O_TRUNC | O_BINARY, 0666); #else - fd = open(filename, O_WRONLY | O_CREAT | O_TRUNC, 0644); + fd = open(filename, O_WRONLY | O_CREAT | O_TRUNC, 0666); #endif + if(fd < 0) return NULL; rz = razf_open_w(fd); } else return NULL; return rz; @@ -435,9 +536,17 @@ int razf_get_data_size(RAZF *rz, int64_t *u_size, int64_t *c_size){ switch(rz->file_type){ case FILE_TYPE_PLAIN: if(rz->end == 0x7fffffffffffffffLL){ +#ifdef _USE_KNETFILE + if(knet_seek(rz->x.fpr, 0, SEEK_CUR) == -1) return 0; + n = knet_tell(rz->x.fpr); + knet_seek(rz->x.fpr, 0, SEEK_END); + rz->end = knet_tell(rz->x.fpr); + knet_seek(rz->x.fpr, n, SEEK_SET); +#else if((n = lseek(rz->filedes, 0, SEEK_CUR)) == -1) return 0; rz->end = lseek(rz->filedes, 0, SEEK_END); lseek(rz->filedes, n, SEEK_SET); +#endif } *u_size = *c_size = rz->end; return 1; @@ -457,7 +566,11 @@ static int _razf_read(RAZF* rz, void *data, int size){ int ret, tin; if(rz->z_eof || rz->z_err) return 0; if (rz->file_type == FILE_TYPE_PLAIN) { +#ifdef _USE_KNETFILE + ret = knet_read(rz->x.fpr, data, size); +#else ret = read(rz->filedes, data, size); +#endif if (ret == 0) rz->z_eof = 1; return ret; } @@ -467,9 +580,17 @@ static int _razf_read(RAZF* rz, void *data, int size){ if(rz->stream->avail_in == 0){ if(rz->in >= rz->end){ rz->z_eof = 1; break; } if(rz->end - rz->in < RZ_BUFFER_SIZE){ +#ifdef _USE_KNETFILE + rz->stream->avail_in = knet_read(rz->x.fpr, rz->inbuf, rz->end -rz->in); +#else rz->stream->avail_in = read(rz->filedes, rz->inbuf, rz->end -rz->in); +#endif } else { +#ifdef _USE_KNETFILE + rz->stream->avail_in = knet_read(rz->x.fpr, rz->inbuf, RZ_BUFFER_SIZE); +#else rz->stream->avail_in = read(rz->filedes, rz->inbuf, RZ_BUFFER_SIZE); +#endif } if(rz->stream->avail_in == 0){ rz->z_eof = 1; @@ -481,7 +602,7 @@ static int _razf_read(RAZF* rz, void *data, int size){ ret = inflate(rz->stream, Z_BLOCK); rz->in += tin - rz->stream->avail_in; if(ret == Z_NEED_DICT || ret == Z_MEM_ERROR || ret == Z_DATA_ERROR){ - fprintf(stderr, "[_razf_read] inflate error: %d (at %s:%d)\n", ret, __FILE__, __LINE__); + fprintf(stderr, "[_razf_read] inflate error: %d %s (at %s:%d)\n", ret, rz->stream->msg ? rz->stream->msg : "", __FILE__, __LINE__); rz->z_err = 1; break; } @@ -566,14 +687,18 @@ int razf_skip(RAZF* rz, int size){ } if(rz->buf_flush) continue; rz->buf_len = _razf_read(rz, rz->outbuf, RZ_BUFFER_SIZE); - if(rz->z_eof) break; + if(rz->z_eof || rz->z_err) break; } rz->out += ori_size - size; return ori_size - size; } static void _razf_reset_read(RAZF *rz, int64_t in, int64_t out){ +#ifdef _USE_KNETFILE + knet_seek(rz->x.fpr, in, SEEK_SET); +#else lseek(rz->filedes, in, SEEK_SET); +#endif rz->in = in; rz->out = out; rz->block_pos = in; @@ -592,7 +717,12 @@ int64_t razf_jump(RAZF *rz, int64_t block_start, int block_offset){ if(rz->file_type == FILE_TYPE_PLAIN){ rz->buf_off = rz->buf_len = 0; pos = block_start + block_offset; +#ifdef _USE_KNETFILE + knet_seek(rz->x.fpr, pos, SEEK_SET); + pos = knet_tell(rz->x.fpr); +#else pos = lseek(rz->filedes, pos, SEEK_SET); +#endif rz->out = rz->in = pos; return pos; } @@ -614,7 +744,12 @@ int64_t razf_seek(RAZF* rz, int64_t pos, int where){ if (where == SEEK_CUR) pos += rz->out; else if (where == SEEK_END) pos += rz->src_end; if(rz->file_type == FILE_TYPE_PLAIN){ +#ifdef _USE_KNETFILE + knet_seek(rz->x.fpr, pos, SEEK_SET); + seek_pos = knet_tell(rz->x.fpr); +#else seek_pos = lseek(rz->filedes, pos, SEEK_SET); +#endif rz->buf_off = rz->buf_len = 0; rz->out = rz->in = seek_pos; return seek_pos; @@ -663,6 +798,18 @@ void razf_close(RAZF *rz){ #ifndef _RZ_READONLY razf_end_flush(rz); deflateEnd(rz->stream); +#ifdef _USE_KNETFILE + save_zindex(rz, rz->x.fpw); + if(is_big_endian()){ + write(rz->x.fpw, &rz->in, sizeof(int64_t)); + write(rz->x.fpw, &rz->out, sizeof(int64_t)); + } else { + uint64_t v64 = byte_swap_8((uint64_t)rz->in); + write(rz->x.fpw, &v64, sizeof(int64_t)); + v64 = byte_swap_8((uint64_t)rz->out); + write(rz->x.fpw, &v64, sizeof(int64_t)); + } +#else save_zindex(rz, rz->filedes); if(is_big_endian()){ write(rz->filedes, &rz->in, sizeof(int64_t)); @@ -673,6 +820,7 @@ void razf_close(RAZF *rz){ v64 = byte_swap_8((uint64_t)rz->out); write(rz->filedes, &v64, sizeof(int64_t)); } +#endif #endif } else if(rz->mode == 'r'){ if(rz->stream) inflateEnd(rz->stream); @@ -691,7 +839,14 @@ void razf_close(RAZF *rz){ free(rz->index); } free(rz->stream); +#ifdef _USE_KNETFILE + if (rz->mode == 'r') + knet_close(rz->x.fpr); + if (rz->mode == 'w') + close(rz->x.fpw); +#else close(rz->filedes); +#endif free(rz); } diff --git a/samtools/razf.h b/samtools/razf.h index f7e5097..60a0c96 100644 --- a/samtools/razf.h +++ b/samtools/razf.h @@ -37,6 +37,10 @@ #include #include "zlib.h" +#ifdef _USE_KNETFILE +#include "knetfile.h" +#endif + #if ZLIB_VERNUM < 0x1221 #define _RZ_READONLY struct _gz_header_s; @@ -76,7 +80,14 @@ typedef struct RandomAccessZFile { char mode; /* 'w' : write mode; 'r' : read mode */ int file_type; /* plain file or rz file, razf_read support plain file as input too, in this case, razf_read work as buffered fread */ +#ifdef _USE_KNETFILE + union { + knetFile *fpr; + int fpw; + } x; +#else int filedes; /* the file descriptor */ +#endif z_stream *stream; ZBlockIndex *index; int64_t in, out, end, src_end; diff --git a/samtools/sam.c b/samtools/sam.c index a74c557..ad4325b 100644 --- a/samtools/sam.c +++ b/samtools/sam.c @@ -12,7 +12,7 @@ bam_header_t *bam_header_dup(const bam_header_t *h0) int i; h = bam_header_init(); *h = *h0; - h->hash = 0; + h->hash = h->dict = h->rg2lib = 0; h->text = (char*)calloc(h->l_text + 1, 1); memcpy(h->text, h0->text, h->l_text); h->target_len = (uint32_t*)calloc(h->n_targets, 4); @@ -21,7 +21,6 @@ bam_header_t *bam_header_dup(const bam_header_t *h0) h->target_len[i] = h0->target_len[i]; h->target_name[i] = strdup(h0->target_name[i]); } - if (h0->rg2lib) h->rg2lib = bam_strmap_dup(h0->rg2lib); return h; } static void append_header_text(bam_header_t *header, char* text, int len) @@ -63,7 +62,6 @@ samfile_t *samopen(const char *fn, const char *mode, const void *aux) fprintf(stderr, "[samopen] no @SQ lines in the header.\n"); } else fprintf(stderr, "[samopen] SAM header is present: %d sequences.\n", fp->header->n_targets); } - sam_header_parse_rg(fp->header); } else if (mode[0] == 'w') { // write fp->header = bam_header_dup((const bam_header_t*)aux); if (mode[1] == 'b') { // binary diff --git a/samtools/sam_header.c b/samtools/sam_header.c new file mode 100644 index 0000000..a119c02 --- /dev/null +++ b/samtools/sam_header.c @@ -0,0 +1,701 @@ +#include "sam_header.h" +#include +#include +#include +#include +#include + +#include "khash.h" +KHASH_MAP_INIT_STR(str, const char *) + +struct _HeaderList +{ + struct _HeaderList *next; + void *data; +}; +typedef struct _HeaderList list_t; +typedef list_t HeaderDict; + +typedef struct +{ + char key[2]; + char *value; +} +HeaderTag; + +typedef struct +{ + char type[2]; + list_t *tags; +} +HeaderLine; + +const char *o_hd_tags[] = {"SO","GO",NULL}; +const char *r_hd_tags[] = {"VN",NULL}; + +const char *o_sq_tags[] = {"AS","M5","UR","SP",NULL}; +const char *r_sq_tags[] = {"SN","LN",NULL}; +const char *u_sq_tags[] = {"SN",NULL}; + +const char *o_rg_tags[] = {"LB","DS","PU","PI","CN","DT","PL",NULL}; +const char *r_rg_tags[] = {"ID",NULL}; +const char *u_rg_tags[] = {"ID",NULL}; + +const char *o_pg_tags[] = {"VN","CL",NULL}; +const char *r_pg_tags[] = {"ID",NULL}; + +const char *types[] = {"HD","SQ","RG","PG","CO",NULL}; +const char **optional_tags[] = {o_hd_tags,o_sq_tags,o_rg_tags,o_pg_tags,NULL,NULL}; +const char **required_tags[] = {r_hd_tags,r_sq_tags,r_rg_tags,r_pg_tags,NULL,NULL}; +const char **unique_tags[] = {NULL, u_sq_tags,u_rg_tags,NULL,NULL,NULL}; + + +static void debug(const char *format, ...) +{ + va_list ap; + va_start(ap, format); + vfprintf(stderr, format, ap); + va_end(ap); +} + +static list_t *list_append(list_t *root, void *data) +{ + list_t *l = root; + while (l && l->next) + l = l->next; + if ( l ) + { + l->next = malloc(sizeof(list_t)); + l = l->next; + } + else + { + l = malloc(sizeof(list_t)); + root = l; + } + l->data = data; + l->next = NULL; + return root; +} + +static void list_free(list_t *root) +{ + list_t *l = root; + while (root) + { + l = root; + root = root->next; + free(l); + } +} + + + +// Look for a tag "XY" in a predefined const char *[] array. +static int tag_exists(const char *tag, const char **tags) +{ + int itag=0; + if ( !tags ) return -1; + while ( tags[itag] ) + { + if ( tags[itag][0]==tag[0] && tags[itag][1]==tag[1] ) return itag; + itag++; + } + return -1; +} + + + +// Mimics the behaviour of getline, except it returns pointer to the next chunk of the text +// or NULL if everything has been read. The lineptr should be freed by the caller. The +// newline character is stripped. +static const char *nextline(char **lineptr, size_t *n, const char *text) +{ + int len; + const char *to = text; + + if ( !*to ) return NULL; + + while ( *to && *to!='\n' && *to!='\r' ) to++; + len = to - text + 1; + + if ( *to ) + { + // Advance the pointer for the next call + if ( *to=='\n' ) to++; + else if ( *to=='\r' && *(to+1)=='\n' ) to+=2; + } + if ( !len ) + return to; + + if ( !*lineptr ) + { + *lineptr = malloc(len); + *n = len; + } + else if ( *nkey[0] = name[0]; + tag->key[1] = name[1]; + tag->value = malloc(len+1); + memcpy(tag->value,value_from,len+1); + tag->value[len] = 0; + return tag; +} + +static HeaderTag *header_line_has_tag(HeaderLine *hline, const char *key) +{ + list_t *tags = hline->tags; + while (tags) + { + HeaderTag *tag = tags->data; + if ( tag->key[0]==key[0] && tag->key[1]==key[1] ) return tag; + tags = tags->next; + } + return NULL; +} + + +// Return codes: +// 0 .. different types or unique tags differ or conflicting tags, cannot be merged +// 1 .. all tags identical -> no need to merge, drop one +// 2 .. the unique tags match and there are some conflicting tags (same tag, different value) -> error, cannot be merged nor duplicated +// 3 .. there are some missing complementary tags and no unique conflict -> can be merged into a single line +static int sam_header_compare_lines(HeaderLine *hline1, HeaderLine *hline2) +{ + HeaderTag *t1, *t2; + + if ( hline1->type[0]!=hline2->type[0] || hline1->type[1]!=hline2->type[1] ) + return 0; + + int itype = tag_exists(hline1->type,types); + if ( itype==-1 ) { + debug("[sam_header_compare_lines] Unknown type [%c%c]\n", hline1->type[0],hline1->type[1]); + return -1; // FIXME (lh3): error; I do not know how this will be handled in Petr's code + } + + if ( unique_tags[itype] ) + { + t1 = header_line_has_tag(hline1,unique_tags[itype][0]); + t2 = header_line_has_tag(hline2,unique_tags[itype][0]); + if ( !t1 || !t2 ) // this should never happen, the unique tags are required + return 2; + + if ( strcmp(t1->value,t2->value) ) + return 0; // the unique tags differ, cannot be merged + } + if ( !required_tags[itype] && !optional_tags[itype] ) + { + t1 = hline1->tags->data; + t2 = hline2->tags->data; + if ( !strcmp(t1->value,t2->value) ) return 1; // identical comments + return 0; + } + + int missing=0, itag=0; + while ( required_tags[itype] && required_tags[itype][itag] ) + { + t1 = header_line_has_tag(hline1,required_tags[itype][itag]); + t2 = header_line_has_tag(hline2,required_tags[itype][itag]); + if ( !t1 && !t2 ) + return 2; // this should never happen + else if ( !t1 || !t2 ) + missing = 1; // there is some tag missing in one of the hlines + else if ( strcmp(t1->value,t2->value) ) + { + if ( unique_tags[itype] ) + return 2; // the lines have a matching unique tag but have a conflicting tag + + return 0; // the lines contain conflicting tags, cannot be merged + } + itag++; + } + itag = 0; + while ( optional_tags[itype] && optional_tags[itype][itag] ) + { + t1 = header_line_has_tag(hline1,optional_tags[itype][itag]); + t2 = header_line_has_tag(hline2,optional_tags[itype][itag]); + if ( !t1 && !t2 ) + { + itag++; + continue; + } + if ( !t1 || !t2 ) + missing = 1; // there is some tag missing in one of the hlines + else if ( strcmp(t1->value,t2->value) ) + { + if ( unique_tags[itype] ) + return 2; // the lines have a matching unique tag but have a conflicting tag + + return 0; // the lines contain conflicting tags, cannot be merged + } + itag++; + } + if ( missing ) return 3; // there are some missing complementary tags with no conflicts, can be merged + return 1; +} + + +static HeaderLine *sam_header_line_clone(const HeaderLine *hline) +{ + list_t *tags; + HeaderLine *out = malloc(sizeof(HeaderLine)); + out->type[0] = hline->type[0]; + out->type[1] = hline->type[1]; + out->tags = NULL; + + tags = hline->tags; + while (tags) + { + HeaderTag *old = tags->data; + + HeaderTag *new = malloc(sizeof(HeaderTag)); + new->key[0] = old->key[0]; + new->key[1] = old->key[1]; + new->value = strdup(old->value); + out->tags = list_append(out->tags, new); + + tags = tags->next; + } + return out; +} + +static int sam_header_line_merge_with(HeaderLine *out_hline, const HeaderLine *tmpl_hline) +{ + list_t *tmpl_tags; + + if ( out_hline->type[0]!=tmpl_hline->type[0] || out_hline->type[1]!=tmpl_hline->type[1] ) + return 0; + + tmpl_tags = tmpl_hline->tags; + while (tmpl_tags) + { + HeaderTag *tmpl_tag = tmpl_tags->data; + HeaderTag *out_tag = header_line_has_tag(out_hline, tmpl_tag->key); + if ( !out_tag ) + { + HeaderTag *tag = malloc(sizeof(HeaderTag)); + tag->key[0] = tmpl_tag->key[0]; + tag->key[1] = tmpl_tag->key[1]; + tag->value = strdup(tmpl_tag->value); + out_hline->tags = list_append(out_hline->tags,tag); + } + tmpl_tags = tmpl_tags->next; + } + return 1; +} + + +static HeaderLine *sam_header_line_parse(const char *headerLine) +{ + HeaderLine *hline; + HeaderTag *tag; + const char *from, *to; + from = headerLine; + + if ( *from != '@' ) { + debug("[sam_header_line_parse] expected '@', got [%s]\n", headerLine); + return 0; + } + to = ++from; + + while (*to && *to!='\t') to++; + if ( to-from != 2 ) { + debug("[sam_header_line_parse] expected '@XY', got [%s]\n", headerLine); + return 0; + } + + hline = malloc(sizeof(HeaderLine)); + hline->type[0] = from[0]; + hline->type[1] = from[1]; + hline->tags = NULL; + + int itype = tag_exists(hline->type, types); + + from = to; + while (*to && *to=='\t') to++; + if ( to-from != 1 ) { + debug("[sam_header_line_parse] multiple tabs on line [%s] (%d)\n", headerLine,(int)(to-from)); + return 0; + } + from = to; + while (*from) + { + while (*to && *to!='\t') to++; + + if ( !required_tags[itype] && !optional_tags[itype] ) + tag = new_tag(" ",from,to-1); + else + tag = new_tag(from,from+3,to-1); + + if ( header_line_has_tag(hline,tag->key) ) + debug("The tag '%c%c' present (at least) twice on line [%s]\n", tag->key[0],tag->key[1], headerLine); + hline->tags = list_append(hline->tags, tag); + + from = to; + while (*to && *to=='\t') to++; + if ( *to && to-from != 1 ) { + debug("[sam_header_line_parse] multiple tabs on line [%s] (%d)\n", headerLine,(int)(to-from)); + return 0; + } + + from = to; + } + return hline; +} + + +// Must be of an existing type, all tags must be recognised and all required tags must be present +static int sam_header_line_validate(HeaderLine *hline) +{ + list_t *tags; + HeaderTag *tag; + int itype, itag; + + // Is the type correct? + itype = tag_exists(hline->type, types); + if ( itype==-1 ) + { + debug("The type [%c%c] not recognised.\n", hline->type[0],hline->type[1]); + return 0; + } + + // Has all required tags? + itag = 0; + while ( required_tags[itype] && required_tags[itype][itag] ) + { + if ( !header_line_has_tag(hline,required_tags[itype][itag]) ) + { + debug("The tag [%c%c] required for [%c%c] not present.\n", required_tags[itype][itag][0],required_tags[itype][itag][1], + hline->type[0],hline->type[1]); + return 0; + } + itag++; + } + + // Are all tags recognised? + tags = hline->tags; + while ( tags ) + { + tag = tags->data; + if ( !tag_exists(tag->key,required_tags[itype]) && !tag_exists(tag->key,optional_tags[itype]) ) + { + debug("Unknown tag [%c%c] for [%c%c].\n", tag->key[0],tag->key[1], hline->type[0],hline->type[1]); + return 0; + } + tags = tags->next; + } + + return 1; +} + + +static void print_header_line(FILE *fp, HeaderLine *hline) +{ + list_t *tags = hline->tags; + HeaderTag *tag; + + fprintf(fp, "@%c%c", hline->type[0],hline->type[1]); + while (tags) + { + tag = tags->data; + + fprintf(fp, "\t"); + if ( tag->key[0]!=' ' || tag->key[1]!=' ' ) + fprintf(fp, "%c%c:", tag->key[0],tag->key[1]); + fprintf(fp, "%s", tag->value); + + tags = tags->next; + } + fprintf(fp,"\n"); +} + + +static void sam_header_line_free(HeaderLine *hline) +{ + list_t *tags = hline->tags; + while (tags) + { + HeaderTag *tag = tags->data; + free(tag->value); + free(tag); + tags = tags->next; + } + list_free(hline->tags); + free(hline); +} + +void sam_header_free(void *_header) +{ + HeaderDict *header = (HeaderDict*)_header; + list_t *hlines = header; + while (hlines) + { + sam_header_line_free(hlines->data); + hlines = hlines->next; + } + list_free(header); +} + +HeaderDict *sam_header_clone(const HeaderDict *dict) +{ + HeaderDict *out = NULL; + while (dict) + { + HeaderLine *hline = dict->data; + out = list_append(out, sam_header_line_clone(hline)); + dict = dict->next; + } + return out; +} + +// Returns a newly allocated string +char *sam_header_write(const void *_header) +{ + const HeaderDict *header = (const HeaderDict*)_header; + char *out = NULL; + int len=0, nout=0; + const list_t *hlines; + + // Calculate the length of the string to allocate + hlines = header; + while (hlines) + { + len += 4; // @XY and \n + + HeaderLine *hline = hlines->data; + list_t *tags = hline->tags; + while (tags) + { + HeaderTag *tag = tags->data; + len += strlen(tag->value) + 1; // \t + if ( tag->key[0]!=' ' || tag->key[1]!=' ' ) + len += strlen(tag->value) + 3; // XY: + tags = tags->next; + } + hlines = hlines->next; + } + + nout = 0; + out = malloc(len+1); + hlines = header; + while (hlines) + { + HeaderLine *hline = hlines->data; + + nout += sprintf(out+nout,"@%c%c",hline->type[0],hline->type[1]); + + list_t *tags = hline->tags; + while (tags) + { + HeaderTag *tag = tags->data; + nout += sprintf(out+nout,"\t"); + if ( tag->key[0]!=' ' || tag->key[1]!=' ' ) + nout += sprintf(out+nout,"%c%c:", tag->key[0],tag->key[1]); + nout += sprintf(out+nout,"%s", tag->value); + tags = tags->next; + } + hlines = hlines->next; + nout += sprintf(out+nout,"\n"); + } + out[len] = 0; + return out; +} + +void *sam_header_parse2(const char *headerText) +{ + list_t *hlines = NULL; + HeaderLine *hline; + const char *text; + char *buf=NULL; + size_t nbuf = 0; + + if ( !headerText ) + return 0; + + text = headerText; + while ( (text=nextline(&buf, &nbuf, text)) ) + { + hline = sam_header_line_parse(buf); + if ( hline && sam_header_line_validate(hline) ) + hlines = list_append(hlines, hline); + else + { + if (hline) sam_header_line_free(hline); + sam_header_free(hlines); + if ( buf ) free(buf); + return NULL; + } + } + if ( buf ) free(buf); + + return hlines; +} + +void *sam_header2tbl(const void *_dict, char type[2], char key_tag[2], char value_tag[2]) +{ + const HeaderDict *dict = (const HeaderDict*)_dict; + const list_t *l = dict; + khash_t(str) *tbl = kh_init(str); + khiter_t k; + int ret; + + if (_dict == 0) return tbl; // return an empty (not null) hash table + while (l) + { + HeaderLine *hline = l->data; + if ( hline->type[0]!=type[0] || hline->type[1]!=type[1] ) + { + l = l->next; + continue; + } + + HeaderTag *key, *value; + key = header_line_has_tag(hline,key_tag); + value = header_line_has_tag(hline,value_tag); + if ( !key || !value ) + { + l = l->next; + continue; + } + + k = kh_get(str, tbl, key->value); + if ( k != kh_end(tbl) ) + debug("[sam_header_lookup_table] They key %s not unique.\n", key->value); + k = kh_put(str, tbl, key->value, &ret); + kh_value(tbl, k) = value->value; + + l = l->next; + } + return tbl; +} + +char **sam_header2list(const void *_dict, char type[2], char key_tag[2], int *_n) +{ + const HeaderDict *dict = (const HeaderDict*)_dict; + const list_t *l = dict; + int max, n; + char **ret; + + ret = 0; *_n = max = n = 0; + while (l) + { + HeaderLine *hline = l->data; + if ( hline->type[0]!=type[0] || hline->type[1]!=type[1] ) + { + l = l->next; + continue; + } + + HeaderTag *key; + key = header_line_has_tag(hline,key_tag); + if ( !key ) + { + l = l->next; + continue; + } + + if (n == max) { + max = max? max<<1 : 4; + ret = realloc(ret, max * sizeof(void*)); + } + ret[n++] = key->value; + + l = l->next; + } + *_n = n; + return ret; +} + +const char *sam_tbl_get(void *h, const char *key) +{ + khash_t(str) *tbl = (khash_t(str)*)h; + khint_t k; + k = kh_get(str, tbl, key); + return k == kh_end(tbl)? 0 : kh_val(tbl, k); +} + +int sam_tbl_size(void *h) +{ + khash_t(str) *tbl = (khash_t(str)*)h; + return h? kh_size(tbl) : 0; +} + +void sam_tbl_destroy(void *h) +{ + khash_t(str) *tbl = (khash_t(str)*)h; + kh_destroy(str, tbl); +} + +void *sam_header_merge(int n, const void **_dicts) +{ + const HeaderDict **dicts = (const HeaderDict**)_dicts; + HeaderDict *out_dict; + int idict, status; + + if ( n<2 ) return NULL; + + out_dict = sam_header_clone(dicts[0]); + + for (idict=1; idictdata, out_hlines->data); + if ( status==0 ) + { + out_hlines = out_hlines->next; + continue; + } + + if ( status==2 ) + { + print_header_line(stderr,tmpl_hlines->data); + print_header_line(stderr,out_hlines->data); + debug("Conflicting lines, cannot merge the headers.\n"); + return 0; + } + if ( status==3 ) + sam_header_line_merge_with(out_hlines->data, tmpl_hlines->data); + + inserted = 1; + break; + } + if ( !inserted ) + out_dict = list_append(out_dict, sam_header_line_clone(tmpl_hlines->data)); + + tmpl_hlines = tmpl_hlines->next; + } + } + + return out_dict; +} + + diff --git a/samtools/sam_header.h b/samtools/sam_header.h new file mode 100644 index 0000000..e5c754f --- /dev/null +++ b/samtools/sam_header.h @@ -0,0 +1,24 @@ +#ifndef __SAM_HEADER_H__ +#define __SAM_HEADER_H__ + +#ifdef __cplusplus +extern "C" { +#endif + + void *sam_header_parse2(const char *headerText); + void *sam_header_merge(int n, const void **dicts); + void sam_header_free(void *header); + char *sam_header_write(const void *headerDict); // returns a newly allocated string + + char **sam_header2list(const void *_dict, char type[2], char key_tag[2], int *_n); + + void *sam_header2tbl(const void *dict, char type[2], char key_tag[2], char value_tag[2]); + const char *sam_tbl_get(void *h, const char *key); + int sam_tbl_size(void *h); + void sam_tbl_destroy(void *h); + +#ifdef __cplusplus +} +#endif + +#endif diff --git a/samtools/sam_view.c b/samtools/sam_view.c index 113c6c4..06dd01a 100644 --- a/samtools/sam_view.c +++ b/samtools/sam_view.c @@ -2,26 +2,45 @@ #include #include #include +#include +#include "sam_header.h" #include "sam.h" #include "faidx.h" static int g_min_mapQ = 0, g_flag_on = 0, g_flag_off = 0; static char *g_library, *g_rg; +static int g_sol2sanger_tbl[128]; + +static void sol2sanger(bam1_t *b) +{ + int l; + uint8_t *qual = bam1_qual(b); + if (g_sol2sanger_tbl[30] == 0) { + for (l = 0; l != 128; ++l) { + g_sol2sanger_tbl[l] = (int)(10.0 * log(1.0 + pow(10.0, (l - 64 + 33) / 10.0)) / log(10.0) + .499); + if (g_sol2sanger_tbl[l] >= 93) g_sol2sanger_tbl[l] = 93; + } + } + for (l = 0; l < b->core.l_qseq; ++l) { + int q = qual[l]; + if (q > 127) q = 127; + qual[l] = g_sol2sanger_tbl[q]; + } +} static inline int __g_skip_aln(const bam_header_t *h, const bam1_t *b) { if (b->core.qual < g_min_mapQ || ((b->core.flag & g_flag_on) != g_flag_on) || (b->core.flag & g_flag_off)) return 1; - if (g_library || g_rg) { + if (g_rg) { uint8_t *s = bam_aux_get(b, "RG"); - if (s) { - if (g_rg && strcmp(g_rg, (char*)(s + 1)) == 0) return 0; - if (g_library) { - const char *p = bam_strmap_get(h->rg2lib, (char*)(s + 1)); - return (p && strcmp(p, g_library) == 0)? 0 : 1; - } return 1; - } else return 1; - } else return 0; + if (s && strcmp(g_rg, (char*)(s + 1)) == 0) return 0; + } + if (g_library) { + const char *p = bam_get_library((bam_header_t*)h, b); + return (p && strcmp(p, g_library) == 0)? 0 : 1; + } + return 0; } // callback function for bam_fetch() @@ -36,15 +55,16 @@ static int usage(int is_long_help); int main_samview(int argc, char *argv[]) { - int c, is_header = 0, is_header_only = 0, is_bamin = 1, ret = 0, is_uncompressed = 0, is_bamout = 0; + int c, is_header = 0, is_header_only = 0, is_bamin = 1, ret = 0, is_uncompressed = 0, is_bamout = 0, slx2sngr = 0; int of_type = BAM_OFDEC, is_long_help = 0; samfile_t *in = 0, *out = 0; char in_mode[5], out_mode[5], *fn_out = 0, *fn_list = 0, *fn_ref = 0; /* parse command-line options */ strcpy(in_mode, "r"); strcpy(out_mode, "w"); - while ((c = getopt(argc, argv, "Sbt:hHo:q:f:F:ul:r:xX?T:")) >= 0) { + while ((c = getopt(argc, argv, "Sbt:hHo:q:f:F:ul:r:xX?T:C")) >= 0) { switch (c) { + case 'C': slx2sngr = 1; break; case 'S': is_bamin = 0; break; case 'b': is_bamout = 1; break; case 't': fn_list = strdup(optarg); is_bamin = 0; break; @@ -96,9 +116,12 @@ int main_samview(int argc, char *argv[]) if (argc == optind + 1) { // convert/print the entire file bam1_t *b = bam_init1(); int r; - while ((r = samread(in, b)) >= 0) // read one alignment from `in' - if (!__g_skip_aln(in->header, b)) + while ((r = samread(in, b)) >= 0) { // read one alignment from `in' + if (!__g_skip_aln(in->header, b)) { + if (slx2sngr) sol2sanger(b); samwrite(out, b); // write the alignment to `out' + } + } if (r < -1) fprintf(stderr, "[main_samview] truncated file.\n"); bam_destroy1(b); } else { // retrieve alignments in specified regions diff --git a/setup.py b/setup.py index 57b19e7..098cb7f 100644 --- a/setup.py +++ b/setup.py @@ -9,7 +9,7 @@ pysam import os, sys, glob, shutil name = "pysam" -version = "0.1.2" +version = "0.2" samtools_exclude = ( "bamtk.c", "razip.c", "bgzip.c" ) samtools_dest = os.path.abspath( "samtools" ) @@ -71,7 +71,7 @@ metadata = { 'platforms': "ALL", 'url': "http://code.google.com/p/pysam/", 'py_modules': [ - "pysam/__init__", "pysam/Pileup"], + "pysam/__init__", "pysam/Pileup", "pysam/namedtuple" ], 'ext_modules': [pysam,], 'cmdclass' : {'build_ext': build_ext} } diff --git a/tests/Makefile b/tests/Makefile index 5ce6e7e..5403750 100644 --- a/tests/Makefile +++ b/tests/Makefile @@ -1,17 +1,16 @@ -all: ex1.glf ex1.pileup.gz ex1.bam.bai ex1.glfview.gz ex2.sam.gz ex2.sam ex1.sam ex3.bam ex3.bam.bai ex4.bam ex4.bam.bai ex5.bam ex5.bam.bai - @echo; echo \# You can now launch the viewer with: \'samtools tview ex1.bam ex1.fa\'; echo; +all: ex1.glf ex1.pileup.gz ex1.bam.bai ex1.glfview.gz \ + ex2.sam.gz ex2.sam ex1.sam \ + ex2.bam \ + ex3.bam ex3.bam.bai \ + ex4.bam ex4.bam.bai \ + ex5.bam ex5.bam.bai \ + ex6.bam ex2.sam.gz: ex1.bam ex1.bam.bai samtools view -h ex1.bam | gzip > ex2.sam.gz -ex3.bam: ex3.sam ex1.fa.fai - samtools import ex1.fa.fai ex3.sam ex3.bam - -ex4.bam: ex4.sam ex1.fa.fai - samtools import ex1.fa.fai ex4.sam ex4.bam - -ex5.bam: ex5.sam ex1.fa.fai - samtools import ex1.fa.fai ex5.sam ex5.bam +%.bam: %.sam ex1.fa.fai + samtools import ex1.fa.fai $< $@ %.sam: %.sam.gz gunzip < $< > $@ @@ -24,7 +23,7 @@ ex1.bam:ex1.sam.gz ex1.fa.fai samtools index $< ex1.pileup.gz:ex1.bam ex1.fa samtools pileup -cf ex1.fa ex1.bam | gzip > ex1.pileup.gz -ex1.glf:ex1.bam ex1.fa +ex1.glf:ex1.bam ex1.fa samtools pileup -gf ex1.fa ex1.bam > ex1.glf ex1.glfview.gz:ex1.glf samtools glfview ex1.glf | gzip > ex1.glfview.gz diff --git a/tests/ex3.sam b/tests/ex3.sam index 801d13a..bae2a22 100644 --- a/tests/ex3.sam +++ b/tests/ex3.sam @@ -1,9 +1,13 @@ @HD VN:1.0 @SQ SN:chr1 LN:1575 @SQ SN:chr2 LN:1584 -@RG ID:L1 PU:SC_1_10 LB:SC_1 SM:NA12891 -@RG ID:L2 PU:SC_2_12 LB:SC_2 SM:NA12891 +@RG ID:L1 PU:SC_1_10 LB:SC_1 SM:NA12891 CN:name:with:colon +@RG ID:L2 PU:SC_2_12 LB:SC_2 SM:NA12891 CN:name:with:colon +@PG ID:P1 VN:1.0 +@PG ID:P2 VN:1.1 @CO this is a comment @CO this is another comment -read_28833_29006_6945 99 chr1 33 20 10M1D25M = 200 167 AGCTTAGCTAGCTACCTATATCTTGGTCTTGGCCG <<<<<<<<<<<<<<<<<<<<<:<9/,&,22;;<<< NM:i:1 RG:Z:L1 -read_28701_28881_323b 147 chr2 88 30 35M = 500 412 ACCTATATCTTGGCCTTGGCCGATGCGGCCTTGCA <<<<<;<<<<7;:<<<6;<<<<<<<<<<<<7<<<< MF:i:18 RG:Z:L2 +read_28833_29006_6945 99 chr1 33 20 10M1D25M = 200 167 AGCTTAGCTAGCTACCTATATCTTGGTCTTGGCCG <<<<<<<<<<<<<<<<<<<<<:<9/,&,22;;<<< NM:i:1 RG:Z:L1 PG:Z:P1 XT:A:U +read_28701_28881_323b 147 chr2 88 30 35M = 500 412 ACCTATATCTTGGCCTTGGCCGATGCGGCCTTGCA <<<<<;<<<<7;:<<<6;<<<<<<<<<<<<7<<<< MF:i:18 RG:Z:L2 PG:Z:P2 XT:A:R +read_28701_28881_323c 147 chr2 88 30 35M = 500 412 ACCTATATCTTGGCCTTGGCCGATGCGGCCTTGCA <<<<<;<<<<7;:<<<6;<<<<<<<<<<<<7<<<< + diff --git a/tests/ex5.sam b/tests/ex5.sam new file mode 100644 index 0000000..f1f8aad --- /dev/null +++ b/tests/ex5.sam @@ -0,0 +1,5 @@ +@HD VN:1.0 +@SQ SN:chr1 LN:100 +@SQ SN:chr2 LN:100 +read_28833_29006_6945 0 * * * * * 200 167 AGCTTAGCTAGCTACCTATATCTTGGTCTTGGCCG <<<<<<<<<<<<<<<<<<<<<:<9/,&,22;;<<< +read_28701_28881_323b 0 * * * * * 500 412 ACCTATATCTTGGCCTTGGCCGATGCGGCCTTGCA <<<<<;<<<<7;:<<<6;<<<<<<<<<<<<7<<<< diff --git a/tests/ex6.sam b/tests/ex6.sam new file mode 100644 index 0000000..7ae90f3 --- /dev/null +++ b/tests/ex6.sam @@ -0,0 +1,5 @@ +@HD VN:1.0 +@SQ SN:chr1 LN:1575 +@SQ SN:chr2 LN:1584 +read_28833_29006_6945 99 chr1 33 20 10M1D25M = 200 167 AGCTTAGCTAGCTACCTATATCTTGGTCTTGGCCG <<<<<<<<<<<<<<<<<<<<<:<9/,&,22;;<<< NM:i:1 RG:Z:L1 +read_28701_28881_323b 147 chr2 88 30 35M = 500 412 ACCTATATCTTGGCCTTGGCCGATGCGGCCTTGCA <<<<<;<<<<7;:<<<6;<<<<<<<<<<<<7<<<< MF:i:18 RG:Z:L2 diff --git a/tests/ex7.sam b/tests/ex7.sam new file mode 100644 index 0000000..12befae --- /dev/null +++ b/tests/ex7.sam @@ -0,0 +1,2 @@ +read_28833_29006_6945 99 chr1 33 20 10M1D25M = 200 167 AGCTTAGCTAGCTACCTATATCTTGGTCTTGGCCG <<<<<<<<<<<<<<<<<<<<<:<9/,&,22;;<<< NM:i:1 RG:Z:L1 PG:Z:P1 XT:A:U +read_28701_28881_323b 147 chr2 88 30 35M = 500 412 ACCTATATCTTGGCCTTGGCCGATGCGGCCTTGCA <<<<<;<<<<7;:<<<6;<<<<<<<<<<<<7<<<< MF:i:18 RG:Z:L2 PG:Z:P2 XT:A:R diff --git a/tests/example.py b/tests/example.py index 42e7347..a1ca7a0 100644 --- a/tests/example.py +++ b/tests/example.py @@ -1,3 +1,4 @@ +import sys import pysam samfile = pysam.Samfile( "ex1.bam", "rb" ) @@ -5,40 +6,40 @@ samfile = pysam.Samfile( "ex1.bam", "rb" ) print "###################" # check different ways to iterate print len(list(samfile.fetch())) -print len(list(samfile.fetch( "seq1", 10, 200 ))) -print len(list(samfile.fetch( region="seq1:10-200" ))) -print len(list(samfile.fetch( "seq1" ))) -print len(list(samfile.fetch( region="seq1"))) -print len(list(samfile.fetch( "seq2" ))) -print len(list(samfile.fetch( region="seq2"))) +print len(list(samfile.fetch( "chr1", 10, 200 ))) +print len(list(samfile.fetch( region="chr1:10-200" ))) +print len(list(samfile.fetch( "chr1" ))) +print len(list(samfile.fetch( region="chr1"))) +print len(list(samfile.fetch( "chr2" ))) +print len(list(samfile.fetch( region="chr2"))) print len(list(samfile.fetch())) -print len(list(samfile.fetch( "seq1" ))) -print len(list(samfile.fetch( region="seq1"))) +print len(list(samfile.fetch( "chr1" ))) +print len(list(samfile.fetch( region="chr1"))) print len(list(samfile.fetch())) -print len(list(samfile.pileup( "seq1", 10, 200 ))) -print len(list(samfile.pileup( region="seq1:10-200" ))) -print len(list(samfile.pileup( "seq1" ))) -print len(list(samfile.pileup( region="seq1"))) -print len(list(samfile.pileup( "seq2" ))) -print len(list(samfile.pileup( region="seq2"))) +print len(list(samfile.pileup( "chr1", 10, 200 ))) +print len(list(samfile.pileup( region="chr1:10-200" ))) +print len(list(samfile.pileup( "chr1" ))) +print len(list(samfile.pileup( region="chr1"))) +print len(list(samfile.pileup( "chr2" ))) +print len(list(samfile.pileup( region="chr2"))) print len(list(samfile.pileup())) print len(list(samfile.pileup())) print "########### fetch with callback ################" def my_fetch_callback( alignment ): print str(alignment) -samfile.fetch( region="seq1:10-200", callback=my_fetch_callback ) +samfile.fetch( region="chr1:10-200", callback=my_fetch_callback ) print "########## pileup with callback ################" -def my_pileup_callback( pileups ): print str(pileups) -samfile.pileup( region="seq1:10-200", callback=my_pileup_callback ) +def my_pileup_callback( column ): print str(column) +samfile.pileup( region="chr1:10-200", callback=my_pileup_callback ) print "##########iterator row #################" -iter = pysam.IteratorRow( samfile, "seq1:10-200") +iter = pysam.IteratorRow( samfile, 0, 10, 200) for x in iter: print str(x) print "##########iterator col #################" -iter = pysam.IteratorColumn( samfile, "seq1:10-200") +iter = pysam.IteratorColumn( samfile, 0, 10, 200 ) for x in iter: print str(x) print "#########row all##################" @@ -54,9 +55,10 @@ class Counter: self.mCounts += 1 c = Counter() -samfile.fetch( "seq1:10-200", c ) +samfile.fetch( "chr1:10-200", c ) print "counts=", c.mCounts +sys.exit(0) print samfile.getTarget( 0 ) print samfile.getTarget( 1 ) @@ -87,7 +89,7 @@ def my_fetch_callback( alignment ): print str(alignment) try: - samfile.fetch( "seq1:10-20", my_fetch_callback ) + samfile.fetch( "chr1:10-20", my_fetch_callback ) except AssertionError: print "caught fetch exception" @@ -98,7 +100,7 @@ samfile = pysam.Samfile( "ex2.sam.gz", "r" ) def my_pileup_callback( pileups ): print str(pileups) try: - samfile.pileup( "seq1:10-20", my_pileup_callback ) + samfile.pileup( "chr1:10-20", my_pileup_callback ) except NotImplementedError: print "caught pileup exception" diff --git a/tests/pysam_test.py b/tests/pysam_test.py index 6d27cc6..c2ae6fa 100755 --- a/tests/pysam_test.py +++ b/tests/pysam_test.py @@ -1,5 +1,9 @@ #!/usr/bin/env python -'''unit testing code for pysam.''' +'''unit testing code for pysam. + +Execute in the :file:`tests` directory as it requires the Makefile +and data files located there. +''' import pysam import unittest @@ -8,6 +12,7 @@ import itertools import subprocess import shutil + def checkBinaryEqual( filename1, filename2 ): '''return true if the two files are binary equal.''' if os.path.getsize( filename1 ) != os.path.getsize( filename2 ): @@ -42,6 +47,7 @@ def runSamtools( cmd ): except OSError, e: print >>sys.stderr, "Execution failed:", e + class BinaryTest(unittest.TestCase): '''test samtools command line commands and compare against pysam commands. @@ -222,6 +228,32 @@ class IOTest(unittest.TestCase): self.checkEcho( input_filename, reference_filename, output_filename, "r", "w" ) + def testFetchFromClosedFile( self ): + + samfile = pysam.Samfile( "ex1.bam", "rb" ) + samfile.close() + self.assertRaises( ValueError, samfile.fetch, 'chr1', 100, 120) + + def testPileupFromClosedFile( self ): + + samfile = pysam.Samfile( "ex1.bam", "rb" ) + samfile.close() + self.assertRaises( ValueError, samfile.pileup, 'chr1', 100, 120) + + def testBinaryReadFromSamfile( self ): + pass + # needs to re-activated, see issue 19 + #samfile = pysam.Samfile( "ex1.bam", "r" ) + #samfile.fetch().next() + + def testReadingFromFileWithoutIndex( self ): + '''read from bam file without index.''' + + assert not os.path.exists( "ex2.bam.bai" ) + samfile = pysam.Samfile( "ex2.bam", "rb" ) + self.assertRaises( ValueError, samfile.fetch ) + self.assertEqual( len(list( samfile.fetch(until_eof = True) )), 3270 ) + class TestIteratorRow(unittest.TestCase): def setUp(self): @@ -321,7 +353,7 @@ class TestIteratorColumn(unittest.TestCase): def tearDown(self): self.samfile.close() -class TestAlignedReadBam(unittest.TestCase): +class TestAlignedReadFromBam(unittest.TestCase): def setUp(self): self.samfile=pysam.Samfile( "ex3.bam","rb" ) @@ -383,96 +415,72 @@ class TestAlignedReadBam(unittest.TestCase): self.assertEqual( self.reads[0].is_proper_pair, True, "is proper pair mismatch in read 1: %s != %s" % (self.reads[0].is_proper_pair, True) ) self.assertEqual( self.reads[1].is_proper_pair, True, "is proper pair mismatch in read 2: %s != %s" % (self.reads[1].is_proper_pair, True) ) - def tearDown(self): - self.samfile.close() - -class TestAlignedReadSam(unittest.TestCase): - - def setUp(self): - self.samfile=pysam.Samfile( "ex3.sam","r" ) - self.reads=list(self.samfile.fetch()) - - def testARqname(self): - self.assertEqual( self.reads[0].qname, "read_28833_29006_6945", "read name mismatch in read 1: %s != %s" % (self.reads[0].qname, "read_28833_29006_6945") ) - self.assertEqual( self.reads[1].qname, "read_28701_28881_323b", "read name mismatch in read 2: %s != %s" % (self.reads[1].qname, "read_28701_28881_323b") ) - - def testARflag(self): - self.assertEqual( self.reads[0].flag, 99, "flag mismatch in read 1: %s != %s" % (self.reads[0].flag, 99) ) - self.assertEqual( self.reads[1].flag, 147, "flag mismatch in read 2: %s != %s" % (self.reads[1].flag, 147) ) - - def testARrname(self): - self.assertEqual( self.reads[0].rname, 0, "chromosome/target id mismatch in read 1: %s != %s" % (self.reads[0].rname, 0) ) - self.assertEqual( self.reads[1].rname, 1, "chromosome/target id mismatch in read 2: %s != %s" % (self.reads[1].rname, 1) ) - - def testARpos(self): - self.assertEqual( self.reads[0].pos, 33-1, "mapping position mismatch in read 1: %s != %s" % (self.reads[0].pos, 33-1) ) - self.assertEqual( self.reads[1].pos, 88-1, "mapping position mismatch in read 2: %s != %s" % (self.reads[1].pos, 88-1) ) + def testTags( self ): + self.assertEqual( self.reads[0].tags, + [('NM', 1), ('RG', 'L1'), + ('PG', 'P1'), ('XT', 'U')] ) + self.assertEqual( self.reads[1].tags, + [('MF', 18), ('RG', 'L2'), + ('PG', 'P2'),('XT', 'R') ] ) - def testARmapq(self): - self.assertEqual( self.reads[0].mapq, 20, "mapping quality mismatch in read 1: %s != %s" % (self.reads[0].mapq, 20) ) - self.assertEqual( self.reads[1].mapq, 30, "mapping quality mismatch in read 2: %s != %s" % (self.reads[1].mapq, 30) ) + def testOpt( self ): + self.assertEqual( self.reads[0].opt("XT"), "U" ) + self.assertEqual( self.reads[1].opt("XT"), "R" ) - def testARcigar(self): - self.assertEqual( self.reads[0].cigar, [(0, 10), (2, 1), (0, 25)], "read name length mismatch in read 1: %s != %s" % (self.reads[0].cigar, [(0, 10), (2, 1), (0, 25)]) ) - self.assertEqual( self.reads[1].cigar, [(0, 35)], "read name length mismatch in read 2: %s != %s" % (self.reads[1].cigar, [(0, 35)]) ) + def testMissingOpt( self ): + self.assertRaises( KeyError, self.reads[0].opt, "XP" ) - def testARmrnm(self): - self.assertEqual( self.reads[0].mrnm, 0, "mate reference sequence name mismatch in read 1: %s != %s" % (self.reads[0].mrnm, 0) ) - self.assertEqual( self.reads[1].mrnm, 1, "mate reference sequence name mismatch in read 2: %s != %s" % (self.reads[1].mrnm, 1) ) - - def testARmpos(self): - self.assertEqual( self.reads[0].mpos, 200-1, "mate mapping position mismatch in read 1: %s != %s" % (self.reads[0].mpos, 200-1) ) - self.assertEqual( self.reads[1].mpos, 500-1, "mate mapping position mismatch in read 2: %s != %s" % (self.reads[1].mpos, 500-1) ) - - def testARisize(self): - self.assertEqual( self.reads[0].isize, 167, "insert size mismatch in read 1: %s != %s" % (self.reads[0].isize, 167) ) - self.assertEqual( self.reads[1].isize, 412, "insert size mismatch in read 2: %s != %s" % (self.reads[1].isize, 412) ) - - def testARseq(self): - self.assertEqual( self.reads[0].seq, "AGCTTAGCTAGCTACCTATATCTTGGTCTTGGCCG", "sequence mismatch in read 1: %s != %s" % (self.reads[0].seq, "AGCTTAGCTAGCTACCTATATCTTGGTCTTGGCCG") ) - self.assertEqual( self.reads[1].seq, "ACCTATATCTTGGCCTTGGCCGATGCGGCCTTGCA", "sequence size mismatch in read 2: %s != %s" % (self.reads[1].seq, "ACCTATATCTTGGCCTTGGCCGATGCGGCCTTGCA") ) - - def testARqual(self): - self.assertEqual( self.reads[0].qual, "<<<<<<<<<<<<<<<<<<<<<:<9/,&,22;;<<<", "quality string mismatch in read 1: %s != %s" % (self.reads[0].qual, "<<<<<<<<<<<<<<<<<<<<<:<9/,&,22;;<<<") ) - self.assertEqual( self.reads[1].qual, "<<<<<;<<<<7;:<<<6;<<<<<<<<<<<<7<<<<", "quality string mismatch in read 2: %s != %s" % (self.reads[1].qual, "<<<<<;<<<<7;:<<<6;<<<<<<<<<<<<7<<<<") ) - - def testPresentOptionalFields(self): - self.assertEqual( self.reads[0].opt('NM'), 1, "optional field mismatch in read 1, NM: %s != %s" % (self.reads[0].opt('NM'), 1) ) - self.assertEqual( self.reads[0].opt('RG'), 'L1', "optional field mismatch in read 1, RG: %s != %s" % (self.reads[0].opt('RG'), 'L1') ) - self.assertEqual( self.reads[1].opt('RG'), 'L2', "optional field mismatch in read 2, RG: %s != %s" % (self.reads[1].opt('RG'), 'L2') ) - self.assertEqual( self.reads[1].opt('MF'), 18, "optional field mismatch in read 2, MF: %s != %s" % (self.reads[1].opt('MF'), 18) ) - - def testPairedBools(self): - self.assertEqual( self.reads[0].is_paired, True, "is paired mismatch in read 1: %s != %s" % (self.reads[0].is_paired, True) ) - self.assertEqual( self.reads[1].is_paired, True, "is paired mismatch in read 2: %s != %s" % (self.reads[1].is_paired, True) ) - self.assertEqual( self.reads[0].is_proper_pair, True, "is proper pair mismatch in read 1: %s != %s" % (self.reads[0].is_proper_pair, True) ) - self.assertEqual( self.reads[1].is_proper_pair, True, "is proper pair mismatch in read 2: %s != %s" % (self.reads[1].is_proper_pair, True) ) + def testEmptyOpt( self ): + self.assertRaises( KeyError, self.reads[2].opt, "XT" ) def tearDown(self): self.samfile.close() -class TestHeaderSam(unittest.TestCase): +class TestAlignedReadFromSam(TestAlignedReadFromBam): def setUp(self): self.samfile=pysam.Samfile( "ex3.sam","r" ) + self.reads=list(self.samfile.fetch()) - def testHeaders(self): - self.assertEqual( self.samfile.header, {'SQ': [{'LN': 1575, 'SN': 'chr1'}, {'LN': 1584, 'SN': 'chr2'}], 'RG': [{'LB': 'SC_1', 'ID': 'L1', 'SM': 'NA12891', 'PU': 'SC_1_10'}, {'LB': 'SC_2', 'ID': 'L2', 'SM': 'NA12891', 'PU': 'SC_2_12'}], 'CO': ['this is a comment', 'this is another comment'], 'HD': {'VN': '1.0'}}, "mismatch in headers: %s != %s" % (self.samfile.header, {'SQ': [{'LN': 1575, 'SN': 'chr1'}, {'LN': 1584, 'SN': 'chr2'}], 'RG': [{'LB': 'SC_1', 'ID': 'L1', 'SM': 'NA12891', 'PU': 'SC_1_10'}, {'LB': 'SC_2', 'ID': 'L2', 'SM': 'NA12891', 'PU': 'SC_2_12'}], 'CO': ['this is a comment', 'this is another comment'], 'HD': {'VN': '1.0'}}) ) +# needs to be implemented +# class TestAlignedReadFromSamWithoutHeader(TestAlignedReadFromBam): +# +# def setUp(self): +# self.samfile=pysam.Samfile( "ex7.sam","r" ) +# self.reads=list(self.samfile.fetch()) - def tearDown(self): - self.samfile.close() +class TestHeaderSam(unittest.TestCase): -class TestHeaderBam(unittest.TestCase): + header = {'SQ': [{'LN': 1575, 'SN': 'chr1'}, + {'LN': 1584, 'SN': 'chr2'}], + 'RG': [{'LB': 'SC_1', 'ID': 'L1', 'SM': 'NA12891', 'PU': 'SC_1_10', "CN":"name:with:colon"}, + {'LB': 'SC_2', 'ID': 'L2', 'SM': 'NA12891', 'PU': 'SC_2_12', "CN":"name:with:colon"}], + 'PG': [{'ID': 'P1', 'VN': '1.0'}, {'ID': 'P2', 'VN': '1.1'}], + 'HD': {'VN': '1.0'}, + 'CO' : [ 'this is a comment', 'this is another comment'], + } + + def compareHeaders( self, a, b ): + '''compare two headers a and b.''' + for ak,av in a.iteritems(): + self.assertTrue( ak in b, "key '%s' not in '%s' " % (ak,b) ) + self.assertEqual( av, b[ak] ) def setUp(self): self.samfile=pysam.Samfile( "ex3.sam","r" ) def testHeaders(self): - self.assertEqual( self.samfile.header, {'SQ': [{'LN': 1575, 'SN': 'chr1'}, {'LN': 1584, 'SN': 'chr2'}], 'RG': [{'LB': 'SC_1', 'ID': 'L1', 'SM': 'NA12891', 'PU': 'SC_1_10'}, {'LB': 'SC_2', 'ID': 'L2', 'SM': 'NA12891', 'PU': 'SC_2_12'}], 'CO': ['this is a comment', 'this is another comment'], 'HD': {'VN': '1.0'}}, "mismatch in headers: %s != %s" % (self.samfile.header, {'SQ': [{'LN': 1575, 'SN': 'chr1'}, {'LN': 1584, 'SN': 'chr2'}], 'RG': [{'LB': 'SC_1', 'ID': 'L1', 'SM': 'NA12891', 'PU': 'SC_1_10'}, {'LB': 'SC_2', 'ID': 'L2', 'SM': 'NA12891', 'PU': 'SC_2_12'}], 'CO': ['this is a comment', 'this is another comment'], 'HD': {'VN': '1.0'}}) ) - + self.compareHeaders( self.header, self.samfile.header ) + self.compareHeaders( self.samfile.header, self.header ) + def tearDown(self): self.samfile.close() +class TestHeaderBam(TestHeaderSam): + + def setUp(self): + self.samfile=pysam.Samfile( "ex3.bam","rb" ) + class TestUnmappedReads(unittest.TestCase): def testSAM(self): @@ -510,17 +518,18 @@ class TestPileupObjects(unittest.TestCase): def tearDown(self): self.samfile.close() - + class TestExceptions(unittest.TestCase): def setUp(self): self.samfile=pysam.Samfile( "ex1.bam","rb" ) - def testBadFile(self): + def testMissingFile(self): + self.assertRaises( IOError, pysam.Samfile, "exdoesntexist.bam", "rb" ) - self.assertRaises( IOError, pysam.Samfile, "exdoesntexist.sam","r" ) + self.assertRaises( IOError, pysam.Samfile, "exdoesntexist.sam", "r" ) self.assertRaises( IOError, pysam.Samfile, "exdoesntexist.bam", "r" ) - self.assertRaises( IOError, pysam.Samfile, "exdoesntexist.sam","rb" ) + self.assertRaises( IOError, pysam.Samfile, "exdoesntexist.sam", "rb" ) def testBadContig(self): self.assertRaises( ValueError, self.samfile.fetch, "chr88" ) @@ -556,9 +565,277 @@ class TestExceptions(unittest.TestCase): def tearDown(self): self.samfile.close() +class TestFastaFile(unittest.TestCase): + + mSequences = { 'chr1' : + "CACTAGTGGCTCATTGTAAATGTGTGGTTTAACTCGTCCATGGCCCAGCATTAGGGAGCTGTGGACCCTGCAGCCTGGCTGTGGGGGCCGCAGTGGCTGAGGGGTGCAGAGCCGAGTCACGGGGTTGCCAGCACAGGGGCTTAACCTCTGGTGACTGCCAGAGCTGCTGGCAAGCTAGAGTCCCATTTGGAGCCCCTCTAAGCCGTTCTATTTGTAATGAAAACTATATTTATGCTATTCAGTTCTAAATATAGAAATTGAAACAGCTGTGTTTAGTGCCTTTGTTCAACCCCCTTGCAACAACCTTGAGAACCCCAGGGAATTTGTCAATGTCAGGGAAGGAGCATTTTGTCAGTTACCAAATGTGTTTATTACCAGAGGGATGGAGGGAAGAGGGACGCTGAAGAACTTTGATGCCCTCTTCTTCCAAAGATGAAACGCGTAACTGCGCTCTCATTCACTCCAGCTCCCTGTCACCCAATGGACCTGTGATATCTGGATTCTGGGAAATTCTTCATCCTGGACCCTGAGAGATTCTGCAGCCCAGCTCCAGATTGCTTGTGGTCTGACAGGCTGCAACTGTGAGCCATCACAATGAACAACAGGAAGAAAAGGTCTTTCAAAAGGTGATGTGTGTTCTCATCAACCTCATACACACACATGGTTTAGGGGTATAATACCTCTACATGGCTGATTATGAAAACAATGTTCCCCAGATACCATCCCTGTCTTACTTCCAGCTCCCCAGAGGGAAAGCTTTCAACGCTTCTAGCCATTTCTTTTGGCATTTGCCTTCAGACCCTACACGAATGCGTCTCTACCACAGGGGGCTGCGCGGTTTCCCATCATGAAGCACTGAACTTCCACGTCTCATCTAGGGGAACAGGGAGGTGCACTAATGCGCTCCACGCCCAAGCCCTTCTCACAGTTTCTGCCCCCAGCATGGTTGTACTGGGCAATACATGAGATTATTAGGAAATGCTTTACTGTCATAACTATGAAGAGACTATTGCCAGATGAACCACACATTAATACTATGTTTCTTATCTGCACATTACTACCCTGCAATTAATATAATTGTGTCCATGTACACACGCTGTCCTATGTACTTATCATGACTCTATCCCAAATTCCCAATTACGTCCTATCTTCTTCTTAGGGAAGAACAGCTTAGGTATCAATTTGGTGTTCTGTGTAAAGTCTCAGGGAGCCGTCCGTGTCCTCCCATCTGGCCTCGTCCACACTGGTTCTCTTGAAAGCTTGGGCTGTAATGATGCCCCTTGGCCATCACCCAGTCCCTGCCCCATCTCTTGTAATCTCTCTCCTTTTTGCTGCATCCCTGTCTTCCTCTGTCTTGATTTACTTGTTGTTGGTTTTCTGTTTCTTTGTTTGATTTGGTGGAAGACATAATCCCACGCTTCCTATGGAAAGGTTGTTGGGAGATTTTTAATGATTCCTCAATGTTAAAATGTCTATTTTTGTCTTGACACCCAACTAATATTTGTCTGAGCAAAACAGTCTAGATGAGAGAGAACTTCCCTGGAGGTCTGATGGCGTTTCTCCCTCGTCTTCTTA", + 'chr2' : + "TTCAAATGAACTTCTGTAATTGAAAAATTCATTTAAGAAATTACAAAATATAGTTGAAAGCTCTAACAATAGACTAAACCAAGCAGAAGAAAGAGGTTCAGAACTTGAAGACAAGTCTCTTATGAATTAACCCAGTCAGACAAAAATAAAGAAAAAAATTTTAAAAATGAACAGAGCTTTCAAGAAGTATGAGATTATGTAAAGTAACTGAACCTATGAGTCACAGGTATTCCTGAGGAAAAAGAAAAAGTGAGAAGTTTGGAAAAACTATTTGAGGAAGTAATTGGGGAAAACCTCTTTAGTCTTGCTAGAGATTTAGACATCTAAATGAAAGAGGCTCAAAGAATGCCAGGAAGATACATTGCAAGACAGACTTCATCAAGATATGTAGTCATCAGACTATCTAAAGTCAACATGAAGGAAAAAAATTCTAAAATCAGCAAGAGAAAAGCATACAGTCATCTATAAAGGAAATCCCATCAGAATAACAATGGGCTTCTCAGCAGAAACCTTACAAGCCAGAAGAGATTGGATCTAATTTTTGGACTTCTTAAAGAAAAAAAAACCTGTCAAACACGAATGTTATGCCCTGCTAAACTAAGCATCATAAATGAAGGGGAAATAAAGTCAAGTCTTTCCTGACAAGCAAATGCTAAGATAATTCATCATCACTAAACCAGTCCTATAAGAAATGCTCAAAAGAATTGTAAAAGTCAAAATTAAAGTTCAATACTCACCATCATAAATACACACAAAAGTACAAAACTCACAGGTTTTATAAAACAATTGAGACTACAGAGCAACTAGGTAAAAAATTAACATTACAACAGGAACAAAACCTCATATATCAATATTAACTTTGAATAAAAAGGGATTAAATTCCCCCACTTAAGAGATATAGATTGGCAGAACAGATTTAAAAACATGAACTAACTATATGCTGTTTACAAGAAACTCATTAATAAAGACATGAGTTCAGGTAAAGGGGTGGAAAAAGATGTTCTACGCAAACAGAAACCAAATGAGAGAAGGAGTAGCTATACTTATATCAGATAAAGCACACTTTAAATCAACAACAGTAAAATAAAACAAAGGAGGTCATCATACAATGATAAAAAGATCAATTCAGCAAGAAGATATAACCATCCTACTAAATACATATGCACCTAACACAAGACTACCCAGATTCATAAAACAAATACTACTAGACCTAAGAGGGATGAGAAATTACCTAATTGGTACAATGTACAATATTCTGATGATGGTTACACTAAAAGCCCATACTTTACTGCTACTCAATATATCCATGTAACAAATCTGCGCTTGTACTTCTAAATCTATAAAAAAATTAAAATTTAACAAAAGTAAATAAAACACATAGCTAAAACTAAAAAAGCAAAAACAAAAACTATGCTAAGTATTGGTAAAGATGTGGGGAAAAAAGTAAACTCTCAAATATTGCTAGTGGGAGTATAAATTGTTTTCCACTTTGGAAAACAATTTGGTAATTTCGTTTTTTTTTTTTTCTTTTCTCTTTTTTTTTTTTTTTTTTTTGCATGCCAGAAAAAAATATTTACAGTAACT", + } + + def setUp(self): + self.file=pysam.Fastafile( "ex1.fa" ) + + def testFetch(self): + for id, seq in self.mSequences.items(): + self.assertEqual( seq, self.file.fetch( id ) ) + for x in range( 0, len(seq), 10): + self.assertEqual( seq[x:x+10], self.file.fetch( id, x, x+10) ) + + def testFetchErrors( self ): + self.assertRaises( ValueError, self.file.fetch ) + self.assertRaises( ValueError, self.file.fetch, "chr1", 0 ) + self.assertRaises( ValueError, self.file.fetch, "chr1", -1, 10 ) + self.assertRaises( ValueError, self.file.fetch, "chr1", 20, 10 ) + # the following segfaults: + # self.assertRaises( IndexError, self.file.fetch, "chr12", ) + pass + + def tearDown(self): + self.file.close() + + +class TestAlignedRead(unittest.TestCase): + '''tests to check if aligned read can be constructed + and manipulated. + ''' + + def checkFieldEqual( self, read1, read2, exclude = []): + '''check if two reads are equal by comparing each field.''' + + for x in ("qname", "seq", "flag", + "rname", "pos", "mapq", "cigar", + "mrnm", "mpos", "isize", "qual", + "is_paired", "is_proper_pair", + "is_unmapped", "mate_is_unmapped", + "is_reverse", "mate_is_reverse", + "is_read1", "is_read2", + "is_secondary", "is_qcfail", + "is_duplicate", "bin"): + if x in exclude: continue + self.assertEqual( getattr(read1, x), getattr(read2,x), "attribute mismatch for %s: %s != %s" % + (x, getattr(read1, x), getattr(read2,x))) + + def testEmpty( self ): + a = pysam.AlignedRead() + self.assertEqual( a.qname, None ) + self.assertEqual( a.seq, None ) + self.assertEqual( a.qual, None ) + self.assertEqual( a.flag, 0 ) + self.assertEqual( a.rname, 0 ) + self.assertEqual( a.mapq, 0 ) + self.assertEqual( a.cigar, None ) + self.assertEqual( a.tags, None ) + self.assertEqual( a.mrnm, 0 ) + self.assertEqual( a.mpos, 0 ) + self.assertEqual( a.isize, 0 ) + + def buildRead( self ): + '''build an example read.''' + + a = pysam.AlignedRead() + a.qname = "read_12345" + a.seq="ACGT" * 3 + a.flag = 0 + a.rname = 0 + a.pos = 33 + a.mapq = 20 + a.cigar = ( (0,10), (2,1), (0,25) ) + a.mrnm = 0 + a.mpos=200 + a.isize=167 + a.qual="1234" * 3 + + return a + + def testUpdate( self ): + '''check if updating fields affects other variable length data + ''' + a = self.buildRead() + b = self.buildRead() + + # check qname + b.qname = "read_123" + self.checkFieldEqual( a, b, "qname" ) + b.qname = "read_12345678" + self.checkFieldEqual( a, b, "qname" ) + b.qname = "read_12345" + self.checkFieldEqual( a, b) + + # check cigar + b.cigar = ( (0,10), ) + self.checkFieldEqual( a, b, "cigar" ) + b.cigar = ( (0,10), (2,1), (0,25), (2,1), (0,25) ) + self.checkFieldEqual( a, b, "cigar" ) + b.cigar = ( (0,10), (2,1), (0,25) ) + self.checkFieldEqual( a, b) + + # check seq + b.seq = "ACGT" + self.checkFieldEqual( a, b, ("seq", "qual") ) + b.seq = "ACGT" * 10 + self.checkFieldEqual( a, b, ("seq", "qual") ) + b.seq = "ACGT" * 3 + self.checkFieldEqual( a, b, ("qual",)) + + # reset qual + b = self.buildRead() + + # check flags: + for x in ( + "is_paired", "is_proper_pair", + "is_unmapped", "mate_is_unmapped", + "is_reverse", "mate_is_reverse", + "is_read1", "is_read2", + "is_secondary", "is_qcfail", + "is_duplicate"): + setattr( b, x, True ) + self.assertEqual( getattr(b, x), True ) + self.checkFieldEqual( a, b, ("flag", x,) ) + setattr( b, x, False ) + self.assertEqual( getattr(b, x), False ) + self.checkFieldEqual( a, b ) + + def testLargeRead( self ): + '''build an example read.''' + + a = pysam.AlignedRead() + a.qname = "read_12345" + a.seq="ACGT" * 200 + a.flag = 0 + a.rname = 0 + a.pos = 33 + a.mapq = 20 + a.cigar = ( (0,10), (2,1), (0,25) ) + a.mrnm = 0 + a.mpos=200 + a.isize=167 + a.qual="1234" * 200 + + return a + +class TestDeNovoConstruction(unittest.TestCase): + '''check BAM/SAM file construction using ex3.sam + + (note these are +1 coordinates): + + read_28833_29006_6945 99 chr1 33 20 10M1D25M = 200 167 AGCTTAGCTAGCTACCTATATCTTGGTCTTGGCCG <<<<<<<<<<<<<<<<<<<<<:<9/,&,22;;<<< NM:i:1 RG:Z:L1 + read_28701_28881_323b 147 chr2 88 30 35M = 500 412 ACCTATATCTTGGCCTTGGCCGATGCGGCCTTGCA <<<<<;<<<<7;:<<<6;<<<<<<<<<<<<7<<<< MF:i:18 RG:Z:L2 + ''' + + header = { 'HD': {'VN': '1.0'}, + 'SQ': [{'LN': 1575, 'SN': 'chr1'}, + {'LN': 1584, 'SN': 'chr2'}], } + + bamfile = "ex6.bam" + samfile = "ex6.sam" + + def checkFieldEqual( self, read1, read2, exclude = []): + '''check if two reads are equal by comparing each field.''' + + for x in ("qname", "seq", "flag", + "rname", "pos", "mapq", "cigar", + "mrnm", "mpos", "isize", "qual", + "bin", + "is_paired", "is_proper_pair", + "is_unmapped", "mate_is_unmapped", + "is_reverse", "mate_is_reverse", + "is_read1", "is_read2", + "is_secondary", "is_qcfail", + "is_duplicate"): + if x in exclude: continue + self.assertEqual( getattr(read1, x), getattr(read2,x), "attribute mismatch for %s: %s != %s" % + (x, getattr(read1, x), getattr(read2,x))) + + def setUp( self ): + + + a = pysam.AlignedRead() + a.qname = "read_28833_29006_6945" + a.seq="AGCTTAGCTAGCTACCTATATCTTGGTCTTGGCCG" + a.flag = 99 + a.rname = 0 + a.pos = 32 + a.mapq = 20 + a.cigar = ( (0,10), (2,1), (0,25) ) + a.mrnm = 0 + a.mpos=199 + a.isize=167 + a.qual="<<<<<<<<<<<<<<<<<<<<<:<9/,&,22;;<<<" + a.tags = ( ("NM", 1), + ("RG", "L1") ) + + b = pysam.AlignedRead() + b.qname = "read_28701_28881_323b" + b.seq="ACCTATATCTTGGCCTTGGCCGATGCGGCCTTGCA" + b.flag = 147 + b.rname = 1 + b.pos = 87 + b.mapq = 30 + b.cigar = ( (0,35), ) + b.mrnm = 1 + b.mpos=499 + b.isize=412 + b.qual="<<<<<;<<<<7;:<<<6;<<<<<<<<<<<<7<<<<" + b.tags = ( ("MF", 18), + ("RG", "L2") ) + + self.reads = (a,b) + + def testSAMWholeFile( self ): + + tmpfilename = "tmp_%i.sam" % id(self) + + outfile = pysam.Samfile( tmpfilename, "wh", header = self.header ) + + for x in self.reads: outfile.write( x ) + outfile.close() + + self.assertTrue( checkBinaryEqual( tmpfilename, self.samfile ), + "mismatch when construction SAM file, see %s %s" % (tmpfilename, self.samfile)) + + os.unlink( tmpfilename ) + + def testBAMPerRead( self ): + '''check if individual reads are binary equal.''' + infile = pysam.Samfile( self.bamfile, "rb") + + others = list(infile) + for denovo, other in zip( others, self.reads): + self.checkFieldEqual( other, denovo ) + self.assertEqual( other, denovo) + + def testSAMPerRead( self ): + '''check if individual reads are binary equal.''' + infile = pysam.Samfile( self.samfile, "r") + + others = list(infile) + for denovo, other in zip( others, self.reads): + self.checkFieldEqual( other, denovo ) + self.assertEqual( other, denovo) + + def testBAMWholeFile( self ): + + tmpfilename = "tmp_%i.bam" % id(self) + + outfile = pysam.Samfile( tmpfilename, "wb", header = self.header ) + + for x in self.reads: outfile.write( x ) + outfile.close() + + self.assertTrue( checkBinaryEqual( tmpfilename, self.bamfile ), + "mismatch when construction BAM file, see %s %s" % (tmpfilename, self.bamfile)) + + os.unlink( tmpfilename ) + + # TODOS # 1. finish testing all properties within pileup objects # 2. check exceptions and bad input problems (missing files, optional fields that aren't present, etc...) if __name__ == "__main__": + # build data files + print "building data files" + subprocess.call( "make", shell=True) + print "starting tests" unittest.main()