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
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
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
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)
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
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)
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":
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)
# 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)
# 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
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.'''
a = makeAlignedRead( alignment )
(<object>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.
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()
+
######################################################################
######################################################################
######################################################################
VALID_HEADER_TYPES = { "HD" : dict,
"SQ" : list,
"RG" : list,
- "PG" : dict,
+ "PG" : list,
"CO" : list }
# order of records within sam headers
'''*(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``.
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
# 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 = <bam1_t*>calloc(1, sizeof(bam1_t))
+
def _isOpen( self ):
'''return true if samfile has been opened.'''
return self.samfile != NULL
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.
# 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'
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
# 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:
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 )
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.
cdef int rtid
cdef int rstart
cdef int rend
+ cdef int max_pos
+ max_pos = 2 << 29
rtid = rstart = rend = 0
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,
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, <void*>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, <void*>callback, fetch_callback )
else:
if region:
return IteratorRow( self, rtid, rstart, rend )
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
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, <void*>callback )
+ buf = bam_plbuf_init( <bam_pileup_f>pileup_callback, <void*>callback )
bam_fetch(self.samfile.x.bam,
self.index, rtid, rstart, rend,
buf, pileup_fetch_callback )
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 )
# 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)
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
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
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
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 = <bam1_t*>calloc( sizeof( bam1_t), 1 )
+ # see bam_init1
+ self._delegate = <bam1_t*>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 = <uint8_t *>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):
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 = <unsigned char*>(&t.core)
+ # oo = <unsigned char*>(&o.core)
+ # for x from 0 <= x < sizeof( bam1_core_t): print x, tt[x], oo[x]
+ # tt = <unsigned char*>(t.data)
+ # oo = <unsigned char*>(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 <char *>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,
+ <uint8_t*>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] = <uint8_t>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 = <char*>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 = <int>bam_aux2i(s)
+ s += 2
+ elif tpe == 'I'[0]:
+ value = <int>bam_aux2i(s)
+ s += 4
+ elif tpe == 'F'[0]:
+ value = <float>bam_aux2f(s)
+ s += 4
+ elif tpe == 'D'[0]:
+ value = <double>bam_aux2d(s)
+ s += 8
+ elif tpe == 'C'[0]:
+ value = <int>bam_aux2i(s)
+ s += 1
+ elif tpe == 'A'[0]:
+ # there might a more efficient way
+ # to convert a char into a string
+ value = "%c" % <char>bam_aux2A(s)
+ s += 1
+ elif tpe == 'Z'[0]:
+ value = <char*>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 = "<cccf"
+ elif t == types.IntType:
+ if value < 0:
+ if value >= -127: fmt, pytype = "<cccb", 'c'
+ elif value >= -32767: fmt, pytype = "<ccch", 's'
+ elif value < -2147483648: raise ValueError( "integer %i out of range of BAM/SAM specification" % value )
+ else: fmt, ctype = "<ccci", 'i'[0]
+ else:
+ if value <= 255: fmt, pytype = "<cccB", 'C'
+ elif value <= 65535: fmt, pytype = "<cccH", 'S'
+ elif value > 4294967295: raise ValueError( "integer %i out of range of BAM/SAM specification" % value )
+ else: fmt, pytype = "<cccI", 'I'
+ else:
+ # Note: hex strings (H) are not supported yet
+ if len(value) == 1:
+ fmt, pytype = "<cccc", 'A'
+ else:
+ fmt, pytype = "<ccc%is" % (len(value)+1), 'Z'
+
+ size = struct.calcsize(fmt)
+ if offset + size > 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 <int>bam_aux2i(v)
elif type == 'f':
- return bam_aux2f(v)
+ return <float>bam_aux2f(v)
elif type == 'd':
- return bam_aux2d(v)
+ return <double>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' % <char>bam_aux2A(v)
elif type == 'Z':
- return bam_aux2Z(v)
+ return <char*>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.
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.
'''
def __get__(self):
return self._level
-
class Outs:
'''http://mail.python.org/pipermail/python-list/2000-June/038406.html'''
def __init__(self, id = 1):
cdef char ** cargs
cdef int i, n, retval
-
n = len(args)
# allocate two more for first (dummy) argument (contains command)
cargs = <char**>calloc( n+2, sizeof( char *) )
return retval, out_stderr, out_stdout
-__all__ = ["Samfile", "IteratorRow", "IteratorRowAll", "IteratorColumn", "AlignedRead", "PileupColumn", "PileupRead" ]
+__all__ = ["Samfile",
+ "Fastafile",
+ "IteratorRow",
+ "IteratorRowAll",
+ "IteratorColumn",
+ "AlignedRead",
+ "PileupColumn",
+ "PileupProxy",
+ "PileupRead" ]
--- /dev/null
+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
+
+
+
+
+
// 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 {
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).
// 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[] )
{
#endif
#endif
+ extern int optind;
+
// reset getop
optind = 1;
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);
}
// 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);
+}
+
+
+
+
#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);
// 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
#include <stdio.h>
#include <ctype.h>
+#include <errno.h>
#include <assert.h>
#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";
* 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;
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);
}
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;
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));
+}
@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
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;
@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)
/*!
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
*/
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.
*/
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
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
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);
-}
#endif
#include "kstring.h"
#include "bam.h"
+#include "sam_header.h"
#include "kseq.h"
#include "khash.h"
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)
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;
-}
-
#include <math.h>
+#include <assert.h>
#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;
/*
P(<b1,b2>) = \theta \sum_{i=1}^{N-1} 1/i
P(D|<b1,b2>) = \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));
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);
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;
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);
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);
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
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) {
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)
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
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
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) {
}
}
}
- 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);
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;
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);
}
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;
// 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);
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;
fprintf(stderr, "Usage: samtools pileup [options] <in.bam>|<in.sam>\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);
#include <string.h>
#include <stdio.h>
#include <zlib.h>
-#include "bam.h"
+#include <unistd.h>
+#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;
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)
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 <input.srt.bam> <output.bam>\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] <input.srt.bam> <output.bam>\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;
}
#include <math.h>
#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 <in.bam> <out.bam>\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);
}
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)
@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);
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;
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");
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);
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] <out.bam> <in1.bam> <in2.bam> [...]\n\n");
+ fprintf(stderr, "Usage: samtools merge [-nr] [-h inh.sam] <out.bam> <in1.bam> <in2.bam> [...]\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 <out.bam> [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;
}
}
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;
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)
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;
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
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]);
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 <maxMem>] <in.bam> <out.prefix>\n");
+ fprintf(stderr, "Usage: samtools sort [-on] [-m <maxMem>] <in.bam> <out.prefix>\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;
}
#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);
}
#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;
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;
}
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;
}
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;
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) {
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;
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
*/
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
--- /dev/null
+/* The MIT License
+
+ Copyright (c) 2003-2006, 2008, 2009, by Heng Li <lh3lh3@gmail.com>
+
+ 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 <stdlib.h>
+#include <stdio.h>
+#include <string.h>
+#include <stdint.h>
+#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;
+}
--- /dev/null
+/* The MIT License
+
+ Copyright (c) 2003-2006, 2008, 2009 by Heng Li <lh3@live.co.uk>
+
+ 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 <stdint.h>
+
+#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
--- /dev/null
+#ifndef _LH3_KLIST_H
+#define _LH3_KLIST_H
+
+#include <stdlib.h>
+
+#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
#include <ctype.h>
#include <stdlib.h>
#include <string.h>
+#include <errno.h>
#include <unistd.h>
#include <sys/types.h>
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;
}
}
#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()
{
* 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;
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;
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;
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);
return 0;
}
+
/**************************
* HTTP specific routines *
**************************/
}
fp->type = KNF_TYPE_HTTP;
fp->ctrl_fd = fp->fd = -1;
- fp->seek_offset = -1;
+ fp->seek_offset = 0;
return 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;
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;
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) {
* 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;
}
#define netwrite(fd, ptr, len) write(fd, ptr, len)
#define netclose(fd) close(fd)
#else
-#include <winsock.h>
+#include <winsock2.h>
#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)
// 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;
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
#include <unistd.h>
#include "razf.h"
+
#if ZLIB_VERNUM < 0x1221
struct _gz_header_s {
int text;
}
#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;i<v32;i++) rz->index->bin_offsets[i] = byte_swap_8((uint64_t)rz->index->bin_offsets[i]);
for(i=0;i<rz->index->size;i++) rz->index->cell_offsets[i] = byte_swap_4((uint32_t)rz->index->cell_offsets[i]);
#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);
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;
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;
}
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;
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;
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:
}
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;
} 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){
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);
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;
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;
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;
}
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;
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;
}
}
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;
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;
}
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;
#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));
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);
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);
}
#include <stdio.h>
#include "zlib.h"
+#ifdef _USE_KNETFILE
+#include "knetfile.h"
+#endif
+
#if ZLIB_VERNUM < 0x1221
#define _RZ_READONLY
struct _gz_header_s;
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;
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);
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)
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
--- /dev/null
+#include "sam_header.h"
+#include <stdio.h>
+#include <string.h>
+#include <ctype.h>
+#include <stdlib.h>
+#include <stdarg.h>
+
+#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 ( *n<len )
+ {
+ *lineptr = realloc(*lineptr, len);
+ *n = len;
+ }
+ if ( !*lineptr ) {
+ debug("[nextline] Insufficient memory!\n");
+ return 0;
+ }
+
+ memcpy(*lineptr,text,len);
+ (*lineptr)[len-1] = 0;
+
+ return to;
+}
+
+// name points to "XY", value_from points to the first character of the value string and
+// value_to points to the last character of the value string.
+static HeaderTag *new_tag(const char *name, const char *value_from, const char *value_to)
+{
+ HeaderTag *tag = malloc(sizeof(HeaderTag));
+ int len = value_to-value_from+1;
+
+ tag->key[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; idict<n; idict++)
+ {
+ const list_t *tmpl_hlines = dicts[idict];
+
+ while ( tmpl_hlines )
+ {
+ list_t *out_hlines = out_dict;
+ int inserted = 0;
+ while ( out_hlines )
+ {
+ status = sam_header_compare_lines(tmpl_hlines->data, 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;
+}
+
+
--- /dev/null
+#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
#include <string.h>
#include <stdio.h>
#include <unistd.h>
+#include <math.h>
+#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()
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;
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
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" )
'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} }
-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 < $< > $@
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
@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<<<<
+
--- /dev/null
+@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<<<<
--- /dev/null
+@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
--- /dev/null
+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
+import sys
import pysam
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##################"
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 )
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"
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"
#!/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
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 ):
except OSError, e:
print >>sys.stderr, "Execution failed:", e
+
class BinaryTest(unittest.TestCase):
'''test samtools command line commands and compare
against pysam commands.
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):
def tearDown(self):
self.samfile.close()
-class TestAlignedReadBam(unittest.TestCase):
+class TestAlignedReadFromBam(unittest.TestCase):
def setUp(self):
self.samfile=pysam.Samfile( "ex3.bam","rb" )
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):
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" )
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()
projects / python-pysam.git / commitdiff